Geographical separation and physiology drive differentiation of microbial communities of two discrete populations of the bat Leptonycteris yerbabuenae

Abstract In this paper, we explore how two discrete and geographically separated populations of the lesser long‐nosed bat (Leptonycteris yerbabuenae)—one in central and the other in the Pacific region of Mexico—differ in their fecal microbiota composition. Considering the microbiota–host as a unity, in which extrinsic (as food availability and geography) or intrinsic factors (as physiology) play an important role in the microbiota composition, we would expect differentiation in the microbiota of two geographically separated populations. The Amplicon Sequences Variants (ASVs) of the V4 region of the 16s rRNA gene from 68 individuals were analyzed using alpha and beta diversity metrics. We obtained a total of 11 566 (ASVs). The bacterial communities in the Central and Pacific populations had a diversity of 6,939 and 4,088 ASVs, respectively, sharing a core microbiota of 539 ASVs accounting for 75% of the relative abundance, suggesting stability over evolutionary time. The Weighted UniFrac metrics tested by a PERMANOVA showed that lactating and pregnant females had significant beta diversity differences in the two populations compared with other reproductive stages. This could be a consequence of the increased energy requirements of these physiological stages, more than the variation due to geographical separation. In contrast, a positive correlation of the observed ASVs of fecal microbiota with the observed ASVs of plastids related to the diet was observed in the juveniles and adults, suggesting that in these physiological stages an extrinsic factor as the diet shapes the microbiota composition. The results provide a baseline for future studies of the microbiome in these two wild populations of the lesser long‐nosed bat, the main pollinator of the Agaves from which the beverages tequila and mezcal are made.

and geography) or intrinsic factors (as physiology) play an important role in the microbiota composition, we would expect differentiation in the microbiota of two geographically separated populations. The Amplicon Sequences Variants (ASVs) of the V4 region of the 16s rRNA gene from 68 individuals were analyzed using alpha and beta diversity metrics. We obtained a total of 11 566 (ASVs). The bacterial communities in the Central and Pacific populations had a diversity of 6,939 and 4,088 ASVs, respectively, sharing a core microbiota of 539 ASVs accounting for 75% of the relative abundance, suggesting stability over evolutionary time. The Weighted UniFrac metrics tested by a PERMANOVA showed that lactating and pregnant females had significant beta diversity differences in the two populations compared with other reproductive stages. This could be a consequence of the increased energy requirements of these physiological stages, more than the variation due to geographical separation. In contrast, a positive correlation of the observed ASVs of fecal microbiota with the observed ASVs of plastids related to the diet was observed in the juveniles and adults, suggesting that in these physiological stages an extrinsic factor as the diet shapes the microbiota composition. The results provide a baseline for future studies of the microbiome in these two wild populations of the lesser long-nosed bat, the main pollinator of the Agaves from which the beverages tequila and mezcal are made.

K E Y W O R D S
geographical separation, holobiont, populations, reproductive stages

