key: cord-0053703-r8vqu9mu
authors: Bogale, Haikel N.; Cannon, Matthew V.; Keita, Kalil; Camara, Denka; Barry, Yaya; Keita, Moussa; Coulibaly, Drissa; Kone, Abdoulaye K.; Doumbo, Ogobara K.; Thera, Mahamadou A.; Plowe, Christopher V.; Travassos, Mark; Irish, Seth; Serre, David
title: Relative contributions of various endogenous and exogenous factors to the mosquito microbiota
date: 2020-12-10
journal: Parasit Vectors
DOI: 10.1186/s13071-020-04491-7
sha: 1d48fcf0ca9ed4740c08a771be0582f708ba99e6
doc_id: 53703
cord_uid: r8vqu9mu

BACKGROUND: The commensal microbiota of mosquitoes impacts their development, immunity, and competency, and could provide a target for alternative entomological control approaches. However, despite the importance of the mosquito/microbiota interactions, little is known about the relative contribution of endogenous and exogenous factors in shaping the bacterial communities of mosquitoes. METHODS: We used a high-throughput sequencing-based assay to characterize the bacterial composition and diversity of 665 individual field-caught mosquitoes, as well as their species, genotype at an insecticide resistance locus, blood-meal composition, and the eukaryotic parasites and viruses they carry. We then used these data to rigorously estimate the individual effect of each parameter on the bacterial diversity as well as the relative contribution of each parameter to the microbial composition. RESULTS: Overall, multivariate analyses did not reveal any significant contribution of the mosquito species, insecticide resistance, or blood meal to the bacterial composition of the mosquitoes surveyed, and infection with parasites and viruses only contributed very marginally. The main driver of the bacterial diversity was the location at which each mosquito was collected, which explained roughly 20% of the variance observed. CONCLUSIONS: This analysis shows that when confounding factors are taken into account, the site at which the mosquitoes are collected is the main driver of the bacterial diversity of wild-caught mosquitoes, although further studies will be needed to determine which specific components of the local environment affect bacterial composition. GRAPHIC ABSTRACT: [Image: see text]

Aedes, Culex, and Anopheles mosquitoes can transmit eukaryotic parasites and viruses to humans and are responsible for devastating diseases, such as malaria, lymphatic filariasis, dengue, West Nile, and Zika, that affect hundreds of millions of people worldwide and cause hundreds of thousands of deaths annually [1] . Successful disease elimination campaigns have focused on interrupting transmission of these diseases by targeting their mosquito vectors. Thus, in the early twentieth century, vector control approaches primarily relied on the use of chemicals and larviciding tools (i.e., petroleum oils and larvivorous fish) to eliminate larval and adult stages of the mosquito [2, 3] . Environmental and human health concerns brought by the persistent use of chemicals like

Parasites & Vectors *Correspondence: dserre@som.umaryland.edu 1 Institute for Genome Sciences, University of Maryland School of Medicine, Baltimore, MD, USA Full list of author information is available at the end of the article dichlorodiphenyltrichloroethane (DDT) over the years have partly led to the increased use of pyrethroids, which are safer alternatives, in vector control measures [3] . In the last three decades, entomological control strategies based on pyrethroid-treated bed nets and indoor residual spraying have been extensively used in the fight against vector-borne diseases with considerable success in reducing disease burden [4] [5] [6] . However, wide-spread use of these approaches, combined with exposure to agricultural pest control chemicals, have led to the emergence and rapid spread of insecticide resistance alleles in many areas [7, 8] . In addition, several populations of mosquitoes have modified their behaviors (e.g., their host biting time [9, 10] and location [11, 12] or their host species preference [13, 14] ) upon exposure to insecticides. This acquisition of chemical and behavioral resistance to insecticides threatens the advances made in control of mosquito-borne diseases and highlights the need for alternative measures.

One alternative to chemicals is to leverage biological agents to control mosquito populations (sometimes referred to as biological controls). For example, Bacillus thuringiensis, a spore-forming bacterium with larvicidal characteristics, is used extensively against disease-transmitting insects as well as agricultural pests [15] . In recent years, the use of Wolbachia strains as biological control agents has also gained momentum, with studies revealing an ability of these bacteria to promote pathogen interference and to reduce the life span of mosquitoes [16] [17] [18] . Furthermore, recent findings have highlighted how modifications of the bacterial communities present in the midgut of mosquitoes could decrease or inhibit transmission of pathogens. For example, studies have demonstrated the importance of gut microbiota in mosquito larval development and shown that bacteria are required for Aedes mosquitoes to survive to the adult stages [19, 20] . Similarly, it has been shown that elimination of native microbiota resulted in delayed growth in Anopheles larvae [21] . Functional studies on adult stage mosquitoes showed that the gut microbiota can increase the resistance of mosquitoes to human pathogens by modulating the mosquito innate immune response [22, 23] or directly through production of anti-pathogen molecules from specific microbial species [24, 25] . Overall, these studies demonstrate the potential of microbiota manipulations for inhibiting pathogen transmission and/or reducing vector competence.

