Characterization of the Bacterial Microbiome in Natural Populations of Barley Stem Gall Midge, Mayetiola hordei, in Morocco

Simple Summary Barley stem gall midge, Mayetiola hordei (Kieffer), is one of the most devastating barley pests in many regions of the world, causing significant losses to agricultural output. In this work, we explore the presence of reproductive symbionts, including Wolbachia, Spiroplasma, Cardinium, and Arsenophonus, which can be employed as a substitute for chemical pesticides in pest management, in addition to the impact of their geographic origin on bacterial composition, as well as the diversity of the natural populations of barley stem gall midges detected in four Moroccan barley-growing regions. The screening revealed the presence of Wolbachia; in contrast, none of the samples analyzed had Spiroplasma, Cardinium, or Arsenophonus infections. Overall, 5 phyla, 7 classes, and 42 genera were detected throughout all samples, with a significant environmental influence on the taxonomic assortment of the microbiota. Abstract Mayetiola hordei (Kieffer), known as barley stem gall midge, is one of the most destructive barley pests in many areas around the world, inflicting significant qualitative and quantitative damage to crop production. In this study, we investigate the presence of reproductive symbionts, the effect of geographical origin on the bacterial microbiome’s structure, and the diversity associated with natural populations of M. hordei located in four barley-producing areas in Morocco. Wolbachia infection was discovered in 9% of the natural populations using a precise 16S rDNA PCR assay. High-throughput sequencing of the V3-V4 region of the bacterial 16S rRNA gene indicated that the native environments of samples had a substantial environmental impact on the microbiota taxonomic assortment. Briefly, 5 phyla, 7 classes, and 42 genera were identified across all the samples. To our knowledge, this is the first report on the bacterial composition of M. hordei natural populations. The presence of Wolbachia infection may assist in the diagnosis of ideal natural populations, providing a new insight into the employment of Wolbachia in the control of barley midge populations, in the context of the sterile insect technique or other biological control methods.


Introduction
Mayetiola hordei (Kieffer) (Diptera: Cecidomyiidae), the barley stem gall midge (BSGM), is a small fly (2-4 mm in length), and one of the most devastating pests of cereal crops. This insect is native to the Mediterranean regions in North Africa, including Morocco, Algeria, To the best of our knowledge, this is the first study about the BSGM bacterial microbiome, although a few studies have recently been published on the characterization of the bacterial community associated with the Moroccan Hessian fly (Mayetiola destructor), an insect that shares a similar distribution to BSGM, similar morphology, and comparable symptoms concerning host plants [33]. According to Bel Mokhtar et al. (2020), the sequencing of the 16S rRNA bacterial gene (V3-V4 region) revealed various taxa at the phylum level, including Proteobacteria, Bacteroidetes, Actinobacteria, Cyanobacteria, Deinococcus-Thermus, and Firmicutes. At the genus level, different frequencies were observed between collection regions in Morocco for Empedobacter, Ralstonia, Afipia, and Pseudomonas. Regarding Wolbachia-infection, it was detected in all Moroccan regions with a different relative abundance [34]. Additionally, Bansal et al. (2011) identified four phyla in Hessian fly adults based on classical sequencing: Proteobacteria, Firmicutes, Actinobacteria, and Bacteroidetes. At the genus level, various types of bacteria were detected, including Enterobacter, Pantoea, Stenotrophomonas, Pseudomonas, Bacillus, Ochrobactrum, Acinetobacter, Alcaligenes, Nitrosomonas, Arcanobacterium, Microbacterium, Paenibacillus, and Klebsiella [35].
The aim of this study is to screen natural populations of BSGM midge for the presence of reproductive symbionts, Wolbachia, Spiroplasma, Arsenophonus, and Cardinium, using the 16S rRNA gene. In addition, a high-throughput sequencing (HTS) approach is performed using the Illumina-MiSeq technology to investigate the bacterial symbionts associated with natural populations of BSGM larvae to elucidate the relationship between the insects' bacterial community profile and their habitat.

