Mosquito vector‐associated microbiota: Metabarcoding bacteria and eukaryotic symbionts across habitat types in Thailand endemic for dengue and other arthropod‐borne diseases

Abstract Vector‐borne diseases are a major health burden, yet factors affecting their spread are only partially understood. For example, microbial symbionts can impact mosquito reproduction, survival, and vectorial capacity, and hence affect disease transmission. Nonetheless, current knowledge of mosquito‐associated microbial communities is limited. To characterize the bacterial and eukaryotic microbial communities of multiple vector species collected from different habitat types in disease endemic areas, we employed next‐generation 454 pyrosequencing of 16S and 18S rRNA amplicon libraries, also known as metabarcoding. We investigated pooled whole adult mosquitoes of three medically important vectors, Aedes aegypti, Ae. albopictus, and Culex quinquefasciatus, collected from different habitats across central Thailand where we previously characterized mosquito diversity. Our results indicate that diversity within the mosquito microbiota is low, with the majority of microbes assigned to one or a few taxa. Two of the most common eukaryotic and bacterial genera recovered (Ascogregarina and Wolbachia, respectively) are known mosquito endosymbionts with potentially parasitic and long evolutionary relationships with their hosts. Patterns of microbial composition and diversity appeared to differ by both vector species and habitat for a given species, although high variability between samples suggests a strong stochastic element to microbiota assembly. In general, our findings suggest that multiple factors, such as habitat condition and mosquito species identity, may influence overall microbial community composition, and thus provide a basis for further investigations into the interactions between vectors, their microbial communities, and human‐impacted landscapes that may ultimately affect vector‐borne disease risk.