However, while these laboratory and field studies highlight the role of the mosquito midgut microbiota in regulating the development and transmission of human pathogens, very little is known about the factors that shape the diversity of the bacterial composition in wild mosquitoes. Some studies have suggested that the location of mosquito collection is associated with the composition of the microbiota [26, 27] , while others have shown that the microbial composition differed between mosquito species, even when they are closely related and collected at the same location [28] , or reared under the same conditions [20] . In addition, it is possible that genetic resistance to insecticides may also influence the microbial composition of field-collected adult mosquitoes (e.g., resistance could result in differential exposure of the gut microbiota to insecticides, which can alter the composition of the microbiota, see, for example, [29] ). Lastly, the source of the blood meal has also been shown to impact gut microbiota composition, with the mammalian blood-meal source altering the gut bacterial composition of adult mosquitoes [30] . However, since these studies typically addressed only one of those factors at a time (without correcting for confounding effects), the relative contribution of each of these parameters to shaping the midgut microbiota composition of wild-caught mosquitoes remains unclear.

We report here the results of our analysis of 665 individual Anopheles mosquitoes collected in Guinea and Mali. We characterized the microbial composition of these mosquitoes and screened them for a large variety of eukaryotic parasites and viruses. We also characterized the species of all mosquitoes, genotyped them at a major associated insecticide resistance locus, and examined the source of their last blood meal. We then tested how these endogenous and exogenous factors influenced the bacterial diversity and simultaneously estimated the relative contribution of these factors to the mosquito microbiota.

Mosquitoes were collected from six sites in Guinea by human landing catches (HLC) and pyrethrum spray catches (PSC) between August and November 2016 and in August 2017 (Table 1 ; Additional file 1: Figure S1 ). The captured mosquitoes were placed in Eppendorf tubes containing ethanol and shipped to the University of Maryland School of Medicine for analysis. Mosquitoes were also collected in Mali from homes in Bandiagara using PSC in July 2011. The captured mosquitoes were placed in Carnoy's solution (volume ratio acetic acid:ethanol, 1:3) and shipped to the molecular biology laboratory of the Malaria Research and Training Center (MRTC) in Bamako, Mali, for DNA extraction. See Table 1 for details on the collection sites and the collection method used.

DNA was extracted from each mosquito using a modified version of the Qiagen 96-well extraction protocol (Qiagen N.V., Hilden, Germany). Briefly, whole individual mosquitoes were randomly placed in each well of a 96-well plate with five 1-mm RNASE free oxide beads for homogenization; 11-14 extraction controls were also included on each plate. Each mosquito was homogenized using a TissueLyser for 6 min at 20 m/s in the lysis buffers provided by the Qiagen DNeasy 96 Blood & Tissue Kit (Qiagen N.V.). The plates were then centrifuged at 1500 g for 3 min. The homogenates were then incubated for 1 h at 55 °C, centrifuged again, and incubated overnight at 55 °C. After a final centrifugation step, the supernatant was transferred to a 96-well DNeasy plate to bind the DNA. The columns were washed twice before elution of the DNA with 100 μl of Qiagen buffer AE (Qiagen N.V.). A NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) was used to determine DNA concentration.

DNA was extracted from the body section of each mosquito using the Chelex-100 chelating resin (Bio-Rad Laboratories, Hercules, CA, USA) protocol. Briefly, a dissecting needle was used to separate out the thorax and abdomen sections of each individual mosquito, each of which was then placed in 1.5-ml tubes containing deionized water. Pipette tips were used to grind each sample in the tube, following which each sample was further homogenized for 20 min in phosphate buffered saline (PBS; 1×)/1% saponin solution, with gentle shaking. The homogenates were then incubated at room temperature (25 °C) overnight, following which the tubes were centrifuged at 20,000 g for 2 min and the supernatants discarded. After the pellets were washed with PBS (1×), each pellet was resuspended in 75 μl of deionized water and 25 μl of 20% Chelex-100 chelating resin solution. This mixture was placed on a heating block for at least 10 min and stirred every 5 min. DNA was transferred into a new tube after a final centrifugation step for 1 min.

DNA extracted from each mosquito was amplified using primers targeting bacterial 16S rRNA primers for microbiota analysis [31] , mosquito knockdown resistance west (kdr-west L1014F mutation variant) for insecticide resistance genotyping [32] , cytochrome c oxidase 1 (COX1) and S200X6.1 [33] loci for mosquito species determination, mammalian mitochondrial 16S rRNA sequences for blood-meal analysis, and eukaryotic parasite and virus primers (targeting 18S rRNA and NS5 loci, respectively) [34] for identification of parasite and virus species (Additional file 2: Table S1 ). The rationale for targeting these specific loci and the primer design are described in detail in the respective references. Briefly, all primers were designed or selected to simultaneously fulfil the following conditions: (i) allow amplification of all sequences of the taxon targeted, (ii) avoid off-target amplification, (iii) provide sufficient sequence information for resolving the genotype or taxonomy, and (iv) be sufficiently short to be sequenced using Illumina technology and with paired end overlaps to allow for sequencing error-correction.