Barley Stem Gall Midge Collection and DNA Isolation
Plants infected with BSGM were obtained from four major barley-producing areas in Morocco: Abda, Rabat-Salé-Zemmour-Zaër, Fès-Meknes, and Doukkala. All infested plant samples were kept at the International Center of Agricultural Research in Dry Areas (ICARDA) in Rabat, Morocco (Table 1), where the BSGM larvae and pupae were collected during the winter of 2020. Larvae were stored in absolute ethanol at −20 • C, whereas pupae protected in infested plant tissues were reared in a growth room under constant conditions (temperature 23 ± 2 • C, relative humidity 65 ± 5%). Once the adults emerged, they were stored in 100% ethanol at −20 • C. Before DNA extraction, each sample was surface sterilized with a 70% v/v ethanol solution, rinsed with sterile deionized water to remove traces of ethanol, left to dry on a sterile surface, and placed in 1.5 mL centrifuge tubes. After surface sterilization, whole-fly DNA was isolated using a modified CTAB (cetyl trimethyl ammonium bromide) method [36]. The quantity and quality of the DNA preparations and the concentration of double-stranded DNA were both analyzed by a Q5000 micro-volume UV-Vis spectrophotometer (Quawell Technology, San Jose, CA, USA). DNA samples were stored in Eppendorf tubes at −20 • C until further analysis.

Screening and Identification of Bacterial Symbionts
From the BSGM samples collected, 152 insects were selected (Table 1), including 108 that were acquired from 4 different areas, whereas only 2 regions (Rabat and Fez Meknes) were the source of the 44 adult specimens. However, the larvae from Doukkala and Abda were unable to develop into adults in the laboratory setting. The detection of Wolbachia, Spiroplasma, Cardinium, and Arsenophonus was performed by PCR assays using 16S rDNA genus specific primers. Blank and negative controls were incorporated during DNA extractions, and the PCRs were conducted using identical conditions. However, no amplicons were produced from these samples. The quality of the extracted DNA was tested by amplifying part of the mitochondrial 12S rRNA gene (Table S2). Amplification was performed in 25 µL reaction mixtures containing 2.5 µL KAPA Taq buffer (KAPA BioSystems, Wilmington, MA, USA) 10×, 0.25 µL dNTPs (25 mM), 0.25 µL of KAPA Taq DNA polymerase, 0.5 µL of the forward primer (25 µM), 0.5 µL of the reverse primer (25 µM), 1 µL of template DNA solution, and was finalized with 20 µL sterile deionized water. The PCR temperature profile was 95 • C for 5 min followed by 35 cycles of 95 • C for 30 s, 30 s at the optimum annealing temperature for each pair of primers, 1 min at 72 • C, and a final extension step of 72 • C for 5 min. PCR products were electrophoresed on a 1.5% agarose gel in order to examine the presence and size of the amplified fragments. The primer sequences used in this study along with the product size and annealing temperature are summarized in Table S2. The PCR-positive products were purified using polyethylene glycol (20% PEG, 2.5 M NaCl) [37] and resuspended in 15 µL of water. Sanger sequencing was performed on purified products with the BigDye Terminator v3.1 Cycle Sequencing Kit following the manufacturer's recommendations (Applied Biosystems, Waltham, MA, USA). Reaction products were purified using an ethanol/EDTA protocol according to the manufacturer's instructions (Applied Biosystems, Waltham, MA, USA) and were analyzed in an ABI PRISM 3500 Genetic Analyzer (Applied Biosystems, Waltham, MA, USA).

Amplification of the V3-V4 Region, Index PCR, Purification, and Illumina Sequencing
For the high-throughput amplicon sequencing analysis, 78 samples of larvae were selected. However, 3 of these larvae (1 from Abda, 1 from Doukkala, and 1 from Rabat) were disregarded because they had low read counts (less than 1000) to prevent issues related to infinite values that may arise when applying the logarithmic transformation. Therefore, the final number of larvae analyzed was 75, including 19 larvae from Rabat, 17 from Abda, 19 from Doukkala, and 20 from Fes-Meknes (Table 1). The hypervariable V3-V4 region of the bacterial 16S rRNA gene was amplified using MiSeq universal primers 341F and 805R (Table S2). The PCR amplification was performed in two steps. The first amplification was performed using the bacterial universal primers 27F-1492R (Table S2)   From all samples examined, the indexed amplicons were mixed in an equimolar ratio (8 nM) and the sequencing was performed by Macrogen using a 2 × 300 bp pair-end kit on a MiSeq platform [38]. Raw sequencing reads were demultiplexed, converted to FASTQ, and the Illumina adapters were trimmed using Illumina standard algorithms.