| INTRODUC TI ON
The microbial communities that inhabit the guts of mammals are complex, dynamic, and critical to the health of the host (Ingala, Becker, Bak Holm, Kristiansen, & Simmons, 2019;Ochman et al., 2010). Microbiota composition and abundance are further determined by several factors such as diet, geography, physiology and health status (Knutie, 2020;Ley et al., 2008;Muegge et al., 2011;Ochman et al., 2010), and even phylogeny. It has been demonstrated, for example, that phylogeny in hominids influences the microbiota at an evolutionary level, suggesting that the microbiota is more than just what the host eats (Ochman et al., 2010).
Microbial symbionts have a variety of roles in the nutrition, immunity, development, reproduction, and speciation of their eukaryote hosts, making this symbiosis a major component of eukaryotic fitness and evolution (Brucker & Bordenstein, 2012). Interactions between host and microbiome have coevolved in mutual adaptations into a superorganism association (MacColl, 2011) and are key to biological adaptations (Brockhurst & Koskella, 2013).
The hologenome concept of evolution states that in addition to the host genome, a significant proportion of the microbiome is also transmitted from one host generation to the next and can thus propagate unique properties of the holobiont (Rosenberg & Zilber-Rosenberg, 2019). Furthermore, genetic variation can occur by changes in the host and/or microbiome genomes; the microbiome genome is even expected to be able to adjust to environmental dynamics faster and through more processes that the host, and thus may be playing fundamental roles in the adaptation and evolution of holobionts (Rosenberg & Zilber-Rosenberg, 2019).
Genetic variation in holobionts can occur by mutation and DNA rearrangement, amplification or reduction of specific microbes, acquisition of novel microbes from the environment, and horizontal gene transfer from microbe to microbe or from microbe to host (Zilber-Rosenberg & Rosenberg, 2008). The multilayered structure of the holobiont consists of a core of host-adapted microbiota assembled from diverse environments and determined by genetic factors, as well as a flexible pool of microbes that depend on environmental diversity and external conditions (Shapira, 2016).
Vertical transmission of the microbiota from parent to progeny allows its maintenance between generations and is likely favored by selection when those microbes are beneficial to the host (Shapira, 2016). The functional profile of the microbiome is probably more constrained by evolution/vertical inheritance than individual bacterial taxa (Martiny et al., 2006;Phillips et al., 2017;Shapira, 2016).
The mutualistic and adaptive potential of microbiota offers flexibility and the ability to adapt to a range of ecological niches (Alberdi, Aizpurua, Bohmann, Zepeda-Mendoza, & Gilbert, 2016).
Speciation-the splitting of a population into two reproductively incompatible populations, each taking a distinct evolutionary path-is a defining process in evolution. Mutualist microbes can drive all modes of isolation, including ecological, behavioral, and developmental asynchronization and genetic divergence, but, as facilitators of niche adaptation, they are frequently associated with ecological isolation (Brucker & Bordenstein, 2012;Shapira, 2016).
Early evidence suggested that in the order Chiroptera, microbiome composition was influenced by the host phylogeny and life history (Phillips et al., 2012). In contrast, recent studies find no evidence of phylosymbiosis in Afrotropical bats (Lutz et al., 2019). Due to their species-specific feeding strategy specialization, Phyllostomid bats could remain as a model clade to study the relationship of the microbiome-host composition, phylogeny, and coevolution (Carrillo-Araujo et al., 2015). Phyllostomid bats are found from the southern USA and northern Mexico to Argentina and show a great evolutionary diversification of species, in which patterns are dependent on geographical and ecological interactions, resulting in a great diversity of dietary strategies and the most ecologically diverse family within the order Chiroptera (Carrillo-Araujo et al., 2015). They show a remarkable degree of evolutionary diversification of dietary strategies, from insectivory (the ancestral trait) (Monteiro & Nogueira, 2011) to feeding on blood, small vertebrates, nectar, fruit, and complex omnivorous diets (Gardner, 1976).
A study of the fecal microbiota of different reproductive stages within the central population of L. yerbabuenae showed that microbiota diversity is related to the reproductive stage rather than geographical distribution within this population (Gaona, Gómez-Acata, Cerqueda-García, Neri-Barrios, & Falcón, 2019). Microbiota composition is consistent in juveniles and nonreproductive females and males, regardless of the roost. Pregnant and lactating females' microbiotas were similar and more diverse than juveniles and nonreproductive adults. One explanation for this is that microbiota evolved with its host to be flexible enough to shift from a specialized diet to a more generalist diet to cope with the increased energy requirements during pregnancy and lactation (Gaona et al., 2019).
The aim of this research is to evaluate the composition of fecal microbiota in two discrete, geographically separated and genetically differentiated lesser long-nosed bat (L. yerbabuenae) populations in Mexico (Morales- Garza et al., 2007). Our hypothesis is that the geographical separation and reproductive isolation of the two L. yerbabuenae populations will lead to significant differences in their fecal microbial composition, making it more similar within populations than between them. We also evaluate whether the heterogeneity within the Pacific population is due to differences between the reproductive stages, as occurs in the central population (Gaona et al., 2019).