DNA extracted from 665 individual mosquitoes and from the 95 extraction control samples was amplified separately with each primer pair using the Promega GoTaq DNA Polymerase (Promega Corp., Madison, WI, USA) under the following conditions: initial denaturing step at 95 °C for 2 min, followed by 40 cycles of 95 °C for 30 s, 50 °C for 30 s and 72 °C for 30 s, and a final extension of 5 min. Only 35 cycles were used to amplify bacterial 16S rRNA. 

Viral sequences were analyzed by first synthesizing cDNA from carry-over RNA present in DNeasy extracted samples using M-MLV Reverse Transcriptase (Promega Corp.). Briefly, 2 µl of mosquito DNA extract from each sample was incubated with 1 µl of random hexamers (0.5 µg) and 12 µl of RNA-free water at 70 °C for 5 min. After this denaturation step, 1.25 mM of each dNTP, 25 units of recombinant RNasin ® Ribonuclease Inhibitor, and 200 units of M-MLV RT were added and the cDNA synthesis carried out at 37 °C for 60 min. The cDNA products were amplified using the same conditions as described for DNA.

The PCR products generated from one DNA sample (i.e., individual mosquito) were then pooled together: the products from all parasite and viral amplifications were pooled at equal molar concentration (pool1) and the products of the remaining amplifications were pooled together (pool2) in various proportions reflecting the sequence diversity expected at each locus (75% Bacteria_16S rRNA, 5% S200X6, 5% CulicCox1, 10% KDR, 5% mammalian_16S rRNA). Of these, 384 pools (i.e., products from a single mosquito or water control) were then randomly assigned to a well of a 384-well plate and reamplified in a second PCR to add a unique barcode and the Illumina sequencing adaptors [34] . Finally, the barcoded product pools from all mosquitoes were combined (pool1:pool2, 1:3 ratio) and sequenced simultaneously on Illumina HiSeq 2500 sequencing platform (Illumina Inc., San Diego, CA, USA) using a protocol allowing the generation of 300-bp paired-end reads [35] .

For the bioinformatic analyses, we first used the sequences of the barcodes incorporated in the second PCR to assign each read to an individual sample. Next, we used the first 18-27 nucleotides of each read to identify the sequence of the PCR primers used to amplify a given sequence and separated the reads by loci. Further analyses were performed for each locus separately as indicated below.

All reads carrying bacterial 16S rRNA primers (see above) were analyzed in DADA2 [36] (v1.6.0) by first trimming low-quality bases at the end of each read pair (> 250 bp for forward reads, > 210 for reverse reads) using the following parameters: maxN = 0, maxEE = 2, truncQ = 2. Dereplication was done by combining identical reads and assigning the number of reads belonging to each unique read (derepFastq). Next, the dereplicated data were analyzed with the dada core sample inference algorithm followed by the merging of read pairs that overlap by at least 12 bases (mergePairs). An amplicon sequence variant table (ASV) table was constructed with makeSe-quenceTable for all samples, and chimeras were removed using removeBimeraDenovo. For taxonomic assignment, Silva (v128) [37] was used as a training set (using assign-Taxonomy and addSpecies) to create taxonomy data. Finally, the R (v3.4.0) package phyloseq (v1.25.2) [38] was used to combine the ASV table, taxonomy data, and sample metadata for downstream microbiome data analyses. Samples with < 5000 reads were discarded from further analysis, as they likely represent low-level cross contamination. Principal coordinates analysis (PCoA) with Bray-Curtis dissimilarity and weighted UniFrac distance matrices were calculated in phyloseq. Finally, analyses using the Adonis function were performed to simultaneously evaluate the contribution of each factor to the bacterial composition: this multivariate analysis provided a statistical assessment of the association of each factor (i.e., collection site, mosquito species, kdr-w genotype, blood-meal status, and infection status; see following sections for details) with bacterial diversity, as well as an estimate of the proportion of the variance explained.

First, reads that did not contain the exact barcode and primer sequences were discarded and the remaining reads were assigned to a given sample according to their unique barcode sequence. In order to remove low-quality bases and sequences that were likely primer dimers, each sequence was searched for the forward and reverse primers and trimmed after the reverse primer (if both primers are found). Those sequences where the forward primer was found but the reverse primer was missing were left untrimmed. Untrimmed sequences that belonged to a primer with an expected amplicon length of < 300 bp were trimmed from 50 bp from the 5′ end for further quality filtering. Afterwards, filtered pairedend sequences from each read pair were merged using FLASH [39] to generate a consensus sequence of the overlapping region. Each correctly merged sequence was trimmed of both amplification primers (forward and reverse) and kept only if the sequence was longer than 90 bp. Then all concatenated sequences amplified with the same primer pair (from all samples) were compared to each other, and a single copy of each unique sequence was kept (while the number of times it was observed in each sample was recorded). Unique sequences that were observed less than ten times across all samples were removed as they likely represent instances of sequences carrying errors [34, 40] . The remaining unique sequences were compared against all DNA sequences annotated on the NCBI nr database using BLAST [34] . We then retrieved the taxonomic information of the most similar sequence(s) if there was at least 70% identity over the entire sequence length. Finally, we summarized the parasite and virus species identified, the percentage identity (i.e., similarity to the most similar sequence[s] on NCBI), and the number of reads observed in each sample. When evaluating the blood-meal composition of each mosquito, the same procedure as just described was applied, but only samples with at least 1000 reads were considered (to avoid including possible cross-contamination or sequences mis-assigned to one sample due sequencing errors in the barcode) [41] . The distribution of reads generated from the mammalian_16S rRNA primers for all mosquito and negative control samples is presented in Additional file 3: Figure S2 .