NGS Sequencing and Bioinformatics Analysis
Bioinformatics analysis of raw sequencing reads was performed using USEARCH v. 11 and QIIME2 distribution [39,40]. Paired-end reads were assembled and trimmed by length using the usearch -fastq_mergepairs command. The quality of assembled sequences was improved using -fastq_filter, followed by the -fastx_uniques command to detect unique read sequences and their frequencies. Sequences were clustered into operational taxonomic units (OTUs) using the -cluster_otus command at 90% OTU clustering based on the UPARSE algorithm [41]. Crosstalk errors were identified and filtered using the -uncross command based on the UNCROSS2 algorithm [42]. The taxonomy was assigned with Qiime2 using the BLAST + algorithm against the SILVA 138 release database [43,44]. The diversity of individual samples was identified by calculating alpha diversity metrics (evenness, Richness, Shannon, and Simpson indexes) using the "diversity" function of the R package "vegan". Alpha diversity indexes were plotted with the ggplot function from the "ggplot2" package. Statistical differences in bacterial compositions between populations were tested using the nonparametric Kruskal and Wallis and Wilcoxon rank-sum tests [45]. To assess the similarity of bacterial communities in different locations, a beta diversity analysis was performed using the generalized UniFrac distance and visualized by non-metric multidimensional scale (NMDS) maps [46]. Permutational multivariate analysis of variance (PERMANOVA) using distance matrices was performed using the "adonis" function from the R package "vegan" to identify significant differences between isolated groups. A p-value < 0.05 was considered indicative of statistical significance.
A network analysis of OTUs was performed to investigate and visualize the interactions between microorganisms, which could be explained by similar or complementary functions and/or sharing similar environmental conditions, but not necessarily having physical interactions [47,48]. Networks were obtained using the CoNet plugin in Cytoscape 3.8.2 (Institute for System Biology, Seattle, WA, USA) and visualized with Gephi 0.9.2 (Institute for System Biology, Seattle, WA, USA) (Gephi, WebAtlas, Paris, France) [49]. The network was built using Pearson's and Spearman's rank correlation coefficients, mutual information, and the Bray-Curtis and Kullback-Leibler dissimilarity indices. The statistical significance of copresence/mutual exclusion, edge-specific permutation, and bootstrap score distributions were calculated using 1000 iterations. Edges with original scores outside the 0.95 range of their bootstrap distribution were rejected, and p-values were corrected using the Benjamini-Hochberg method.

Phylogenetic Analysis
Phylogenetic analysis was based on partial 16S rRNA gene sequences obtained from Wolbachia-infected specimens and sequences representing Wolbachia-associated OTUs obtained by Illumina MiSeq sequencing. Multiple alignments were performed with MUSCLE, implemented in MEGA 7 software, using standard algorithm parameters [50,51]. Alignments were manually edited and trimmed to adjust the sequence length. Phylogenetic tree reconstruction was based on the maximum likelihood statistical model and was performed using MEGA 7 software. The (GTR + G + I) substitution model was used to estimate nucleotide evolution [52,53]. A bootstrap test with 1000 iterations was used to assess the reliability of the produced phylogeny [54]. The entirety of the 16S rRNA gene sequences produced in this work were submitted to the GenBank database of NCBI with accession numbers OQ189903-OQ189915.

Core Microbiome
Core microbiome analyses were performed using microbiome R package based on the M. hordei dataset [55]. The analysis was performed to detect major bacterial taxa that were recognized in 75% of the samples with a relative abundance threshold value greater than 0.01%, as well as to determine the taxa that were common between the different population regions, and to differentiate them from those which were specific to some regions.

