Wild specimens of sand fly phlebotomine Lutzomyia evansi, vector of leishmaniasis, show high abundance of Methylobacterium and natural carriage of Wolbachia and Cardinium types in the midgut microbiome

Phlebotomine sand flies are remarkable vectors of several etiologic agents (virus, bacterial, trypanosomatid Leishmania), posing a heavy health burden for human populations mainly located at developing countries. Their intestinal microbiota is involved in a wide range of biological and physiological processes, and could exclude or facilitate such transmission of pathogens. In this study, we investigated the Eubacterial microbiome from digestive tracts of Lu. evansi adults structure using 16S rRNA gene sequence amplicon high throughput sequencing (Illumina MiSeq) obtained from digestive tracts of Lu. evansi adults. The samples were collected at two locations with high incidence of the disease in humans: peri-urban and forest ecosystems from the department of Sucre, Colombia. 289,068 quality-filtered reads of V4 region of 16S rRNA gene were obtained and clustered into 1,762 operational taxonomic units (OTUs) with 97% similarity. Regarding eubacterial diversity, 14 bacterial phyla and 2 new candidate phyla were found to be consistently associated with the gut microbiome content. Proteobacteria, Firmicutes, and Bacteroidetes were the most abundant phyla in all the samples and the core microbiome was particularly dominated by Methylobacterium genus. Methylobacterium species, are known to have mutualistic relationships with some plants and are involved in shaping the microbial community in the phyllosphere. As a remarkable feature, OTUs classified as Wolbachia spp. were found abundant on peri-urban ecosystem samples, in adult male (OTUs n = 776) and unfed female (OTUs n = 324). Furthermore, our results provide evidence of OTUs classified as Cardinium endosymbiont in relative abundance, notably higher with respect to Wolbachia. The variation in insect gut microbiota may be determined by the environment as also for the type of feeding. Our findings increase the richness of the microbiota associated with Lu. evansi. In this study, OTUs of Methylobacterium found in Lu. evansi was higher in engorged females, suggesting that there are interactions between microbes from plant sources, blood nutrients and the parasites they transmit during the blood intake.