Reads amplified from the kdr and S200X6.1 primer pairs were processed as described above, with sequences assigned to their specific sample, filtered for quality, and merged with FLASH.

For the L1014F locus, the top two most abundant unique sequences in each sample (Seq1 and Seq2) were considered for further analysis. If one sample is homozygous, the second most abundant sequence will be a read carrying sequencing error(s) and should account for a small fraction of the Seq1 reads (i.e., Seq2/[Seq1+Seq2] ~ 0). Alternatively, if one sample is heterozygous, we would expect the number of Seq2 reads to be very close to Seq1 and Seq2/(Seq1+Seq2) ~ 0.5. We calculated the ratio Seq2/(Seq1+Seq2) for all samples and, based on the distribution (Additional file 4: Figure S3 ) determined cutoffs for homozygous and heterozygous genotypes (taking into account small deviation from expectation due to sequencing errors). Only samples with at least 1000 reads were considered for further analysis.

For the S200X6.1 locus [33] , the most abundant sequences for each sample were compared against DNA sequences on the NCBI nr database as described above. For samples with ≥ 1000 reads, the species level taxonomy was retrieved and the Anopheles species as well as the total read count were summarized. Sequences generated from the primers targeting the COX1 gene (Culic-Cox1) were used to identify An. nili species in samples as this species is not successfully amplified with the S200X6.1 primers [33, 42] .

DNA was extracted from 665 individual Anopheles mosquitoes collected in Guinea (n = 584) and Mali (n = 81). To characterize the microbiota of each mosquito, we amplified and sequenced the V2 variable region of the bacterial 16S ribosomal RNA genes (see Methods for details). We obtained a total of 8,467,703 sequences derived from 760 samples (665 mosquitoes and 95 extraction controls). On average, each mosquito sample yielded 11,730 sequences (minimum = 259, maximum = 29,908) compared to 5984 sequences on average per extraction control (minimum = 31, maximum = 15,946). We assigned these sequences to 21,527 ASVs (analogue of operational taxonomic units [36] ), representing 37 phyla, including Proteobacteria (6692 ASVs accounting for 64% of all reads), Firmicutes (26%), Actinobacteria (6%), and Bacteroidetes (2%) ( Fig. 1 ; Additional file 5: Figure S4 ).

To investigate differences in bacterial composition among mosquitoes, we calculated β-diversity estimates using weighted UniFrac and Bray-Curtis dissimilarity matrices. PCoA performed using Bray-Curtis dissimilarity or weighted UniFrac distances showed that the Fig. 1 Average relative abundance of bacterial phyla from each mosquito collection site in Guinea and Mali. Bacterial species from the Proteobacteria phylum are the most abundant, followed by Firmicutes and Actinobacteria. Less than 1% abund. represents the aggregate of all phyla that make up < 1% of all bacteria microbial composition separates mosquitoes into distinct clusters ( Fig. 2a ; Additional file 6: Figure S5a ). These clusters appeared to group mosquitoes collected in the same sites (Fig. 2) , and this observation held true when we restricted our analyses to mosquitoes only collected from sites in Guinea ( Fig. 2b ; Additional file 6: Figure S5b ). The details of all the ASVs identified and their taxonomy are provided in Additional file 7: Table S2 .

We simultaneously genotyped the same mosquitoes at loci informative of their species and insecticide resistance status by high-throughput sequencing (see Methods).

Of 665 mosquitoes, 551 (82.9%) were successfully genotyped for the S200X6.1 [33] and cox1 [42] loci. We primarily used the S200X6.1 locus to identify the species of each mosquito as this locus: (i) was more robustly amplified and sequenced than the cox1 locus (with an average read count of 2917 and 1181 per mosquito, respectively) and (ii) provided clearer taxonomic resolution (with, for example, 233-234 [mean 233.67] nucleotides differentiating the sequences from Anopheles gambiae (s.s.) from those of Anopheles coluzzii, compared to 0-4 [mean 1.71] nucleotide differences using the cox1 locus) (Additional file 8: Table S3 ). However, the S200X6.1 locus systematically failed to yield sequences for some mosquitoes that were identified at Anopheles nili using the cox1 sequences. Overall, we identified that the mosquitoes belonged to five Anopheles species: 404 mosquitoes (74.5%) were identified as An. gambiae (s.s.), while the remaining mosquitoes consisted of An. coluzzii (61 mosquitoes, 11.3%), An. melas (57 mosquitoes, 10.5%), An. arabiensis (8 mosquitoes, 1.5%), and An. nili (7 mosquitoes, 1.3%) (Fig. 3) . We also identified five mosquitoes that were heterozygous for the S200X6.1 locus and likely represented F1 hybrids of An. gambiae (s.s.) and An. coluzzii species. The species distribution varied extensively between locations, with An. gambiae (s.s.) accounting > 90.00% of all mosquitoes collected in five of six locations in Guinea, while An. melas was the most abundant species (79.2%, 57/72) in Boffa, a coastal region in (Fig. 3) .