Unfortunately, while studies on the application of microbes such as Wolbachia to control human disease have advanced, we still lack basic insight into the natural microbial communities associated with vectors, from viruses to bacteria to single-celled eukaryotes. Given that infectious vector-borne diseases are increasing in humans (Jones et al., 2008), our goal is to characterize the composition and structure of microbiota in mosquito vectors across landscapes altered by human activities, where differences may influence the biology of vectors, their vectorial capacity, and ultimately the outcome of microbe-mediated disease control (Hughes et al., 2014).
Living within and upon mosquitoes are numerous, and diverse microscopic life-forms, including viruses, bacteria, fungi, protozoa and nematodes, collectively referred to as their microbiota. Recent improvements in the exploration of microscopic organisms through the use of metagenomic PCR (Hugenholtz, Goebel, & Pace, 1998;Liles, Manske, Bintrim, Handelsman, & Goodman, 2003) and nextgeneration sequencing (Logares et al., 2012) have allowed more accurate and sensitive documentation of these organisms. Several studies have characterized and compared mosquito bacterial microbiota in natural habitats (Akorli et al., 2016;Buck et al., 2016;Gimonneau et al., 2014;Minard et al., 2012Minard et al., , 2014Minard et al., , 2015Tchioffo et al., 2015;Yadav et al., 2015;Zouache et al., 2010) but few have characterized eukaryotes and in such cases have targeted fungi (Muturi, Bara, Rooney, & Hansen, 2016;Muturi, Ramirez, Rooney, & Dunlap, 2016;Steyn, Roets, & Botha, 2016). We document mosquito microbiota composition and diversity for vector species Aedes aegypti, Aedes albopictus, and Culex quinquefasciatus across ecologically distinct, yet geographically proximate, habitats in Thailand where many mosquitoborne diseases circulate. Significant diseases in the area vectored by Aedes include dengue (Limkittikul, Brett, & L'Azou, 2014;Nisalak et al., 2016), chikungunya (Chusri et al., 2014;Wanlapakorn et al., 2014), and Zika (Musso & Gubler, 2016). In addition, Cx. quinquefasciatus, amongst other Culex spp., locally transmits Japanese encephalitis, an important cause of encephalitis in Thailand (Olsen et al., 2015). In characterizing mosquito species and their diversity across habitats, we have observed that mosquito overall diversity decreased with increasing relative abundance of several medically important invasive mosquito vector species and degree of human habitat modification (Thongsripong et al., 2013). In this study, we describe the microbial communities of these disease vectors in the subset of habitats where they were found in our previous study, that is, in rural, suburban, and urban habitats. We employ 16S and 18S rRNA amplicon sequencing, a technique referred to as metabarcoding (Taberlet, Coissac, Hajibabaei, & Rieseberg, 2012), to simultaneously characterize the bacterial and eukaryotic microbiota of vector species-specific pools of three to 25 individual mosquitoes. This study is an important first step in describing the microbial communities associated with different vector species in different habitats. Our results add to the current knowledge of vector-associated microbiota, especially the eukaryotes, to help reveal microbial species with the potential to control mosquito-borne diseases, and to emphasize the need to understand how vector symbiotic microbial communities vary across habitats and may ultimately affect vector-borne disease transmission.

| Mosquito collection and identification
Adult mosquito collection and habitat characterization was described in a previous study (Thongsripong et al., 2013). In short, we collected adult mosquitoes using a combination of trap types (BG sentinel, Mosquito Magnet, CDC light trap, and CDC backpack aspirator) during the rainy season of 2008 from different habitat types ranging from forest to urban in Nakhon Nayok, a small province situated in microbial community composition, and thus provide a basis for further investigations into the interactions between vectors, their microbial communities, and humanimpacted landscapes that may ultimately affect vector-borne disease risk.

K E Y W O R D S
biodiversity, metabarcoding, microbiota, mosquito, vector, vector-borne diseases central Thailand approximately 60 km northeast of Bangkok. Over 83,000 adults mosquitoes were collected as described in the previous study (Thongsripong et al., 2013). Trap collections from each site were transported on dry ice to a local field station. Upon arrival, mosquitoes were sorted immediately and then stored at −20°C for less than one month before they were transported on dry ice to Mahidol University in Bangkok. The female mosquitoes were identified to species using available morphological keys (Rattanarithikul, Harbach, Harrison, Panthusiri, & Coleman, 2007;Rattanarithikul, Harrison, Panthusiri, Peyton, & Coleman, 2006;Rattanarithikul et al., 2010;Rattanarithikul, Harbach, et al., 2005;Rattanarithikul, Harrison, Harbach, et al., 2006;Thongsripong et al., 2013). Adult females used in this study were not visibly blood-engorged. We focused on adult females because they can transmit pathogens to vertebrates through blood-feeding, and therefore, species identification keys have been better developed for females. For each taxon identified, three to twenty specimens were vouchered as references and housed at the Center of Excellence for Vectors and Vector-Borne Diseases, Faculty of Science, Mahidol University at Salaya, Nakhon Pathom, Thailand.
Female mosquitoes of the same species were grouped according to their site and date of collection. Of the over 100 mosquito taxa collected, whose mosquito community composition varied significantly across habitat (Thongsripong et al., 2013), we subsampled ten mosquito pools of the species Ae. aegypti, Ae. albopictus, and Cx. quinquefasciatus, restricted to rural, suburban, or urban habitats, in pool sizes of three to 25 depending on their abundance (Table 1), to characterize their microbiota (Table 2).

| DNA extraction from mosquito pools
Mosquitoes were homogenized in 250 μl of 1× phosphate-buffered saline using Tissue Lyser II (QIAGEN, CA, USA); 1-2 stainless steel beads (5 mm diameter) were added to each sample and trituration performed at 30 Hz for 3 min with one minute on ice halfway through to cool samples. The homogenate was divided, with half (125 μl) frozen at −80°C for future studies, and the other half mixed with 875 μl of TRIzol LS (Invitrogen, CA, USA). RNA was extracted from samples as per manufacturer's instructions. Following RNA extraction and phase separation, total DNA was isolated from the interphase/organic phase using back extraction buffer (BEB; 4 M guanidine thiocyanate, 50 mM sodium citrate, and 1 M Tris free base). In short, 500 μl of BEB was added to the interphase/organic phase, and the solution was vigorously mixed and then incubated at room temperature for 10 min. The solution was then subjected to centrifugation at 12,000 g and 4°C for 10 min, and the upper aqueous phase was collected. Next, we added 400 μl of isopropanol and incubated the samples at room temperature for 5 min. The DNA pellet was precipitated from the solution by centrifugation at 12,000 g for 10 min at 4°C. Each DNA pellet was washed twice with 75% ethanol, air-dried, and dissolved in 50-100 μl of 8 mM NaOH.
An appropriate amount of 0.1 M HEPES was added to the sample to obtain a neutral solution, and the DNA was stored at −30°C.

| Metabarcoding: 16S and 18S rRNA amplicon library preparation and pyrosequencing
To characterize the noneukaryotic and eukaryotic microbiota of these mosquitoes, we amplified and sequenced 16S and 18S ribosomal RNA (rRNA) genes, respectively, using titanium chemistry on the 454 Genome Sequencer FLX instrument (Roche, Switzerland). The variable region V3 of the 16S rRNA gene was amplified with universal primers U341F (5′ CCT ACG GGR SGC AGC AG 3′) and U800R (5′ CCR GGG TAT CTA ATC C 3′) designed to bind all known bacteria and archaea (Wang & Qian, 2009), and a region of the 18S rRNA gene was amplified with eukaryotic microbial universal primers Euk82F (5′ GAA ACT GCG AAT GGC TC 3′; López-García, Philippe, Gail, & Moreira, 2003) and Euk516R (5′ ACC AGA CTT GCC CTC C 3′; Diez, Pedros-Alio, Marsh, & Massana, 2001). Additionally, the forward primers contain primer A sequence (5′ CGT ATC GCC TCC CTC GCG CCA TCA G 3′) followed by a ten-base molecular identifier (MID) sequence tag (Table   S1), while the reverse primer contained primer B sequence (5′ CTA TGC GCC TTG CCA GCC CGC TCA G 3′) followed by a ten-base MID tag ( Amplicons were purified using PCR purification kits following the manufacturer's protocol (NucleoSpin ExtractII; Macherey-Nagel, Germany). We measured the purified DNA concentration and distribution of amplicon length using the Agilent Bioanalyzer 2100 instrument and Agilent DNA 1000 reagent and kit (Agilent Technologies, CA, USA). Amplicons were combined based on DNA concentration in equimolar ratios. Pyrosequencing was carried out on one-fourth of a PicoTiterPlate (PTP) using primer B (producing reverse reads; R2).

| Sequence quality control and generation of operational taxonomic units
To characterize the number of distinct microbial taxa present in each sample, that is, species richness, we grouped sequence reads into the same operational taxonomic unit (OTU) if they were equal to or less than 3% divergent, a standard threshold which buffers against overestimation of microbial diversity (Kunin, Engelbrektson, Ochman, & Hugenholtz, 2010). Despite a relatively low error rate, the large number of reads resulting from next-generation sequencing can include a substantial number of erroneous sequences and thus an overestimation of total OTUs (Gilles et al., 2011). Consequently, it is critical to distinguish true sequence diversity from noise introduced by experimental procedures. We adopted strict filtering and quality control procedures prior to the analysis steps in order to limit noise due to The estimation of tree cover and amount of trash and clutter were averaged between two same observers in all sites.
c The number of humans and cars that travelled passed/into the study sites were counted at around noon during the weekdays for half an hour.
sequencing errors (described below). We utilized a sequence analysis pipeline: QIIME, http://qiime.org (Caporaso, Kuczynski, Stombaugh, & Bittinger, 2011) for sequence processing and OTU generation. The command lines used in the analyses are given in Figure S1. Default settings were used with the exception of minor changes described below.
Analysis began with assignment of reads to samples based on MIDs followed by initial sequence quality check. Suitable quality includes (i) no mismatches to either MID and reverse primer sequences, (ii) meeting the mean quality score for sequences, (iii) sequence length of at least 200 bp, (iv) no ambiguous called bases, and (v) no homopolymers over 9 bp. We did not discard entire sequence reads containing low-quality score sections, which would have resulted in too few remaining sequences; instead, those sequences were truncated to exclude the poor-quality section. After initial quality check, we then employed QIIME's algorithms that approximate PyroNoise (Quince et al., 2009;Reeder & Knight, 2010) to cluster the flowgram (analogous to the trace data in Sanger sequencing) and denoise the data. Next, we removed forward primer sequences from the end of reads, allowing at most one mismatch to the primer sequence, again to minimize the number of discards, where reads with more than one mismatch to the forward primer sequences were removed. Base calling also becomes less accurate toward the end of the reads (Ledergerber & Dessimoz, 2011). After the quality filtering processes, sequences were aligned (Caporaso et al., 2010) and grouped into OTUs based on ≤3% divergence. Specifically, within this step, we employed "pick_otus.py" command in QIIME to cluster reads de novo and maintained the default options, including the use of UCLUST as the OTU picking method (Edgar, 2010). Next, the OTUs identified as chimeric by the Chimera Slayer algorithm (Haas et al., 2011) were removed. There was no chloroplast or mitochondrial sequences present after this step.

| Alpha-diversity analysis
As diversity indices are sensitive to sample size, in this case number of reads analyzed (Magurran, 1988), we used two approaches. First, rarefaction analysis was used to project the rate of accumulation of OTUs with the addition of more sequence reads, resulting in rarefaction curves. Second, equal numbers of reads were randomly subsampled from each original set of reads prior to OTU estimation. Samples with fewer than 25% of the average number of reads were excluded, resulting in approximately 7,100 and 11,000 reads for analysis of 16S and 18S rRNA, respectively. This set of OTUs was used to calculate the species diversity indices Shannon entropy (Shannon, 1949) and diversity estimator Chao1 (Chao, 1984) for across-sample comparisons.

| Determination of mosquito microbiota composition
To determine the taxonomic identities of OTUs, we used QIIME's command line "assign_taxonomy.py" and maintained the default options (Wang, Garrity, Tiedje, & Cole, 2007). We used the Greengenes database for 16S rRNA sequence analysis (DeSantis et al., 2006) and SILVA taxonomy sequence database for 18S rRNA analysis (Yilmaz et al., 2013). Note that although 16S primers were designed to amplify both bacteria and archaea, no archaea were recovered. Archaea have no known associations with mosquitoes (Clemens, 2012), associating primarily with wood and detritus feeders such as termites and beetles (Engel & Moran, 2013). Results based on 16S will herein address bacterial diversity.

| Sequence processing
A total of 273,590 DNA sequence reads of the 16S and 18S rRNA genes were generated from the 10 mosquito-pool samples, varying in number by host species and habitat type (

| Richness of microbiota associated with mosquito vector species
The total number of bacterial OTUs found in all samples was 90, and the median estimates of bacterial OTUs per sample were 18. For eukaryotic microbiota, the total estimated OTUs were 53, and the median estimate of OTUs was 12.5. Rarefaction curves for 16S rRNA and  Table S2.
F I G U R E 1 16S and 18S rRNA OTU richness of mosquito-associated microbial communities collected from rural habitats (a and b, respectively). Rarefaction curves by vector species in the rural habitat where all species were collected for pool sizes of 25. Gray shading around lines represents the 95% confidence interval of the estimate

| Richness of microbiota associated with mosquito vectors by habitat type
In this study, the bacterial microbiota of Ae. aegypti had a higher richness in rural than in suburban habitats based on rarefaction curves of OTUs from 16S rRNA reads for the same pool size of 25 ( Figure 2a). OTU richness appears to be lowest in the urban habitat according to pool-of-10 estimates, but direct comparison for the same pool size cannot be made as the pool-of-25 urban estimate is uncertain due to a lack of reads ( Table 2). The eukaryotic microbiota of Ae. aegypti also appears to be slightly more OTU-rich for a given pool size (25) in rural habitats than in suburban and urban ( Figure 2b). Across-habitat comparison of OTU richness for Cx. quinquefasciatus samples was not possible due to insufficient sequence recovery, nor for Ae. albopictus samples which were only collected in the rural habitat.
Pool size did affect OTU richness estimates of bacteria by rarefaction analysis, where the larger the pool size, the higher the richness estimated. Rarefaction curves based on 16S rRNA for Ae. albopictus samples all from the rural habitat but differing in pool sizes (3, 9, and 25 mosquitoes) projected higher OTU richness with increasing pool size ( Figure S2a). This trend was not consistently observed for estimates of eukaryote richness ( Figure S2b).

| Diversity of microbiota associated with mosquito vectors
Diversity, a function of not only species richness (number of species) but also species relative abundance, here measured with the diversity indices Shannon entropy (Shannon, 1949), and Chao1 (Chao, 1984), varied greatly by sample (Figure 3). Chao indices of diversity ranged from 11 to 41 with an average of 19 OTUs across samples.
The Shannon indices of diversity ranged from 0.13 to 2.30 with an average of 0.73.
In this study, diversity indices suggest that bacteria increase in diversity from the samples collected in the urban to rural sites for a given vector species, in this case Ae. aegypti (Figure 3a,b). Bacterial diversity associated with Ae. albopictus from rural sites also increased with pool size (Figure 3a,b). Diversity indices based 18S rRNA were inconsistent, except for the Shannon index, where eukaryote microbial diversity was again highest in the rural site for a given vector species, in this case Ae. aegypti, and pool size (25; Figure 3c,d).

| Changes in microbiota community across vector species and habitat type
Comparison of mosquito microbiota or beta-diversity across samples using cluster analysis in this study indicated significant differences to a certain extent by both species and the habitats characterized in this study ( Figure 4). The grouping of samples by host species was also observed somewhat in Nonmetric Multidimensional Scaling (NMDS) plots at least for bacterial community composition ( Figure S3). Bacterial community composition was similar within a given vector species for Ae. aegypti and Ae. albopictus, but not so in the case of Cx. quinquefasciatus, whose microbiota shared similarities with both Aedes species rather than being Culex-specific ( Figure 4a). Interestingly, sequences identified as Wolbachia spp. constituted the majority of Ae. Albopictusand Cx. quinquefasciatus-derived sequences ( Figure 5); when these F I G U R E 2 16S and 18S rRNA OTU richness of Aedes aegypti-associated microbial communities across habitat type. The expected numbers of bacterial (a) and eukaryotic (b) OTUs associated with Ae. aegypti collected from different habitats were plotted as a function of sampling effort (in this case, number of sequences were excluded, there was no consistent clustering by either vector species or habitat type ( Figure S4).
Eukaryotic microbiota clustered by vector genus, with Cx. quinquefasciatus exhibiting a more distinct eukaryotic microbiota than the Aedes spp., regardless of habitat. When the most common sequences, identified as Ascogregarina spp. (Figure 6), were excluded from Ae. aegypti and Ae. albopictus sample analysis, no clustering pattern based on vector genus, species, or habitat type was observed ( Figure S4).

| Microbiota composition of the mosquito vectors
The identities of bacterial and eukaryotic OTUs as determined by our database search algorithms implemented in QIIME are summarized to major taxonomic group in Figures 5 and 6 (see Tables S3 and S4 for (Minard, Mavingui, & Moro, 2013).
Those that applied the culture-independent deep sequencing approach were mainly restricted to gut bacterial communities (Minard et al., 2012;Wang, Gilbreath, Kukutla, Yan, & Xu, 2011) or other isolated adult arthropod parts (Damiani et al., 2010;Sharma et al., 2014), although one study investigated bacteria diversity from whole Ae. albopictus (Minard et al., 2014). Studies using shotgun metagenomic sequencing (which is agnostic to the genetic material sequenced, as opposed to the gene-specific amplicon-based method used here) have identified both bacterial and eukaryotic microbes associated with mosquitoes, however with insufficient depth to perform diversity analyses (Chandler, Liu, & Bennett, 2015). Thus, our study contributes a uniquely comprehensive view of mosquito-associated microbial communities, unrestricted to bacteria, whose dynamics across vectors and landscapes could ultimately inform the distribution of infectious diseases.
Overall richness of bacterial OTUs in our study, estimated as 18 OTUs, was much lower than comparable studies using the same OTU cutoff (3% genetic divergence). For example, a median of 42 bacterial OTUs was estimated from the guts of eight mosquito species (Minard et al., 2012); higher still by at least 200 OTUs was the bacterial richness in the gut of a single adult Anopheles gambiae (Wang et al., 2011).
This enormous variability amongst studies, including ours, likely reflects inherent differences in microbiota diversity, composition, and Likewise, the most common eukaryotic microbes identified in this study have also been previously found associated with mosquitoes.
For example, an OTU most likely classified to Cladosporium sp., which was found in nine of the ten samples, has been found in Culex pipiens from California (Chandler et al., 2015).

Interestingly, our results indicate infection with both a eukaryotic
and a bacterial pathogen in the urban Cx. quinquefasciatus sample.
An OTU most likely classified as Pandora/Furia, which comprises the majority of the eukaryotic reads in this sample, is a known and exclusive insect pathogen (Gryganskyi et al., 2012). Likewise, Novispirillum, which comprises over 96% of the total bacterial reads in this study, has been isolated from unhealthy-looking Cx. pipiens larvae in another study (Tranchida et al., 2012). Notably, experimental inoculation of mosquito larvae with Novispirillum only caused mortality when the cuticle was damaged by co-inoculation with a parasitic nematode or by needle puncture (Tranchida et al., 2012). As entomopathogenic fungi infect the host by penetrating the cuticle (Shah & Pell, 2003), this suggests the interesting scenario that an initial infection by Pandora/Furia overlapping ranges and habitats of these two mosquito species (Coon et al., 2016). Novel Wolbachia-mosquito associations have also been reported in An. gambiae collected in West Africa (Baldini et al., 2014).
Given that Wolbachia can spread to new populations amongst arthropods within and between species through horizontal transmission (Ahmed, Breinholt, & Kawahara, 2016;Boyle, Oneill, Robertson, & Karr, 1993), its resulting distribution could conceivably be quite heterogeneous across different populations within a species. Despite the preliminary nature of our findings, the implications that Wolbachia may naturally occur in some Ae. aegypti populations are worth briefly discussing. Wolbachia can affect vector competence (Jupatanakul, Sim, & Dimopoulos, 2014) and vectorial capacity by reducing mosquito populations (Atyame et al., 2016). In the former case, dengue virus is blocked in experimentally transfected Ae. aegypti with wMel and wMel-Pop strains of Wolbachia (Hoffmann, Ross, & Rai, 2015): These mosquitoes are being released to suppress dengue transmission (Frentiu et al., 2014). In the latter case, experimentally transfected Ae. albopictus have been developed and released into Ae. albopictus populations with naturally occurring Wolbachia for population suppression (Atyame et al., 2016;Mains, Brelsfoard, Rose, & Dobson, 2016) by way of a mechanism known as cytoplasmic incompatibility (CI), where matings between males and females with different strains of Wolbachia experience mortality of the developing embryos (Zabalou et al., 2004).
Thus, naturally occurring Wolbachia incompatible with introduced transfected strains could interfere with virus suppression and could also lead to local Ae. aegypti population suppression. Further study is warranted, not only to confirm results from this study but also, given other records of naturally occurring Wolbachia in Ae. aegypti (Coon et al., 2016) findings, to delineate the specificity, frequency, geographic breadth, and the disease suppression and vector control implications of this potentially novel vector-microbe association We recognize that some limitations exist in our study, which affect the generality of our findings and their interpretations. First, because we wanted to maximize the breadth of locations and species in this study for a broad survey, we did not include replications sufficient to statistically test associations by site and species. Similarly, we employed multiple trap types to attract a range of mosquito species, but these also target mosquitoes at different physiological states (e.g., recently bloodfed, resting, or host-seeking). Although some studies suggest that physiological state may affect the composition and diversity of midgut microbiota (Wang et al., 2011;Yadav et al., 2016), others have demonstrated that Ae. aegypti harbors a stable bacterial community during its adult life (David, dos Santos, Vicente, & Maciel-de-Freitas, 2016). Nevertheless, to minimize the impact of bloodfed status on our microbiota characterization, we excluded visibly bloodfed females from this study, and thus feel that the contribution of trap type to error variation in vector microbiota is minimal. Finally, in this study, we wanted to characterize the breadth of microbes associated with a species, and therefore, pooling of individuals was appropriate. Although pooling resulted in the loss of information at the individual mosquito level, others have observed individuals to be highly variable in field-caught specimens (Osei-Poku et al., 2012).
In characterizing microbiota diversity across vectors, we acknowledge that sample preparatory steps including the use of phenolic extraction and differing PCR protocols may have an impact on our data. A previous study indicated that the PCR conditions that we used could introduce biases in the next-generation sequencing, including a higher risk of errors and chimera production (Ahn, Kim, Song, & Weon, 2012). To limit biases introduced by PCR, we performed a strict sequence quality control. In the future, the use of mock community control, although not common practice at the time of this study, is recommended to quantify sequence contamination (Lusk, 2014;Salter et al., 2014). However, the latter study also show that contamination is more pronounced with samples containing low microbial biomass. As mosquitoes are known to harbor high titers of bacteria and eukaryotic microbes (Minard et al., 2013), we believe that contamination contributed minimally to our data.
In this study, the most important factor affecting mosquito microbiota was vector species, and secondarily habitat. In spite of variation in microbial diversity and composition between samples, we observed  (Minard et al., 2014).
However, we expect that richness indices, rather than diversity indices that measure relative abundance, are more robust to the "saturation" effect of Wolbachia or Ascogregarina read numbers.
Thus, vector microbiota includes common obligate endosymbionts that tend to be vector-specific, as well as rare and/or facultative microbes whose distribution is heterogeneous across vectors but may depend in part on habitat. Indeed, for a given vector species, microbiota varied across different habitats. Mosquitoes are thought to acquire their microbiota from the environment in which the species develops and lives (although some bacteria, such as Wolbachia, are also vertically transmitted). For example, mosquitoes acquire Ascogregarina sp.
from their larval habitat (Lantova & Volf, 2014). Many other bacteria are shared within and between species through the larval aquatic habitat (Coon et al., 2016). As a result, the diversity and composition of a vector's microbiota are determined by, at least partially, the diversity and composition of microbes in their habitat, which can in turn be highly variable (Coon et al., 2016). In our study, composition and diversity in terms of rarefaction asymptote estimates of OTU richness, indices of diversity, and number of identified taxa, particularly for bacteria, showed an increasing trend in diversity from urban to rural sites.
This trend may be due to greater complexity in terms of heterogeneity of rural habitats relative to urban sites, supporting more microbe species, by various potential mechanisms for species distributions within and across landscapes (Tscharntke et al., 2012). In contrast, a study of field-caught mosquitoes from Ghana found that midgut microbiota was more diverse in urban mosquitoes compared to rural ones (Akorli et al., 2016). However, Akorli et al. (2016) sampled sites that were much farther apart and in markedly different ecological zones (forestsavannah transitional zone vs. coastal savannah zone) compared to our study. Thus, differences they observed in mosquito microbiota could be as much a function of ecological zone as urbanization.
The results of this study suggest that in addition to biodiversity loss occurring with land-use change and habitat degradation in terms of mosquito diversity (Thongsripong et al., 2013), similar patterns may exist at the microbial level, in terms of the microbiota of mosquitoes. Loss of microbial diversity has been observed in other habitats modified by humans and in some cases associated with increased disease risk. A shift in soil microbial communities has been observed when forest is transformed to pasture (Nüsslein & Tiedje, 1999); reduced soil microbial diversity has also been reported in arable land (Torsvik, Ovreas, & Thingstad, 2002). Agricultural transformation of land also favors pathogenic Burkholderia pseudomallei in soil communities (Kaestli et al., 2009); this can include higher abundance at sites in agricultural versus unfarmed land (Limmathurotsakul et al., 2010). Notably, ecological differences in the distributions of certain mosquitoes, such as Ae. aegypti, an invasive species which increases relative to Ae. albopictus as habitats become increasingly modified by humans (Thongsripong et al., 2013), will also affect the distribution of specific microbial agents such as Wolbachia, species of which occur natively in Ae. albopictus and are known to interfere with infectious agents such as dengue, chikungunya, and Zika viruses. Thus, shifts in the communities of nonpathogenic vector-associated microbes in response to ecological changes and distribution of vectors could affect the distribution of pathogens across human-modified landscapes.
The relationship between increased threats to human health via infectious microbes and land-use change and/or biodiversity loss has been reported in several systems, linked to urbanization, agricultural intensification, and deforestation (McFarlane, Sleigh, & McMichael, 2013;Olson, Gangnon, Silveira, & Patz, 2010;Vitter et al., 2006;Wilcox & Colwell, 2005). Our results highlight the importance of subsequent investigation to determine how symbiotic microbial diversity of mosquito vectors responds to habitat change and the implications for human health.
Because most female mosquitoes feed on vertebrates as a requisite for egg maturation, we wanted to explore the impact of vertebrate diversity on microbiota diversity and composition. Others have shown that mosquito midgut bacterial communities are dynamic and vary significantly with respect to bloodfed state (Wang et al., 2011).
Host DNA is detectable by sequencing even when engorgement is no longer visible (Davey, Casey, Burgess, & Cable, 2007;Lee et al., 2002), and quite interestingly, we identified 18S rRNA sequences belonging to phylum Chordata, class Mammalia, in a few of our mosquito samples. Although their presence was not frequent enough across pools to conduct further analyses, with the everimproving depth of sequencing provided by technological advances, we envision future studies being able to simultaneously characterize both microbial community composition and vertebrate host utilization from field-collected mosquitoes.

| CONCLUSIONS
This field study characterizes both bacterial and eukaryotic microbiota associated with naturally occurring mosquito vector species using culture-independent methods. The diversity of microbiota associated with mosquitoes collected in our study differed both between vector species and across environments within a given species. Importantly, for Ae. aegypti, microbiota diversity decreased along a gradient of increasing habitat modification from rural to urban, warranting further study to confirm this finding and to investigate the role of microbial symbiont diversity in vector-borne disease transmission. We found the composition of mosquito microbiota to be largely vector-specific and often dominated by a few taxa, such as Ascogregarina and Wolbachia.
Our multidomain metabarcoding approach also revealed potential bacterial and microbial eukaryote interactions within the vector host.
This study is a first step in improving our understanding of mosquito microbiota, elements of which (e.g., Wolbachia) are known to influence disease transmission. Further research is required to understand how the dynamics of mosquito microbiota across vectors and landscapes can affect mosquito-borne disease transmission, and ultimately the impact of land use and associated habitat change on infectious disease emergence.

ACKNOWLEDGMENTS
We would like to thank R. Rattanarithikul

DATA ACCESSIBILITY
Raw unprocessed sequencing reads are available through the NCBI Short Read Database as part of BioProject PRJNA406988, and the representative sequences of the bacterial and eukaryotic datasets are found in Tables S5 and S6, respectively.

CONFLICT OF INTEREST
None declared.

AUTHORS' CONTRIBUTIONS
PT participated in study design, mosquito identification, molecular laboratory experiment, and sequence analysis and carried out fieldwork. JAC participated in sequence analysis. ABG participated in study design and mosquito identification and carried out fieldwork.
PK participated in study design and coordinated and accommodated field experiments and parts of molecular laboratory processes. BAW, DDK, and SNB conceived of the study and participated in study design and coordination. PT, JAC, and SNB wrote the manuscript. DDK made substantial edits. All authors read, gave input, and approved the final manuscript.