Infection Prevalence
PCR screening was used to assay the presence of four reproductive symbionts (Wolbachia, Spiroplasma, Cardinium, and Arsenophonus) in four natural BSGM populations. In total, 152 BSGMs were screened for reproductive symbionts (Table S1). The screening results indicate that the barley midge is infected only by Wolbachia with a prevalence of 9%. Interestingly, the Wolbachia percentage infection in natural populations did not reveal an even distribution among the different locations. In fact, only BSGMs from Rabat and Fes-Meknes were found to be infected with Wolbachia. In total, 13 flies were infected, five larvae, two females, and three males, out of 44 samples examined from Rabat (22%), and one larva, one female, and one male, out of 52 samples examined from Fes-Meknes (6%) (Table S1, Figure 1). On the contrary, none of the BSGM populations examined were infected with Spiroplasma, Cardinium, or Arsenophonus.

Phylogenetic Analysis of Wolbachia Sequences in BSGM Populations
The Wolbachia phylogenetic analysis was performed on the thirteen Wolbachia-infected samples based on the partial 16S rRNA gene sequences, using a total of 325 bp (File S1) remaining in the thirteen Wolbachia sequences after manually removing the low-quality bases. The 325 bp still contained enough informative data to establish a distinct phylogenetic signal for Wolbachia found in insects, thus allowing an accurate inference of their evolutionary relationships [56][57][58]. According to the results, the Wolbachia strains detected in BSGM populations belonged to supergroup A, exhibiting a high sequence similarity of pairwise distances (98.5%) with Wolbachia sequences isolated from Glossina and Drosophila species (Figure 1).

Phylogenetic Analysis of Wolbachia Sequences in BSGM Populations
The Wolbachia phylogenetic analysis was performed on the thirteen Wolbachia-infected samples based on the partial 16S rRNA gene sequences, using a total of 325 bp (File S1) remaining in the thirteen Wolbachia sequences after manually removing the low-quality bases. The 325 bp still contained enough informative data to establish a distinct phylogenetic signal for Wolbachia found in insects, thus allowing an accurate inference of their evolutionary relationships [56][57][58]. According to the results, the Wolbachia strains detected in BSGM populations belonged to supergroup A, exhibiting a high sequence similarity of pairwise distances (98.5%) with Wolbachia sequences isolated from Glossina and Drosophila species (Figure 1).

16S rRNA Amplicon Sequencing
Applying Illumina high-throughput sequencing of 16S rRNA gene amplicons, the bacterial community composition and diversity of 75 natural barley stem gall midge samples from Rabat, Doukkala, Fes-Meknes, and Abda locations were examined. After sequencing and quality filtering, a total of 1,757,130 qualified paired-end reads with an average count per sample of 44,504 reads were divided into 638 OTUs. Based on a 90% sequence similarity, 61 OTUs were classified in five phyla, that is, for Proteobacteria, Firmicutes, Dependentiae (a recently characterized phylum, Table S3), Bacteroidetes, and Actinobacteria. Seven classes, Gammaproteobacteria, Alphaproteobacteria, Bacilli, Bacterodia, Actinobacteria, Negativicutes (a new class, Table S4), and Babeliae, plus 42 genera were procured across all the samples (Table S5).

Bacterial Diversity among BSGM Natural Populations
The bacterial communities of BSGM were clustered according to their geographic origins. The NMDS plots based on the generalized UniFrac distance revealed that samples from Rabat, Fes-Meknes, Abda, and Doukkala were significantly separated (PERMANOVA, p < 0.001, Figure 2), with similar results also obtained with pairwise analyses (Supplementary Figure S1).

16S rRNA Amplicon Sequencing
Applying Illumina high-throughput sequencing of 16S rRNA gene amplicon bacterial community composition and diversity of 75 natural barley stem gall midge ples from Rabat, Doukkala, Fes-Meknes, and Abda locations were examined. Aft quencing and quality filtering, a total of 1,757,130 qualified paired-end reads with a erage count per sample of 44,504 reads were divided into 638 OTUs. Based on a 90 quence similarity, 61 OTUs were classified in five phyla, that is, for Proteobacteria, Fi utes, Dependentiae (a recently characterized phylum, Table S3), Bacteroidetes, an tinobacteria. Seven classes, Gammaproteobacteria, Alphaproteobacteria, Bacilli, Ba dia, Actinobacteria, Negativicutes (a new class, Table S4), and Babeliae, plus 42 g were procured across all the samples (Table S5).