Pyrethroid resistance is often due to a point mutation in the voltage gated sodium channel gene, described as knockdown resistance (kdr) [32] . 550 (82.7%) of the mosquitoes were successfully genotyped at this locus (with an average coverage of 2,436 reads per mosquito). In Guinea, with the exception of mosquitoes collected in Boffa, most mosquitoes (>92.6%) were homozygous for the kdr-w (L1014F) alleles that is associated with resistance to pyrethroids [32] (Fig. 4) . In Boffa, where An. melas is the predominant species, most mosquitoes were homozygous for the wild-type allele (L1014L). In Mali, the distribution was more heterogeneous, with roughly equal proportions of mosquitoes homozygous for the wild-type, resistant allele or heterozygous. Across mosquitoes, the genotype at the kdr-w locus correlated almost perfectly with the mosquito species, with An. gambiae carrying primarily L1014F alleles while An. arabiensis and An. melas were essentially wild-type. Only An. coluzzii showed high proportion of both alleles (Additional file 9: Figure S6 ). The details of all the genotypes and sequences amplified from each mosquito are provided in Additional file 10: Table S4 .

To characterize the composition of the last blood meal of each of these mosquitoes, we used the same DNA extract to amplify and sequence a short fragment of the mammalian mitochondrial 16S rRNA gene. Of these, 133 mosquitoes yielded > 1000 reads and were considered blood fed in later analyses. A total of 126 mosquitoes carried human DNA, 14 carried cow DNA, and two carried sheep DNA (Fig. 5 ; Additional file 11: Table S5 ). Nine mosquitoes fed on more than one mammalian host species (Fig. 5) . The blood-meal composition differed between sites with, for example, 12 mosquitoes (20.1%) from Kankan that fed, at least partially, on cow blood while mosquitoes from all other sites in both Mali and Guinea fed almost exclusively on human blood.

A sequencing-based assay recently developed in our laboratory [34] was then used to determine whether each mosquito carried a eukaryotic parasite and/or arbovirus. After PCR amplification and subsequent sequencing of DNA extracted from individual mosquitoes, we identified DNA sequences from eukaryotes and viruses from 127 mosquitoes, with an average of 1876 reads supporting each identification in each mosquito. Of the 95 extraction controls (water samples that were processed simultaneously), nine were also positive for one or more parasites, but with an average of 54 reads per parasite (Additional file 12: Table S6 ). The low read counts in those water controls could be explained by low-level cross-contamination during extraction or by sequence mis-assignment due to sequencing errors in the Illumina index sequences [34] .

Eight mosquitoes (1.2%) from three sites, six of which were identified as An. gambiae, yielded DNA sequences identical to Plasmodium falciparum, the primary cause of human malaria in Africa. We detected DNA belonging to Theileria species in a relatively high number of mosquitoes (27/665, 4.1%) . Theileria species are protozoan parasites that can be infective to domestic (i.e., cattle) and wild (buffalo) animals and are transmitted by ticks [43] . Interestingly, from the 14 samples for which Bos indicus (cow) was identified with > 500 reads, Theileria species were detected in 12 (86%), suggesting the tick-transmitted parasite may have been ingested by these mosquitoes during their last blood meal. Seven mosquitoes yielded a DNA sequence identical to Loa loa and to several other filarial worms (and were further characterized as deriving from Mansonella perstans after sequencing of longer amplification products; M. Cannon, personal communication). Several DNA sequences were closely related to known parasites of mosquitoes themselves, belonging to Table 2 ). The details of all parasite sequences amplified and their taxonomic information for each mosquito are provided in Additional file 12: Table S6 .

The characterization of the mosquito species, insecticide resistance associated genotype, blood-meal status, and infection from the same mosquitoes that have been examined for their microbial diversity enables a rigorous assessment of the relative contribution of these factors to the microbial composition. Note that to avoid possible biases introduced by sample storage or DNA extraction, we restricted this analysis to mosquitoes collected in Guinea that were all processed identically and simultaneously. We implemented an analysis of variance [48] to simultaneously test the contribution of each factor, while accounting for the others (multivariate analysis). The geographical location of the samples explained most of the variation in microbial composition (R 2 = 0.200, ρ = 0.001), whereas the mosquito species (R 2 = 0.004, ρ = 0.208), kdr-w genotype (R 2 = 0.006 ρ = 0.438), and feeding status (R 2 = 0.003, ρ = 0.173) were not statistically associated with the bacterial composition. The mosquito infection status was significantly associated with the microbial composition but had a very marginal effect (R 2 = 0.015, ρ = 0.001) ( Table 3 ). To further investigate the relative roles of geography and species that are partially confounded in this dataset, we examined PCoAs of the bacteria composition, restricting the analysis to (i) all An. gambiae collected in seven sites (Additional file 13: Figure S7a ) and, separately, (ii) mosquitoes from all Anopheles species collected in Boffa (Additional file 13: Figure S7b ). Together, these analyses validated the results of the statistical testing and confirmed that geographical location of the mosquitoes had a much greater influence on the bacterial composition than did the species of the mosquitoes.