| Bat fecal microbiome sampling
Bats were captured using 12-m mist nets (Avinet) at the entrance of the caves using Kunz's technique (Kunz, Betke, Hristov, & Vonhof, 2009) between 18:30 and 7:00 hr. Standard measurements were taken in order to confirm identification and assess the age and reproductive stage of individuals (Anthony, 1988;Kunz, Wemmer, & Hayssen, 1996). All fecal samples were obtained directly from individuals and frozen in sterile Eppendorf tubes with liquid nitrogen Dewar and stored at a −20°C temperature until processed for DNA extraction (see Appendix 2).

| Analysis of the sequence data
The paired-end 2 × 250 reads were processed in QIIME2 (Bolyen et al., 2019). The reads were denoised with the DADA2 (Callahan et al., 2016) plugin to resolve the amplicon sequence variants (ASVs).
Both forward and reverse reads were truncated at 200 bp, and chimeric sequences were removed using the "consensus" method.
Representative ASV sequences were taxonomically assigned using the "classify-consensus-vsearch plugin" with default arguments, using the SILVA 132 database as a reference (Quast et al., 2013). An alignment was performed with the MAFFT algorithm (Katoh, 2002).
After masking positional conservations and gap filtering, a phylogeny was built with the FastTree2 algorithm (Price, Dehal, & Arkin, 2010 (Oksanen, 2015), and ggplot2 (Wilkinson, 2011) packages. Plastidic ASVs were filtered out of the samples (for subsequent separate analysis, see below), and then, the samples were rarefied to a minimum sequencing effort of 10,000.
A PCoA ordination was performed with the weighted UniFrac distance ( Figure 2). The Shannon alpha diversity index was calculated ( Figure 4). A PERMANOVA with a weighted UniFrac distance was performed to assess whether there were significant differences among groups of samples (region and reproductive stages as predictor variables) with 1,000 permutations (Table 2).
To explore the differential abundance of taxa, a LEfSe (linear discriminant analysis of effect size) analysis (Segata et al., 2011) was performed at the family level, first with all samples using the populations as categories, and then within each population, using the reproductive stages as categories, using an LDA cutoff >2 and p value < .05.
Counts of plastidic ASVs (separated from prokaryotic ASVs before rarefaction) were normalized with the cumulative sum scaling (CSS) method with the metagenomeSeq package (Paulson, Stine, Bravo, & Pop, 2013), and a DPCoA (Double Principle Coordinate Analysis) and the Shannon index were calculated to assess the likely variation in the diet of the bats in the different stages within the regions. A PERMANOVA was carried out with the DPCoA distance matrix to assess the beta diversity between populations. A Pearson correlation analysis was performed with the observed ASVs of plastid and prokaryotic 16S to assess the relationship between microbiota diversity (prokaryotic 16S ASVs) and diet diversity (plastidic ASVs).

| Microbiome composition in the two populations of Leptonycteris yerbabuenae
Of the total samples collected, 68 were positively PCR amplified: 30 for the Pacific population (12 juvenile, 7 adult, 6 pregnant, and 5 lactating) and 38 for the central population (6 juvenile, 20 adult, 6 pregnant, and 6 lactating). We obtained a total of 11 566 ASVs after rarefying to a sampling depth of 10,000 reads (see the accumulation curve in Figure A1 of Appendix 1).  There was a trend toward higher diversity among pregnant and lactating females in the central population, but higher diversity among juveniles in the Pacific population (Figures 3 and 4). The Shannon index shows that alpha diversity was more homogeneous within the Pacific population ( Figure 4) than in the central population (Gaona et al., 2019).
Comparison using a PERMANOVA test showed differences among reproductive stages; there were significant differences in pairwise comparisons of lactating and pregnant females versus juveniles and adults, and a significant effect of reproductive stage nested within the population (Table 2)

| Bat diet variability between the two regions
The amplification of plastid 16S genes is incidental in fecal microbiome studies. Since chloroplasts are plant organelles acquired by endosymbiosis, their amplification is due to contamination when the host feeds on plant material. Since different plant species' chloroplasts differ genetically, plastidic 16S gene diversity may provide a relative index of plant diversity in the diet (Knight et al., 2018). The number of plastidic reads in the 68 samples ranged from 17 to 76,979, with a mean of 5,589 (see Table A1). We used the entire plastid dataset with the rationale that although some samples had very low plastid reads, the total reads per sample including bacterial 16S had more than 10,000 reads. After normalization by the CCS method to avoid the effect of differences in the libraries' sample sizes, we calculated the Shannon index and performed a DPCoA and a PERMANOVA ( Figure 6 and Table 2).
Similar to the microbiota of the two populations, the beta diversity of chloroplast ASVs was higher in the central population TA B L E 2 PERMANOVA analysis. In the pairwise mode, p value was adjusted using the Bonferroni method for multiple comparisons. In the nested mode, stages was nested within the populations in both, for prokaryotes and plastids (variance within the population), but alpha diversity was similar ( Figure 6). Chloroplast diversity also differed among reproductive stages, with pregnant and lactating females having the highest beta chloroplast diversity ( Figure 6). . It remains to be tested whether these interpopulation microbiome differences could indicate differential success if the changes in microbial composition are related to losses or gains of metabolic capabilities that contribute to host fitness (Hooper, Littman, & Macpherson, 2012).

| D ISCUSS I ON
We found that the two populations share a core microbiota of only 539 ASVs and that these are the most abundant bacterial groups for both populations, accounting for about 75% of the relative abundance ( Figure 7). This suggests stability in the microbiome within the populations over evolutionary time, and consistent with findings in hominids, the microbiome is more related to the species than to diet (Ochman et al., 2010). Each population has six families with differential abundance detected by the LeFSe analysis ( Figure 5a) and a difference of 2,851 ASVs. In addition, ontological characteristics of bats, such as the rearrangement and growth of organs, including the stomach and the enlargement of the intestine, during pregnancy (Speakman, 2008), constrain differences and similarities in the microbiome (Gaona et al., 2019). Despite the changes in beta diversity due to geographical separation, the reproductive stages are more homogenize (Xiao et al., 2019). Our results partially support this result, since the populations we have studied are more similar within populations. However, we also found a core microbiome that also supports microbiota stability (Coyte, Schluter, & Foster, 2015) and ontogenetic differences in both populations (Gaona et al., 2019). The correlation between microbiota ASVs and plastid ASVs was positive F I G U R E 5 LeFSe analysis of samples with an LDA cutoff > 2 at the family level.
(a) Families with differential abundance between populations; (b) Venn diagram showing the ASV shared between populations; (c) Families with differential abundance detected by the LeFSe analysis between reproductive stages within each population-for the central population, just the 10 families with the highest LDA score are shown for each stage (full list provided in Appendix 1) tions. This suggests that in the demanding stages of pregnancy and lactation, the alpha diversity of microbiota increases more steeply than the plastids (diet), suggesting that the physiological state is driving the changes more than the diet. Thus, the core microbiota is present throughout the life cycle, but in pregnancy and lactation, the microbiota is more divergent (Figure 7). The PERMANOVA analysis of only pregnant and lactating females between the two populations showed that these stages have a higher coefficient of variance (R 2 ~ .13) than the adult and juvenile stages (R 2 ~ .10) ( Figures A4   and A5).

F I G U R E 6
The second and third requirements-differential success and evidence of inheritance of the microbiome-remain to be evaluated.
One of the main predictions of microbiome heritability is that the offspring's microbiome will resemble parental microbiome. Thus, a progenitor holobiont (species A with its microbiome) will resemble its holobiont offspring (species B with its microbiome). We are unable to directly evaluate that statement here since we do not have samples of a parental-progeny holobiont and this study solely focused on one bat species with geographical separation.
The term phylosymbiosis is used to describe the congruence between the differences in bacterial communities and the phylogenetic divergences among species (Brooks, Kohl, Brucker, Opstal, & Bordenstein, 2016;Kohl, Varner, Wilkening, & Dearing, 2018;Mazel et al., 2018), which, in a broad sense, corresponds to the inheritance of the microbiome. While the two populations we sampled in this study are the same species, our data show a relationship between differences in host bacterial communities and genetic differences among populations (Morales-Garza et al., 2007), which could suggest an early process of phylosymbiosis which has not been explored before in a geographical separation model. As described above, the differences between the two populations of L. yerbabuenae show higher variation (beta diversity) in the microbiota the central population (Figures 2 and 4), which also has a higher genetic diversity and is a more stable population (Morales-Garza et al., 2007).
Macroecological theory indicates that the highest genetic diversity will be found in areas with better access to resources and the best ecological conditions for the species. For this species, the available resources have determined two main divergence routes and very likely the segregation of the species, as confirmed by stable isotope analysis of mitochondrial DNA (Burke et al., 2019), and these changes correspond with orography (Burke et al., 2019). L. yerbabuenae are reliant on the availability of flowers, and the central region has year-long availability (Fleming & Nassar, 2002;Rojas-Martínez et al., 1999;Villaseñor et al., 2017) in contrast to the seasonal availability of the Pacific region (Petit, Excoffier, & Mayer, 1999). Thus, the reproductive peaks of bats and the population composition are different between the populations, apparently supported by the beta diversity of the plastid ASVs, and the variance between the two populations is high (R 2 ~ .2 from the PERMANOVA, Table 2).
Our results show a differentiation of the microbiota composition between populations and among reproductive stages. One explanation for these trends is the differentiation of the diet, which has been demonstrated to be the main factor shaping the functionality and diversity of the gut microbiome, resulting in the convergence between microbial communities and their hosts' foraging habits (Muegge et al., 2011). There is evidence that the diet shapes the functionality, diversity, and relative abundance of dominant phyla, as well as the populations of specific bacterial groups; this is influenced in part by F I G U R E 7 Family classification of the 539 shared ASVs among the two populations. The families with abundance below 1% were combined in the "others" category. The full list of families per sample is shown in Figure A1 the composition of macronutrients consumed and in part by the introduction of new microbes from food itself (Muegge et al., 2011;Voreades, Kozil, & Weir, 2014).
In the two populations of L. yerbabuenae, our results show significant differences in the beta diversity of lactating and pregnant females, perhaps associated with the high energy, protein and calcium requirements during these life stages (Gaona et al., 2019;Speakman, 2008). Previous work has shown that the diversity of plants consumed by L. yerbabueanae is 2.4 times higher among pregnant and lactating females than among nonreproductive females and males (Riechers, Martínez-Coronel, & Vidal, 2003). This increased foraging diversity is also reflected in our results, with a higher diversity of gut microbes in pregnant and lactating females in the central population, likely due to the diversified diet due to the physiological conditions of pregnancy and lactation and the high plant diversity in the central region ( Figure 6).
The analysis of chloroplasts showed a difference between the two populations; while there was some overlap between the central population and Pacific population adults, the two populations tended to separate, with lower diversity in the Pacific ( Figure 6).

| CON CLUS ION
The differential fecal microbiota composition of the two L. yerbabuenae populations, central and Pacific, which inhabit geographically disjunct ranges that differ in their availability of nectar and pollen, opens the door to exploring microbiome-bat relationships that may be influenced by natural selection. Both variation and stability are present, suggesting that it is not just the host, but the host-microbiome unit, the "holobiont," which could be subject to natural selection and evolutionary processes (Cerqueda-García & Falcón, 2016;Suzuki, 2017). This study suggests a differentiation between both populations' microbiota as a consequence of geographical separation and food resources availability. Diversity was highest in the central population, where there is also higher genetic diversity, perhaps due to the unrestricted availability of foraging resources. This is in contrast to the Pacific population, which had lower microbiome and genetic variability, probably because resources are only available seasonally. These resource availability interactions influence reproduction and population size, with increased genetic diversity, both in the host and in the microbiome. Interestingly, these two populations share a core microbiota of 539 ASVs accounting for 75% of the relative abundance, being the most abundant microbiota stable over evolutionary time, and these results may support the phylosymbiosis theory.

ACK N OWLED G M ENTS
This study comprises part of the PhD thesis work of O. Gaona

CO N FLI C T O F I NTE R E S T S
None declared. Resources; Supervision.

E TH I C S S TATEM ENT
Sampling procedures strictly followed the guidelines of the American Society of Mammalogists for capture, handling and care of mammals (Gardner, 1976;Sikes & Gannon, 2011). To minimize animal suffering and distress, each animal was processed by experts only, ensuring safe and effective handling. Samples were taken from wild bats that were released in the same area of capture, causing no apparent harm to individuals. Leptonycteris yerbabuenae is not subject to federal pro-