Adult specimens were collected from the two locations using Shannon-type extra-domiciliary white light traps that remained active between 18:00 h and 22:00 h 11 . Adult sand flies collected were transported live to the laboratory in entomological cages to obtain gut and morphological structures (male genitalia and female spermathecae) for taxonomic identification of Lu. evansi using a taxonomic key 11 . Guts of adult Lu. evansi specimens were removed aseptically with sterile stilettos under a stereoscope in 1X PBS buffer 11 . Gut pools were grouped according to locality and physiological stage; these were re-suspended in 100 μL of sterile PBS, before extracting total DNA 11 . The conformation of groups assessed (source, number of intestines, sex, feeding status) and DNA concentrations obtained for sequencing is illustrated in supplementary online material (Table SP1). Samples from each municipality were collected in one single sampling campaign 11 . The number of guts pooled on every group had on average 17.45 ± 2.98 specimens, a retrieval depending of availability in the natural populations at the time of sampling 11 . Bacterial 16S rRNA gene fragment PCR amplification and sequencing. PCR amplicon libraries of V4 region of 16S rRNA genes were prepared. Metagenomic DNA was extracted of gut pools from Lu. evansi males (Colosó = COM; Ovejas = OVM), fed females (Colosó = COFF; Ovejas = OVFF), and unfed females (Colosó = COUF; Ovejas = OVUF) from both locations 30 . Total gut DNA of Lu. evansi was extracted using the Ultra CleanTM Soil DNA Isolation Kit (MO BIO Laboratories, Inc., USA), according to manufacturer's instructions 30 . Final DNA aliquots were quantified using an ND-100 Nanodrop Thermo Scientific spectrophotometer (Thermo Fisher Scientific Inc, MA). The quality of genomic DNA was analyzed on 1% agarose gel, followed by EZ-visionTM DNA 6X STAIN (AMRESCO) 30 .
Intestinal microbiome analysis and statistics. All sequence analyses were performed using Quantitative Insight software package into Microbial Ecology (QIIME 1.9.0 version) 32 and Phyloseq R package 1.5.21 33 . Initially, the complete set of Illumina reads was filtered and split according to barcodes for each sample, chimeric sequences were identified, extracted and excluded out of datasets by Usearch 6.1. A multi-step open-reference OTU picking workflow was performed within QIIME system. Centroid sequence was extracted as a representative sequence for each OTU picked and aligned to SILVA 128 release. In the next steps, a single representative sequence of each OTU was again realigned using PyNAST to build a phylogenetic tree using FastTree 30 . Sequences were deposited in NCBI sequence reads archive (SubmissionID: SUB4855063; BioProject accession: PRJNA507409; BioSample accession: SAMN10492589; Accession: KCPM01000000). All OTUs classified as ascribed to ribosomal sequences from Archaea or from Eukaryotic mitochondria or plastids present and amplified in metagenomic DNA were eliminated and not used in subsequent analyses, as the purpose of this study was to focus on the composition of Eubacterial microbiome 30 .
To estimate species richness and diversity across communities, OTU based analysis was used and implemented using both software package QIIME and Phyloseq R package 32,33 . Based on the OTU assignment, this QIIME workflow script computed measurements of alpha diversity and generated rarefaction plots. Phyloseq R package was used to produce graphics representing abundances, diversity (including Shannon diversity index, observed species, chao1 and Simpson) and distribution of all phyla throughout each sample. Alpha diversity (Shannon_H) and Evenness (Simpson_1-D; Dominance_D; Berger-Parker; Equitability_J) metrics also were estimated in the PAST package, using OTUs associated with core microbiome, with the aim of obtaining a more specific analysis 30 .
Diversity between all the communities (beta diversity) was analyzed by generating Unifrac calculations in both methodologies weighted (presence/absence/abundance matrix) and unweighted (presence/absence matrix) version 30 . This parameter, besides measuring differences between two collections of sequences as the amount of evolutionary history unique to either one, accounts for differences in relative abundances 34 . Using Phyloseq, additional dendrograms were produced out of these calculations. These plots also allowed comparisons of differences and similarities between communities among all the samples 30 . Bray-Curtis dissimilarity matrix was calculated and used it to build Principal Coordinate Analyses and Constrained Principal Coordinate Analysis constrained using Phyloseq R package for your graphically visualized. Differences in the beta diversity of bacterial communities were verified using the non-parametric Permutational Multivariate Analysis of Variance (PERMANOVA) and nonparametric adonis tests, both retrieved from Phyloseq 30 .
Data were analyzed by one-way ANOVA for tested differences in these indices of alpha diversity between groups of guts established, followed by Levene's Test for Homogeneity of Variance, Shapiro-Wilk test for checking normality, and Durbin Watson test for checking independence 30 . A principal component analysis (PCA) was performed using the "XLSTAT-Ecology" package to identify patterns of bacterial taxa distribution between guts at different stages, sex, and among biotypes 30 .