Bacterial Diversity among BSGM Natural Populations
The bacterial communities of BSGM were clustered according to their geograph igins. The NMDS plots based on the generalized UniFrac distance revealed that sam from Rabat, Fes-Meknes, Abda, and Doukkala were significantly separated ( MANOVA, p < 0.001, Figure 2), with similar results also obtained with pairwise ana (Supplementary Figure S1).  Based on the number and relative abundance of OTUs and Simpson and Shannon indices ( Figure S2), the four natural populations of BSGM were characterized by different species' richness and diversity. The samples from Abda, Doukkala, and Fes-Meknes contained most of the bacterial species. Rabat samples, on the other hand, had statistically lower species richness and diversity than the other locations. Furthermore, based on the Shannon score, samples from the Doukkala region had statistically higher diversity than samples from the Rabat region (pairwise ANOVA: p < 0.05).

Wolbachia Phylogenetic Analysis
Using conventional sequence clustering into OTUs, the Wolbachia infection status was determined from normalized libraries (90% sequence identity). Wolbachia-related reads were clustered into a single OTU (OTU 156), with prevalence discrepancies between the two regions investigated. In Rabat, seven samples out of 20 were infected, whereas, in Abda, one sample out of 18 was found to be infected with Wolbachia. The sequence of the Wolbachia-related OTU (443 bp) was positioned within the supergroup A sequences according to the phylogenetic analyses based on the 16S rRNA gene (Figure 4). These results indicate that the BSMG derived from different regions most likely carries the same Wolbachia strain.

Wolbachia Phylogenetic Analysis
Using conventional sequence clustering into OTUs, the Wolbachia infection status was determined from normalized libraries (90% sequence identity). Wolbachia-related reads were clustered into a single OTU (OTU 156), with prevalence discrepancies between the two regions investigated. In Rabat, seven samples out of 20 were infected, whereas, in Abda, one sample out of 18 was found to be infected with Wolbachia. The sequence of the Wolbachia-related OTU (443 bp) was positioned within the supergroup A sequences according to the phylogenetic analyses based on the 16S rRNA gene (Figure 4). These results indicate that the BSMG derived from different regions most likely carries the same Wolbachia strain.