The importance of the mosquito microbiota on vector biology and pathogen transmission has been well recognized, with several studies demonstrating the role of endogenous bacteria on the vector's development [19, 20] , immunity [22, 23] , and competency [49] . However, few studies have examined the factors that shape the bacterial composition of mosquitoes, and most of those used laboratory-reared mosquitoes to assess influences on bacterial communities [20, 30, 50] . This latter limitation could be especially problematic since bacterial diversity of wild-caught Anopheles mosquitoes has been shown to be greater than that of mosquitoes reared in the laboratory [51] . In addition, these studies typically focus on Table 2 Eukaryotic parasites and viruses identified from screening mosquito samples Table 3 Relative contribution of mosquito factors to the microbial composition testing the influence of a single factor on the bacterial composition without accounting for confounding factors.

In the present study, we examined the microbial diversity of 665 individual wild-caught Anopheles mosquitoes collected in six sites in Guinea and one site in Mali. Consistent with previous studies, the bacterial composition of mosquitoes was dominated by bacteria from the phyla Proteobacteria, Firmicutes, and Actinobacteria [26] [27] [28] .

Pseudomonas, a bacterial genus commonly found in mosquito larvae and larval habitats [27] , was one of the most abundant genera detected across all samples (Additional file 5: Figure S4 ; Additional file 7: Table S2 ). Despite these overall similarities, we observed significant differences in bacterial composition among mosquitoes and examined the contribution of various factors to this diversity. For each mosquito, we characterized their species, kdr-w genotype, blood-meal status, and infection with various eukaryotic parasites and viruses. We then simultaneously estimated the relative contribution of each of those endogenous and exogenous factors on the microbial composition of the mosquitoes. In this analysis, the mosquito collection site accounted for approximately 20% of the variation in bacterial composition, whereas the other factors made a marginal or non-significant contribution ( Table 3) .