Results
335,437 sequences of 250 nucleotides on average were obtained by the Illumina-MiSeq platform, covering region V4 of the 16S rRNA gene. Following the quality filter steps and removal of chimeric sequences, 289,068 reads were clustered into 1,762 OTUs at 97% similarity. Regarding Eubacterial diversity, OTUs recovered from all Lu. evansi samples were affiliated with 14 described bacteria phyla and 2 candidate phyla (Fig. 2a).
A rarefaction analysis by Chao1 estimator values revealed a good coverage of sequencing, suggesting that most of the members of all the guts associated community have a description near complete (Fig. 2b). Samples CO10UF, CO7UF and CO9FF were not reaching in this survey a coverage close to saturation; thus, they may still harbor a higher richness of non-detected OTUs. Through this analysis it was possible to determine significant differences in the richness between all samples, being OV5M and OV6UF those with the highest number of OTUs, doubling the number contained of almost all other samples (Figs 2a, 3).
Proteobacteria, Firmicutes and Bacteroidetes were the most abundant phyla in Lu. evansi. Proteobateria composed more than 65% of communities in most of the samples, with the exception of OV5M (30%) and OV3UF (50%), which also have the largest fraction of predominant OTUs (Figs 2a, 3). As a proxy to establish a comparable set of predominant OTUs in the communities associated to all specimens, we defined "abundant OTU" as the one with a number of sequences equal to or greater than 1% of the total number of sequences analyzed in each sample (Fig. 3). The resulted in a set of 46 OTUs classified inside five dominant phyla including Proteobacteria, Firmicutes, Bacteroidetes, Verrucomicrobia and Gemmatimonadetes (Fig. 2).
Distances between communities were measured using ( Fig. SP1) unweighted UniFrac method to see the difference based on lineages they contain but also on the sum of their counts. Comparison of weighted and unweighted ( Figure SP1, Fig. 3) Unifrac analysis supported what can be inferred by the Chao1 estimator; communities from samples OV5M, OV6UF and CO9FF have a significantly higher number of OTUs, but those predominant, have a very lower frequency in that group than in the communities associated with the other specimens. Dendrograms do not show clustering by gender, origin or feeding status, suggesting that these features are not strictly determining the microbiome structure in the set of samples analyzed.
Analyzing the set of shared OTUs between all communities, we defined the core community as the set of OTUs present in all specimens and collectively comprising a major fraction of the communities: We found 14 OTUs with these features that comprised a fraction between 30 to 70% in all communities, most of them belong Scientific RepoRtS | (2019) 9:17746 | https://doi.org/10.1038/s41598-019-53769-z www.nature.com/scientificreports www.nature.com/scientificreports/ to the Proteobacteria phylum. The detailed structure and taxonomy composition of core community are presented in Fig. 4. Three different OTUs from Methylobacterium organophilum are present in this set and one of them it is also the most abundant one (Fig. 4). OTUs from the families Oxalobacteraceae, Pseudomonadaceae, Infraporangiaceae and Staphylococcaceae are close in abundance within the core community (Fig. 4). Rarefaction curve from Chao1 analysis using partial 16S rRNA gene sequences of guts from field-collected Lu. evansi (male/fed females/unfed females). Saturated rarefaction curve indicates that the vastness of microbial diversity was retrieved from each sample. Nomenclature of gut pools: COM: males from Colosó; COFF: fed females from Coloso; COUF: Unfed females from Coloso; OVM: males from Ovejas; OVFF: fed females from Ovejas; OVUF: unfed females from Ovejas.
Frequency patterns of the most abundant OTUs revealed that there are three samples; OV6UF, OV5M and OV3UF with a greater number of predominant OTUs, wherein each one comprised a maximum 15% of the community. In all samples, the most abundant OTU was Methylobacterium organophilum (Fig. 4), ranging between 30 to 56% in all communities except for the three previously mentioned, where it did not exceed 15%.
A conclusive finding is the detection of OTUs of Wolbachia and Cardinium endosymbionts in guts from natural populations of Lu. evansi (Fig. 5). A high relative abundance of OTUs of Cardinium was found mainly in males   www.nature.com/scientificreports www.nature.com/scientificreports/ (OV5M OTUs = 1326) (Fig. 5), and these are related to samples from the Peri-urban (Ovejas) and forest (Coloso) ecosystems (Fig. 5), whereas OTUs of the endosymbiont Wolbachia were only identified in the Peri-urban environment and in the intestines of males and females (Fig. 5), in lower relative abundance with respect to Cardinium (Fig. 5). Both symbionts are present in unfed females.
According to the ANOVA, there are no statistically significant differences between groups of guts established from Lu. evansi according to the sex or nutritional status (one-way ANOVA, df = 2, F = 0.235, P < 0.796). A greater specific richness is distinguished in the peri-urban ecosystem with respect to the forest ecosystem (Fig. 6). The same behavior is related to Chao1 and Shannon diversity indices. The specific richness is higher in males with respect to females. Only the Chao1 index in unfed females is similar (Fig. 6).
Analysis of alpha diversity of the core gut microbiome of Lu. evansi was performed. In terms of species richness, which is directly associated with the number of OTUs per sample, a high number of OTUs was detected in fed females of the two biotypes (OTUs CO9FF = 13655, OTUs OV1FF = 13623) ( Table-SP2), as well as in males (OTUs C011M = 18941; OTUs OV2M = 14909). Dominance revealed a higher index in the C09FF sample, which reflects the relative importance of given OTUs (as Methylobacterium) or a group (Table S1, Fig. 3), in this case, showing the influence of blood intake in Lu. evansi midgut. Simpson diversity index is higher in males of the two biotypes (CO11M, C08M, OV5M) and in unfed females of the peri-urban biotype (Table S1), a greater diversity and changes in shared common species, possibly triggered by the exposure to several microorganism sources. Among the intestine samples of Lu. evansi, we observed a higher uniformity (e ^ H/S = 0.5104) of bacterial OTUs in the C08M group of males (Fig. 4). In all the samples, the equity index does not have values close to 1, suggesting that bacterial OTUs have the same or similar abundances, reflecting highly dynamic bacterial communities within Lu. evansi intestines. Berger-Parker index estimated a decrease in equity and an increase in dominance, which is congruent to the detection of a core gut microbiome of Lu. evansi. Fed females mainly from the forest environment represent Berger-Parker indices of 0.6853 to 0.8295 (Table S1), so they can have members of some OTUs with a representation of more than half of the individuals, being close to the value of 1.
The CAP and PcoA plot (Beta-diversity analysis) show two clustering patterns using Bray-Curtis dissimilarity matrix, which indicate clear differences in the microbiota of males and females mainly (Fig. 7). The statistical analysis (permutation multivariate analysis PERMANOVA and Adonis test) totally reflected PCoA plots results, indicating a significant divergence between groups for both of sex (P = 0.002; R2 = 0.45). The R2 values indicate that ∼45% of the variation in distances can be explained by this grouping.