Bacterial Co-Occurrence/Mutual Exclusion Network Analysis of BSGM Natural Populations
Co-occurrence and mutual exclusion network analysis was performed to examine the potential interactions between bacterial partners in each region separately. The networks for each region were shown at the family ( Figure 5) and genus ( Figure S3

Bacterial Co-Occurrence/Mutual Exclusion Network Analysis of BSGM Natural Populations
Co-occurrence and mutual exclusion network analysis was performed to examine the potential interactions between bacterial partners in each region separately. The networks for each region were shown at the family ( Figure 5) and genus ( Figure S3

Core Microbiome Analysis of BSGM Populations
The core bacterial community was composed of 20 OTUs out of 61 (33%  (Figure 6) (see Table S7 for more details about the overall prevalence of the core microbiome).
[OTU33]) ( Figure 6) (see Table S7 for more details about the overall prevalence of the core microbiome).

Discussion
In the current study, we evaluated the prevalence of reproductive symbionts, the bacterial diversity of natural populations, and the impact of geographic origin on the bacterial community structure of barley stem gall midges, located in four different producing regions in Morocco. The bacterial microbiome analysis was performed using next-generation, high-throughput sequencing of the V3-V4 region of the bacterial 16S rRNA gene.
To analyze the existence of reproductive symbionts, a genus specific 16S rRNA PCR test was used. The screening results revealed that BSGM was only infected with Wolbachia. Surprisingly, Wolbachia infection was only found in BSGM populations from Rabat and Fez-Meknes, five larvae, two females, and three males from 44 samples from Rabat (22%), and one larva, one female, and one male from 52 individuals from Fez-Meknes (6%); in total 13 samples out of 152 BSGMs were infected, with a prevalence of 9%. However, none of the BSGM populations studied contained Spiroplasma, Cardinium, or Arsenophonus. The HTS DNA sequence datasets also revealed the presence of Wolbachia infections, which were detected in two out of four regions studied. Seven samples out of 20 were infected in Rabat (35%), while in Abda, one sample out of 18 (5%) was found to be infected with Wolbachia (Table S6). The discrepancy in prevalence between the two methods may be attributable to the low Wolbachia infection density, as well as to the exclusion of some samples from the HTS sequence datasets due to their low read counts (less than 1000) and the higher sensitivity of HTS compared to conventional PCR screening [54,59]. To the best of our knowledge, Wolbachia had never been identified in BSGM populations. However, it was detected in a recent study of Hessian fly originating from Morocco, with a total of 13 samples out of 244 (5%) carrying Wolbachia infections: three females and ten males. Concerning HTS DNA sequence datasets, Wolbachia infections were also detected: 31 individuals out of 40 (78%) in various Hessian fly populations [34]. The Hessian fly shares with BSGM a similar distribution, similar morphology, and comparable symptoms in host plants [29]. On the other hand, no Wolbachia was found in the few other Hessian fly studies

Discussion
In the current study, we evaluated the prevalence of reproductive symbionts, the bacterial diversity of natural populations, and the impact of geographic origin on the bacterial community structure of barley stem gall midges, located in four different producing regions in Morocco. The bacterial microbiome analysis was performed using next-generation, high-throughput sequencing of the V3-V4 region of the bacterial 16S rRNA gene.
To analyze the existence of reproductive symbionts, a genus specific 16S rRNA PCR test was used. The screening results revealed that BSGM was only infected with Wolbachia. Surprisingly, Wolbachia infection was only found in BSGM populations from Rabat and Fez-Meknes, five larvae, two females, and three males from 44 samples from Rabat (22%), and one larva, one female, and one male from 52 individuals from Fez-Meknes (6%); in total 13 samples out of 152 BSGMs were infected, with a prevalence of 9%. However, none of the BSGM populations studied contained Spiroplasma, Cardinium, or Arsenophonus. The HTS DNA sequence datasets also revealed the presence of Wolbachia infections, which were detected in two out of four regions studied. Seven samples out of 20 were infected in Rabat (35%), while in Abda, one sample out of 18 (5%) was found to be infected with Wolbachia (Table S6). The discrepancy in prevalence between the two methods may be attributable to the low Wolbachia infection density, as well as to the exclusion of some samples from the HTS sequence datasets due to their low read counts (less than 1000) and the higher sensitivity of HTS compared to conventional PCR screening [54,59]. To the best of our knowledge, Wolbachia had never been identified in BSGM populations. However, it was detected in a recent study of Hessian fly originating from Morocco, with a total of 13 samples out of 244 (5%) carrying Wolbachia infections: three females and ten males. Concerning HTS DNA sequence datasets, Wolbachia infections were also detected: 31 individuals out of 40 (78%) in various Hessian fly populations [34]. The Hessian fly shares with BSGM a similar distribution, similar morphology, and comparable symptoms in host plants [29]. On the other hand, no Wolbachia was found in the few other Hessian fly studies [35,60,61]. Wolbachia was also reported in the Asian rice gall midge (Orseolia oryzae) in India, where the sequencing of the 16S rRNA bacterial gene (V3-V4 region) revealed that Wolbachia was the predominant bacterium in the pupae and adults of both genders, with prevalences of 97.3%, 89.8%, and 79.6%, respectively [62]. In addition, Wolbachia is the most widespread intracellular bacterial genus, infecting filarial nematodes and a wide range of arthropod groups, including almost 65% of insect species [63][64][65], such as aphids [36,[66][67][68] and fruit flies (Diptera:Tephritidae) [69][70][71].
Our efforts to identify Wolbachia strains found in natural populations of BSGM were limited by the difficulty in amplifying MLST (multi-locus sequence typing) and wsp genes, which could be attributed to the low infection density of Wolbachia [72,73]. The phylogenetic analyses were limited to 16S rRNA sequences amplified using Wolbachia-specific primers and the Wolbachia-related sequences acquired using HTS. All Wolbachia sequences found in this investigation exhibited the highest homology to strains from supergroup A (Figures 1 and 6), indicating that the provenance of the BSGM samples did not contribute to the diversification of Wolbachia strains. It was previously discovered that the density of Wolbachia in a host influences the level of CI [74]. In light of this, the low infection rate found in our study implies that Wolbachia is unlikely to cause CI in the studied BSGM populations. Additionally, the fact that both males and females carry Wolbachia may indicate that the infection should not have any consequences on the sex ratio in favor of females. Even though Wolbachia are typically transmitted vertically, the population of barley midges in Rabat, Fes-Meknes, and Abda most likely became infected through horizontal transmission, which is common among insects [75,76] and may be mediated by parasitoids or host plants [77,78]. However, our results provide a fresh insight into the BSGM capacity to harbor the Wolbachia infection that may allow us to use the transinfection method to establish a reliable Wolbachia-infected lineage, by directing symbiont strains to be transferred to new hosts either within the same species or between species. Since it has been previously accomplished, the case of Wolbachia transinfections includes transporting strains between Drosophila species [79,80], Aedes aegypti [81][82][83], Culex quinquefasciatus [84], Bactrocera oleae [85], and Ceratitis capitata [86][87][88][89].
In order to investigate the effect of geographical origin on bacterial structure, HT sequencing was applied on 75 natural barley stem gall midge samples, based on a 90% sequence similarity. The results reveal a high number of distinct OTUs from four separate BSGM populations (Table S5 and Figure 3c). Our findings demonstrate that native environments have a major impact on the microbiota. The four BSGM natural populations were distinguished by various levels of species richness and diversity. Most bacterial species were found in the samples from Abda, Doukkala, and Fes-Meknes. In contrast, Rabat samples had statistically lower species richness and diversity. Additionally, samples from the Doukkala region displayed significantly higher diversity than samples from the Rabat region, according to the Shannon index. Since the research on the bacterial microbiome of this insect is still limited, it is the first time that changes in the structure of microbial communities have been observed in populations of BSGM with various geographic distributions. However, it has also been reported in other fly or insect species, including cases of the Hessian fly (M. destructor) [34,61], oriental fruit fly (Bactrocera dorsalis) [90,91], and melon fly (Zeugodacus cucurbitae) [92,93].
Proteobacteria was the predominant phylum detected in the four regions of Abda, Doukkala, Fes-Meknes, and Rabat. This phylum was also the most abundant in other insects, such as the Hessian fly from Morocco (three out of four regions) [34], the Asian rice gall midge (O. oryzae) from India [63], and the melon fly from Bangladesh (Z. cucurbitae) [92,93]. However, the classification of the rest of the phyla, Actinobacteria, Bacteroidetes, and Firmicutes, was different from one region to the other, in agreement with the results for other insects [34,62,92]. Additionally, a recently characterized phylum was detected in our study as the fifth most-common phylum among the four studied regions, with a prevalence of 32% (24 out of 75 samples), that is, candidate phylum Dependentiae (also known as TM6). The identification of this phylogenetic group of bacteria was based on isolates infecting free-living amoebae and on metagenomic studies. This phylum is common in a variety of habitats, including hospital biofilms, soil, waste water, and peat bogs; it is characterized by a relatively limited metabolic potential, but encodes a vast array of transporters, including ATP transporters and genes with bacterial endosymbiontenriched activities [94][95][96][97][98]. At the class level, Gammaproteobacteria and Alphaproteobacteria were the most prevalent members in all samples analyzed, followed by Bacilli in Abda, Doukkala, and Fes-Meknes, and Actinobacteria in Rabat. However, in the Hessian fly populations from Morocco, the most frequent classes were Gammaproteobacteria, Alphaproteobacteria, and Betaproteobacteria [34]. Interestingly, members of Negativicutes were also detected but at a much lesser degree (Figure 3b), with a prevalence of 9.33% (seven out of 75 samples). This class was recently certified in the phylum Firmicutes that has generally been defined by having a low genomic GC content, Gram-positive (monoderm) cells, and an envelope composed of a cytoplasmic membrane and a thick coating of peptidoglycan [99]. However, Negativicutes are diderms, characterized by an inner and outer membrane containing lipopolysaccharides [100,101]. In addition, Negativicutes have colonized a variety of hosts, including insects, animals, soil, ocean, and sediments. They mainly rely on fermentation because they cannot grow in the presence of oxygen. Some members of Negativicutes are known by their ability to ferment lactate, as in the case of Veillonella, the genus that was present in our results, while other Gram-negative members can produce endospores [102,103]. At the genus level, Pseudomonas exhibited the highest relative abundances in Abda, Doukkala, and Rabat; this bacterial genus is known from a wide range of habitats and hosts and exhibits remarkable metabolic diversity [104,105]. The genus has adapted to interact specifically in both useful and harmful ways, mostly with plants [106,107], humans [108], and insects [109,110]. In the same three regions, Pseudomonas was followed by Stenotrophomonas, a genus reported from a variety of habitats that is frequently characterized as a multidrug-resistant, opportunistic human pathogen [111]. However, they are more usually discovered in soils or in close proximity to plants [112], where they can create symbiotic relationships [113]. Stenotrophomonas is also linked to a variety of insect species [114,115]. On the other hand, Pantoea was prevalent in the region of Fes-Meknes; this genus has regularly been isolated from a variety of terrestrial and aquatic ecosystems, in addition to connections with humans, animals, and insects [116,117]. Additionally, Pantoea isolates have been used as immunopotentiators to create supportive medications for melanoma, infections, allergies, and the reversal of immunosuppression [118][119][120]. Furthermore, several Pantoea isolates have been reported to produce antimicrobials and have been converted into commercial biocontrol treatments [121][122][123]. Pantoea was followed by Pseudomonas in the region of Fes-Meknes; this genus was also present in Hessian flies with different frequencies among different populations [34]. The reported variation in the microbiome of BSGM populations could be influenced by barley cultivars, insect pressure, environmental factors, including annual mean temperature, precipitation, and maximum snowfall, as well as the effects of host plants and insect developmental stages and sexes, plus the effect of Wolbachia infection on the microbiome [124,125].

Conclusions
The bacterial composition of populations of BSGM was found to be significantly influenced by the geographic location. With the highest number of OTUs, Proteobacteria were found to be the most abundant phylum in flies from the four sites. To the best of our knowledge, this is the first study of the BSGM bacterial microbiome that reports Wolbachia infection in this insect. The low infection rate identified in our study led us to suspect that Wolbachia is unlikely to be the source of cytoplasmic incompatibility in BSGM reproduction. On the other hand, our findings may offer a new perspective on the application of Wolbachia in the control of BSGM.

Supplementary Materials:
The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/microorganisms11030797/s1, Table S1. Prevalence of bacterial endosymbionts screened in populations of barley stem gall midges. Table S2. List of bacterial primers and annealing temperatures used; Table S3. Distribution of Dependencia family (Otu 45); Table S4. Distribution of Negativicutes class (Otu 106); Table S5. Relative abundance (%) of different bacteria detected in Barley stem gall midge regions of Abda, Doukkala, Fes Meknes, and Rabat; File S1. Wolbachia sequence alignments; Figure S1. Pairwise non-metric multidimensional scaling (NMDS) plot of bacterial communities in barley stem gall midge samples collected from Rabat (purple), Fes-Meknes (blue), Abda (red), and Doukkala (green) (p < 0.001). The 'd' stands for the grid's dissimilarity scale (d = 0.2 implies that the distance between two grid lines represents approximately 20% dissimilarity between the regions); Figure S2. Species richness and diversity indices with significance differences in natural populations of barley stem gall midge samples collected from the regions of Rabat, Fes Meknes, Abda, and Doukkala. The interquartile range (IQR) is represented by boxes, the median is represented by a line within the boxes, and samples are represented by dots (* 0.01 < p < 0.05, ** p ≤ 0.01); Figure S3. Co-occurrence and mutual exclusion networks at the genus level for OTUs that compose the bacterial communities of Mayetiola hordei natural populations from the four regions of (a) Fes-Meknes, (b) Rabat, (c) Doukkala, and (d) Abda. The degree of interaction determines the size of each node. Cases of copresence are represented by green edges, while mutual exclusion is represented by red edges. Table S6. Distribution of Wolbachia genus (Otu 156). Table S7. Overall core prevalence. Funding: This research was conducted with the assistance of the European Union's Erasmus + Programme. Additionally, this study has been co-financed by the European Regional Development Fund of the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship and Innovation, under the call RESEARCH-CREATE-INNOVATE (project code: T2EDK-2020).

Data Availability Statement:
All the data that support the findings of this study are available in NCBI under BioProject PRJNA915357.