Our findings are consistent with previous studies that have shown that collection site is a major contributor to the microbial diversity of field-caught Anopheles mosquitoes [26, [52] [53] [54] . For example, Muturi et al. found that sampling site has a strong effect on microbial composition and diversity, even in their examination of nine different mosquito species [27] . It should be noted, however, that "collection site" in our study represents a summary of many parameters. In particular, the mosquitoes were collected from very different ecoregions (grasslands and canopy forests, mountainous forests, or savanna, Table 1 ); as such, the influence of "collection site" on the microbial composition shown in our study could reflect the effect of differences in larval habitats, flora the mosquitoes rely on for nectar feeding, and/or local population differences. In addition, it is worth noting that the bacterial composition of the adult mosquito has been shown to vary depending on the larval breeding sites and the bacterial composition of these aquatic habitats [55] . Also, sugar source appears to have a pronounced influence on the vectoral capacity of An. sergentii mosquitoes [56] and has been shown to impact the microbial composition of laboratory-reared adult mosquitoes [57] . Future studies using a denser, more local sampling of wildcaught mosquitoes will be required to better understand the individual contribution of these local parameters.

On the other hand, our analyses provide new insights on the role of factors other than collection site on the microbial composition of mosquitoes. We did not observe any significant contribution of feeding status on microbial variation of the wild mosquitoes. This finding is in contrast to the results reported in a previous study which described that the bacterial diversity of Aedes aegypti mosquitoes fed on human, chicken, or rabbit blood was significantly lower than that of newly emerged unfed mosquitoes [30] . This discrepancy could reflect differences between mosquito species/genera or, more likely, differences between wild-caught mosquitoes (that might have had prior blood meals) and laboratoryreared mosquitoes with less variable microbial composition. Similarly, our study did not reveal any significant contribution of genetic factors (i.e., mosquito species, kdr-w genotype) on mosquito microbial variation. These observations contrast with those made in a previous study that described distinct bacterial compositions in two species of Culex mosquitoes collected from the same site and with an identical larval aquatic environment [28] ; it is possible that the different results are due to differences among Culex species in their larval feeding habits [58] . In a previous study, the L1014F kdr allele frequency was reported as being high or near fixation in the Kankan and Kissidougou sites of Guinea and low in the Boffa site [59] , which is consistent with our findings. In theory, insecticide-resistant mosquitoes could display a different microbial composition since this resistance may allow them to survive insecticide exposure that could impact the bacterial populations. In our study, we saw no evidence of the L1014F resistant allele influencing adult mosquito microbial composition, although the lack of information on whether these mosquitoes might have been exposed to insecticides limits the conclusions that can be drawn from this observation. However, given that insecticide resistance alleles in genotypes and mosquito species only represent a small fraction of the genetic factors that could impact the mosquito microbiota and in light of our observation that the collection site is strongly associated with the bacterial composition, it would be interesting to further investigate whether genetic diversity is associated with the microbiota of mosquitoes [60] .

Interestingly, we observed a marginal but statistically significant association between infection status (infected n = 127 vs non-infected n = 513) and the mosquito microbial composition. Modification of insect gut microbiota by parasitic [61] or viral [62] infections has been demonstrated in a few studies. Pathogenic or nonpathogenic (e.g., insect-specific viruses) species could be involved in crosstalk with insect metabolism pathways or immune system to influence the microbiota [63] . It should be noted that, due to the low infection rate with parasitic and viral species we found in the mosquitoes (< 5.0%), we assessed the influence of infection on the microbiota using an aggregate of all the parasite and viruses we detected (as opposed to individual parasite and virus species), and it is possible that the effect of one organism on the microbiota might be diluted down and undetected once analyzed together with other parasites and viruses that have no influence on the bacterial communities. For example, Theileria spp. are transmitted by ticks and unlikely to be viable in mosquitoes and, therefore, they probably contribute little or not at all in terms of influencing infection on the mosquito microbiota. Future studies assessing the direct influence of some of the parasites found in abundance in this study (e.g., Parathelohania sp., Microsporidium sp.,) and the recently discovered virus (Anopheles flavivirus) could further elucidate the tripartite relationship between the mosquito, microbiota, and mosquito-infecting agents. One important caveat of our study is that we screened for RNA virus sequences from DNA extracts and that while this approach successfully detected multiple Flaviviruses, the extraction was not optimized for RNA molecules and many sequences might have been lost, leading to an underestimation of the number of viruses.

Finally, the approach described in this study is easily adaptable to other disease vectors (e.g., ticks and sand flies) or insects important in agriculture (such as bees) and easily customizable to examine specific factors of interest by simply adding or replacing PCR primers.

In summary, we provide a comprehensive assessment of the microbial composition and diversity of 665 wild mosquitoes and a simultaneous examination of the relative contribution of five different mosquito-related factors on microbial variation. This approach enables rigorous estimation of the importance of these factors to shaping the bacterial composition, while correcting for their often confounding effect. Our results highlight the prominent role of the mosquito collection site and, to a lesser extent, parasitic and viral infection, on shaping the bacterial composition of wild-caught mosquitoes. These findings provide a solid foundation to implement further investigations and examine the specific components of the environment (e.g., bacterial communities of the larval habitats, source of nectar, genetic diversity) shaping the microbial composition of wild mosquitoes and the mechanisms mediating these effects.

World Health Organization. Mosquito-borne diseases 2020. Geneva: World Health Organization

Mosquito control then, now, and in the future

Biological control of mosquito-borne diseases: the potential of Wolbachia-based interventions in an IVM framework

Insecticide resistance in malaria vectors: an update at a global scale

Global trends in the use of insecticides to control vector-borne diseases

Use of epidemiological and entomological tools in the control and elimination of malaria in Ethiopia

Evidence of insecticide resistance selection in wild Anopheles coluzzii mosquitoes due to agricultural pesticide use

Impact of environment on mosquito response to pyrethroid insecticides: facts, evidences and prospects

Changes in the feeding behaviour of a malaria vector, Anopheles farauti Lav., following use of DDT as a residual spray in houses in the British Solomon Islands Protectorate

Changes in Anopheles funestus biting behavior following universal coverage of long-lasting insecticidal nets in Benin

Increasing outdoor host-seeking in Anopheles gambiae over 6 years of vector control on Bioko Island

Increased proportions of outdoor feeding among residual malaria vector populations following increased use of insecticide-treated nets in rural Tanzania

Permethrin-impregnated bednet effects on resting and feeding behaviour of lymphatic filariasis vector mosquitoes in Kenya

Increasing the potential for malaria elimination by targeting zoophilic vectors

Bacillus thuringiensis: a story of a successful bioinsecticide

Immune activation by lifeshortening Wolbachia and reduced filarial competence in mosquitoes

Stable introduction of a life-shortening Wolbachia infection into the mosquito Aedes aegypti

Mosquito/ microbiota interactions: from complex relationships to biotechnological perspectives

Mosquitoes host communities of bacteria that are essential for development but vary greatly between local habitats

Mosquitoes rely on their gut microbiota for development

Delayed larval development in Anopheles mosquitoes deprived of Asaia bacterial symbionts

Implication of the mosquito midgut microbiota in the defense against malaria parasites

The Anopheles innate immune system in the defense against malaria infection

Chromobacterium Csp_P reduces malaria and dengue infection in vector mosquitoes and has entomopathogenic and in vitro anti-pathogen activities

Natural microbe-mediated refractoriness to Plasmodium infection in Anopheles gambiae

Bacterial diversity of field-caught mosquitoes, Aedes albopictus and Aedes aegypti, from different geographic regions of Madagascar

Mosquito microbiota cluster by host sampling location

Culex pipiens and Culex restuans mosquitoes harbor distinct microbiota dominated by few bacterial taxa

Whole metagenome sequencing reveals links between mosquito microbiota and insecticide resistance in malaria vectors

Host blood-meal source has a strong impact on gut microbiota of Aedes aegypti

Factors shaping the gut bacterial community assembly in two main Colombian malaria vectors. Microbiome

Detection of the East and West African kdr mutation in Anopheles gambiae and Anopheles arabiensis from Uganda using a new assay based on FRET/Melt curve analysis

Insertion polymorphisms of SINE200 retrotransposons within speciation islands of Anopheles gambiae molecular forms

A high-throughput sequencing assay to comprehensively detect and characterize unicellular eukaryotes and helminths from biological and environmental samples. Microbiome

Ultrahigh-throughput multiplexing and sequencing of > 500-basepair amplicon regions on the Illumina HiSeq 2500 platform. mSystems

DADA2: high-resolution sample inference from Illumina amplicon data

Silva taxonomic training data formatted for DADA2. Silva version 132 ed: Zenodo

phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data

FLASH: fast length adjustment of short reads to improve genome assemblies

In silico assessment of primers for eDNA studies using PrimerTree and application to characterize the biodiversity surrounding the Cuyahoga River

Unbiased characterization of Anopheles mosquito blood meals by targeted high-throughput sequencing

Direct PCR of indigenous and invasive mosquito species: a time-and cost-effective technique of mosquito barcoding

Theileria: intracellular protozoan parasites of wild and domestic ruminants transmitted by ixodid ticks

Insect pathology

Morphological characterization of Setaria worms collected from cattle

Functional analysis of a novel parasitic nematode-specific protein of Setaria digitata larvae in Culex quinquefasciatus by siRNA mediated RNA interference

West African Anopheles gambiae mosquitoes harbor a taxonomically diverse virome including new insect-specific flaviviruses, mononegaviruses, and totiviruses

Community Ecology Package. R package version 2.5-6 ed2019

Mosquito microbiome dynamics, a background for prevalence and seasonality of West Nile Virus

Diverse laboratory colonies of Aedes aegypti harbor the same adult midgut bacterial microbiome

Bacterial diversity analysis of larvae and adult midgut microflora using culture-dependent and culture-independent methods in lab-reared and field-collected Anopheles stephensi-an Asian malarial vector

Diversity of the bacterial microbiota of Anopheles mosquitoes from Binh Phuoc Province

Seasonality and locality affect the diversity of Anopheles gambiae and Anopheles coluzzii midgut microbiota from Ghana

Bacterial associations reveal spatial population dynamics in Anopheles gambiae mosquitoes

Midgut microbiota of the malaria mosquito vector Anopheles gambiae and interactions with Plasmodium falciparum infection

Natural plant sugar sources of Anopheles mosquitoes strongly impact malaria transmission potential

The environment and species affect gut bacteria composition in laboratory co-cultured Anopheles gambiae and Aedes albopictus mosquitoes

Larval behavior of four Culex (Diptera: Culicidae) associated with treatment wetlands in the southwestern United States

Species identification and resistance status of Anopheles gambiae s.l. (Diptera: Culicidae) mosquitoes in Guinea

The mosquito microbiota influences vector competence for human pathogens

Parasites modulate the gut-microbiome in insects: a proofof-concept study

Increase in gut microbiota after immune suppression in Baculovirus-infected larvae

The mosquito holobiont: fresh insight into mosquito-microbiota interactions. Microbiome

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Supplementary information accompanies this paper at https ://doi. org/10.1186/s1307 1-020-04491 -7.Additional file 1: Figure S1 : Geographical locations (green pins) of mosquito collection sites in Guinea and Mali. Map image was prepared using the online ArcGIS ® software by ESRI.Additional file 2: Table S1 : Summary of all primers used in the study. Table shows , for each primer pair, the loci targeted, the base pair length of amplicon, and the forward and reverse primer sequences.Additional file 3: Figure S2 : Distribution of reads counts for the mammalia_16S primer across mosquito samples and negative control samples.Additional file 4: Figure S3 : Distribution of the Seq2/(Seq1+Seq2) ratio across samples with one or more reads for the KDR primer used to determine genotype for the kdr locus. Samples with a ratio < 0. 15 Table S2 : ASV taxonomy. Table shows, for each ASV, the sequence identified, the taxonomic information (Kingdom to Species), and abundance values. NA represents when an ASV is unknown at that taxonomic rank Additional file 8: Table S3 : Pairwise comparison of cox1 and S200X.6 primers for resolution of mosquito species. Table shows, for each pairwisecomparison between two Anopheles species, the mean nucleotide difference and range for cox1 and S200X.6 loci Additional file 9: Figure S6 : Distribution of L1014F mutation (kdr-w) in mosquitoes grouped by Anopheles species. Numbers above each bar represent the total number of mosquitoes that kdr-w genotype is identified, per site. RR Homozygous resistant, SS homozygous sensitive, R/S heterozygous.Additional file 10: Table S4 : Summary of species and kdr-w determination. Table shows, for each sample, the mosquito collection site, the number of reads belonging to kdr alleles and Anopheles species and their sequences. kdr_w Knockdown resistance west (mutant), WT wildtype, H2O water controls. Seq sequence Additional file 11: Table S5 : Summary of host blood-meal composition. Table also gives, for each sequence, the frequency per sample (count), percentage match to NCBI database (% Identity), length of sequence (in bp), and length of the match (in bp) to NCBI sequence