Discussion
Lu. evansi is an important vector insect in the transmission cycle of visceral leishmaniasis in the Caribbean coast of Colombia, its control is a health priority. For this reason, a complete analysis of their intestinal bacterial microbiome was developed in this work, in order to find possible alternatives for its biological control. A detailed understanding about how pathogens interact with their vectors and the resident microbes can lead to the discovery of new tools to block disease transmission 35 . The first and crucial step in this approach is the identification of suitable commensal bacteria in the vector.
Previous reports of gut microbiome in female sandflies from South America have been performed applying classic techniques of molecular ecology such as DGGE and band sequencing of 16S rRNA gene fragments, and some bacterial types (Enterobacter, Bacillus, Serratia) are usually found in plants and soils, including some sequences of potential plant pathogens or species with entomopathogenic potential 8,11,35 . In this study, we investigated the Eubacterial microbiome from digestive tracts of Lu. evansi adults structure using 16S rRNA gene sequence amplicon high throughput sequencing (Illumina MiSeq).
Proteobacteria, Firmicutes, Actinobacteria and Bacteroidetes were the most abundant phyla in all Lu. evansi samples. This result is in good agreement with a previous study on the gut microbiota of 81 insects, showing that Proteobacteria and Firmicutes were the predominant phyla, comprising 57.4% and 21.7% of sequences, respectively 36,37 . It is also consistent with microbiota reported to be associated to different species of Lutzomyia and Phlebotomus genus, as well of mosquitoes (Aedes and Anopheles) although, Bacteroidetes and Actinobacteria phylum showed lower abundance with respect to our study 3,5,19,38 .
Gut microbiomes on Lu. evansi clearly show such abundant phyla shared across all the samples, but also a large fraction of phyla with quite distinct and variable composition even between similar sample groups. The dynamic variation in insect gut microbiota could be explained due to the influence of microhabitats, gut morphology www.nature.com/scientificreports www.nature.com/scientificreports/ and physicochemical conditions, such as pH, oxygen availability in the insect gut, as well as type of diet 39 . These diverse gut conditions may cause variation in host-specific gut microbiota in insects, influencing nutritional and immune system status as well as reproductive fitness 15 .
Regarding the presence of constant OTUs shared among all Lu. evansi samples, notably most dominant Methylobacterium, a common cosmopolitan inhabitant of phyllospheres of plants (endophytic, transporting and rhizosphere tissues), thus suggesting the constant feeding from plants in specimens collected 40,41 . Recently, the family Methylobacteriaceae was reported with low representation in (gravid and fed females) Lu. intermedia, in an endemic area of leishmaniasis in Brazil 5,38 . These findings together allowed us to suggest that Methylobacteriaceae have the capabilities to remain into the midgut during physiological processes carried out during the blood intake, proposal backed up by experimental evidence of the study published by Kelly et al., 2017, showing the continued presence or even increase in abundance of Methylobacteriaceae after the sixth day of feeding and infection with Leishmania infantum parasites in Lu. longipalpis 19 .
The relative abundance of Methylobacterium is significant in salivary glands and midgut of mosquitoes infected and non-infected with P. falciparum, which suggests that interactions occur between microbes and parasites 42 . The above indicates that Methylobacterium can be important in several physiological and metabolic processes in Lu. evansi, which suggests that interactions could occur with Leishmania parasite. More studies are needed to deepen in these interactions.
Studies on distribution of Wolbachia among its arthropod hosts are important both for better understanding why this bacterium is so common 43,44 , as well as for its potential use as a biological control agent. In America, Wolbachia infection has been observed previously in Lu. cruciata, Lu. shannoni, Lu. vespertilionis, Lu. trapidoi, Lu. longipalpis, Lu. whitmani, and Lu. intermedia 19,23,45 . One of the main objectives of this study was to confirm, in a new set of samples and with a higher resolution experimental approach the natural presence of Wolbachia in the microbiome of Lu. evansi specimens wild populations collected in Colombia (department of Sucre). Wolbachia endosymbiont was detected in adult male and unfed female specimens, distinctly found only in samples coming from Ovejas. A high proportion of Wolbachia was determined, which allows to confirm the preliminary results by Recently, in groups of females fed from Lu. intermedia, the family Rickettsiaceae was the most prevalent, with 46.7% of the group represented by Wolbachia 5 . This group of females midgut microbiomes was also analyzed by Illumina 16S rRNA gene amplicon sequencing. Wolbachia is the microbiota genus most enhanced by the presence of blood meal in the midgut 5 . This is in stark contrast with our study of Lu. evansi and for Kelly et al. 2017 about Lu. longipalpis, in which Wolbachia spp., were identified only in sucrose-fed flies and not in blood-fed or infected flies 19,46 .
Despite several reports about the distribution of intracellular symbionts that regulate host reproduction, such as Wolbachia, there are few studies addressing the broad distribution of other specific symbionts on insect hosts 47,48 . Among these symbionts, we would like to remark the case of Cardinium, which is reported for the first time for Lu. evansi, and not reported so far in any other insect vector in Colombia. It is showing a variability between similar groups (ex: Cardinium present in one group of OVM, but not in the other one). Wild populations of Lu. longipalpis from Brazil had been reported as hosts of this bacterial type 37 . 'Candidatus Cardinium' is a recently described bacterium from the Bacteroidetes group involved in diverse reproduction alterations of its arthropod hosts (including cytoplasmic incompatibility, parthenogenesis, and feminization) similar to Wolbachia 49 . Therefore, this finding could have potential and significance for biological control design. Cardinium is also vertically transmitted to progeny but rarely show co-speciation with the host. In sap-feeding insects, as phlebotomine sandflies, plant tissues have been proposed as alternative horizontal routes of interspecific transmission, but experimental evidence is limited 50 .
This remarks the potential of these intracellular bacteria to be studied and considered as possible controlling agents 51,52 ; and allows setting new platforms and perspectives to further expand the observations to correlate its presence and abundance in a larger sampled population where Leishmania content is also defined on each specimen. Given the recent reports by Kelly et al. (2017) and Louradour et al. (2017), that providing experimental evidence about the requirement of a gut microbiome linked to the survival of L. infantum in Lu. longipalpis and of L. major in P. duboscqi sand flies, respectively 19,53 , this kind of conditioning and controlling of parasite hosting, mediated by the insect gut microbiome, could also may be happening in wild populations of Lu. evansi. Further efforts are underway to continue exploring these possibilities.
The presence of OTUs of Oxalabacteraceae and Staphylococcus in Lu. evansi, is also very remarkable, mainly in males and unfed females. Oxalobacteraceae is a family within the order Burkholderiales in the subclass Betaproteobacteria that contains 13 genera. Possibly these OTUs were acquired by Lu. evansi during the larval stage and remained in the intestine until adult. These results correlate with the previous finding of Burkholderia cenocepacia in larvae 11 .
It is notable that Staphylococcus genera has been associated only to the intestinal microbiota of sandflies from Old World (Phlebotomus argentipes, Phlebotomus chiniensis, Phlebotomus papatasi) 6 , but with significant representation of OTUs, which is similar to the results obtained in Lu. evansi. Staphylococcus has a variable consistency between hosts, and is mainly related to plants and sap suggesting that they are acquired through sugar feeding 39 .
Microbiota can vary between sexes, feeding status, origin, and status of Leishmania infection. Sequencing by Illumina MiSeq platform, allowed us to define the composition of the intestinal microbiota of phlebotomines in the wild 54 . The results presented here open up a promising path to continue the studies on these and other vector species of Leishmania in aspects such as enzymatic and antimicrobial activities, influence of microbiota on the resistance to insecticides, and especially, the recalcitrant or favorable effects that specific microbiome components may exert on vector competence.

Data availability
All sequence data generated in this study have been deposited in the European Nucleotide Archive with the accession number (SubmissionID: SUB4855063; BioProject accession: PRJNA507409; BioSample accession: SAMN10492589; Accession: KCPM01000000).