The gut microbiome and metabolome of Himalayan Griffons (Gyps himalayensis): insights into the adaptation to carrion-feeding habits in avian scavengers
Himalayan Griffons (Gyps himalayensis), large scavenging raptors widely distributed in Qinghai-Tibetan Plateau, have evolved a remarkable ability to feed on carcasses without suffering any adverse effects. The gut microbiome plays an important role in animal physiological and pathological processes, and has also been found to play a health protective role in the vulture adaptation to scavenging. However, the microbial taxonomic diversity (including nonculturable and culturable microbes), functions, and metabolites related to Himalayan Griffons have not been fully explored.
Methods
In the present study, the 28 fecal samples of the Himalayan Griffons and 8 carrion samples were collected and sequenced using high-throughput 16S rRNA gene sequencing methods to analyze the composition and functional structures of the microbiomes. Twelve fecal samples of the Himalayan Griffons were analyzed using untargeted Liquid Chromatography Mass Spectroscopy (LC–MS) to identify metabolites. We used different culture conditions to grow Himalayan Griffons gut microbes. Inhibitory effects of gut beneficial bacteria on 5 common pathogenic bacteria were also tested using the Oxford cup method.
Results
According to the results of the culture-independent method, a high abundance of four major phyla in Himalayan Griffons were identified, including Fusobacteria, Firmicutes, Bacteroidetes, and Proteobacteria. The most abundant genera were Fusobacterium, followed by Clostridium_sensu_stricto_1, Cetobacterium, Epulopiscium, and Bacteroides. The predicted primary functional categories of the Himalayan Griffons' gut microbiome were associated with carbohydrate and amino acid metabolism, replication and repair, and membrane transport. LC–MS metabolomic analysis showed a total of 154 metabolites in all the fecal samples. Cultivation yielded 184 bacterial isolates with Escherichia coli, Enterococcus faecium, Enterococcus hirae, and Paeniclostridium sordellii as most common isolates. Moreover, 7 potential beneficial gut bacteria isolated showed certain inhibition to 5 common pathogenic bacteria.
Conclusions
Our findings broaden and deepen the understanding of Himalayan Griffons' gut microbiome, and highlighted the importance of gut microbiome-mediated adaptation to scavenging habits. In particular, our results highlighted the protective role of gut beneficial bacteria in the Himalayan Griffons against pathogenic bacteria that appear in rotten food resources.
Located in southeastern Mexico, the Isthmus of Tehuantepec (IT) has been long recognized as one of the most important geographic features in the diversification of a large number of taxa in northern Mesoamerica (Binford, 1989; Peterson et al., 1999; Barber and Klicka, 2010; Ornelas et al., 2013). Studies in montane-forest birds distributed east and west of the IT have suggested its crucial role as a geographic barrier acting in different times or “pulses” (Barber and Klicka, 2010; Ornelas et al., 2013; Tsai et al., 2019). However, the role that IT has played in the diversification of lowland birds is still poorly understood (see Sosa-López et al., 2013), probably because it has been long assumed that lowland bird species can disperse freely throughout the region (Edwards and Lea, 1955; Binford, 1989), thus preventing complete divergence in these taxa (e.g., Howell and Webb, 1995).
The White-tipped Dove (Leptotila verreauxi) is a widely distributed polytypic species mainly distributed in tropical lowlands from extreme southern USA (the Lower Rio Grande Valley in southern Texas) through Mexico and Central America to northern Argentina (Friedmann et al., 1950). Currently, 13 subspecies are recognized (Gill et al., 2021), four of which are distributed in Mexico (Fig. 1): capitalis, endemic to the Tres Marias Islands; angelica, from southern Texas in the Gulf slope, and southern Sonora in the Pacific slope south to the IT; fulviventris, from the IT on the Gulf slope to the Yucatan Peninsula; and bangsi, from the southern IT and Chiapas (Friedmann et al., 1950) to Honduras and western Nicaragua (Hogan, 1999; Gill et al., 2021). The ranges of the three Mexican mainland L. verreauxi subspecies abut at the IT (Friedmann et al., 1950), although precise boundaries between these are unclear, suggesting a complex biogeographic history in the region. As currently described, the main differentiation occurs in a west-east axis (angelica – bangsi, fulviventris). Such west to east differentiation, apparently follows a geographic pattern opposite to what may be expected given the environmental gradient. A previous study on the biogeography of New World doves supported a similar pattern of genetic differentiation, however details on geographical/ecological limits for this pattern are unclear (Johnson and Weckstein, 2011).
Figure
1.
Geographical representation of the samples used in this study. Color of the circles correspond to detected clades. Green, NA clade; Blue, NCA clade; Yellow, YUC clade. Red and pink: South American clades. (A) Inset 1 shows the distribution of the haplogroups and the number of individuals in each locality in the Isthmus of Tehuantepec (IT). The black star (subdivided in the two-colored pie charts) at the southern IT depicts the secondary contact zone between the NA and NCA clades. Colors depict vegetation types at the IT. Dashed line shows putative boundaries between currently recognized subspecies: A = L. v. angelica, B = L. v. bangsi, C = L. v. fulviventris. (B) Inset 2 depicts the Tres Marías Islands, where the endemic L. v. capitalis is found. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
Here, we described genetic differentiation patterns at the IT using L. verreauxi as a model. Based on the distribution of three subspecies, which abut at the region (Fig. 1), as a starting framework, we expected: a) that genetic differentiation corresponds to the ecological structure, in which main differentiation would occur in a north to south axis across the environmental gradient. A west-east pattern of differentiation would not be expected due to the ecological continuity of the lowland forests in the region which has promoted the recognition of mainly phenotypically-based bird subspecies across a north to south ecological gradient (e.g., Howell and Webb, 1995); b) ecological niche differentiation, which would support the recognized differentiation in current taxonomic treatments; and c) that contact zones between differentiated taxa will be maintained by local environmental differences across the ecological gradient. We performed phylogenetic and population genetics analyses using two mitochondrial and one nuclear markers, and estimated divergence times between lineages to explore genetic variation and historical population dynamics, focusing particularly on the populations at the IT. In addition, we used correlative ecological niche models (ENM) to estimate paleodistributions of lineages and test hypotheses of niche conservatism, as well as evaluate the role of ecological differences as a biogeographic limit in these populations. All of these tests allowed us to establish if these distribution patterns may be the result of secondary contact, which would suggest a vicariant event at the IT. We also evaluated if there are differences in the ecological niche that could be contributing to the maintenance of differentiation.
2.
Materials and methods
2.1
Taxon sampling and laboratory protocols
We obtained 48 tissue samples from the collection of the Museum of Zoology “Alfonso L. Herrera” (MZFC), representing all Mexican subspecies (Appendix Table S1). We supplemented our samples with 20 sequences available in GenBank (Benson et al., 2013; http://www.ncbi.nlm.nih.gov; Appendix Table S1). Outgroups were selected based on a previous phylogenetic hypothesis (Johnson and Weckstein, 2011) and included six species: L. jamaicensis, L. cassini, L. megalura, L. rufaxilla, L. plumbeiceps and Leptotrygon veraguensis. In total, we included 68 individual genetic sequences, considering samples deposited in GenBank.
DNA was extracted using both a Qiagen DNeasy Blood and Tissue Kit following the manufacturer’s protocol, as well as a phenol-chloroform extraction protocol (Sambrook and Russell, 2006). We amplified and sequenced two mitochondrial DNA markers [NADH, dehydrogenase subunit 2 (ND2) and the cytochrome C oxidase subunit 1 (COI)]; and one nuclear marker [the β-fibrinogen intron 7 (β-fibint 7)] via the polymerase chain reaction (PCR) following standardized protocols (Lorenz, 2012). We used primers L5215-H5766 and 487L-H6313 (Sorenson et al., 1999) for ND2, L6626-H7005 (Hafner et al., 1994) for COI, and B17L-B17U (Prychitko and Moore, 1997) for β-fibint 7. These molecular markers provide good resolution when analyzing phylogenetic relationships and biogeographic patterns in Columbiformes, in addition to show high phylogenetic consistency between mtDNA and nDNA in some studies at the supraspecific level (Johnson and Clayton, 2000; Pereira et al., 2007; Johnson and Weckstein, 2011; Sweet and Johnson, 2015; Sweet et al., 2017; Sangster et al., 2018). PCR products were visualized in agarose gels (1.3%) stained with ethidium bromide (0.5 μL/mL). Purification and sequencing processes were performed at Macrogen Inc., South Korea. Chromatograms were edited using Geneious Prime 2019.0.4 (Kearse et al., 2012) and multiple alignments were performed using MAFFT (Katoh and Standley, 2013) with default parameters.
2.2
Model selection and phylogenetic analyses
Partition schemes and sequence evolution models for Bayesian Inference were selected using PartitionFinder v2.1.1 (Lanfear et al., 2017). Molecular evolution models selected for each of the partition schemes are shown in Appendix Table S2. Phylogenies were generated using Maximum Likelihood (ML) and Bayesian Inference (BI). For ML, we used RAxMLGUI (Silvestro and Michalak, 2012), with the selected partition scheme, the GTR substitution model and the algorithm “ML + rapid bootstrap” (Stamatakis et al., 2008), and 10,000 iterations. We used MrBayes for estimating BI trees, which were obtained using the partition schemes and selected models as previously described (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003). We ran 10,000,000 generations in four parallel runs and eight Markov Chain Monte Carlo (MCMC) sampling every 1000 generations. We used Tracer 1.7 (Rambaut et al., 2018) to determine a stationary phase and the burn-in, and used FigTree v1.4.3 (Rambaut, 2009) to visualize the consensus tree.
2.3
Divergence time estimation
Divergence times were estimated in BEAST v.1.10.4 (Suchard et al., 2018) using the partition schemes and substitution models previously selected. We applied a strict molecular clock with a Yule speciation process (Steel and McKenzie, 2001). Since we were analyzing genetic patterns at intraspecific level, we opted for a strict molecular clock, which assumes a single rate of constant evolution throughout the phylogeny (Pybus, 2006). We applied a divergence rate of 1.96% (0.0098 ± 0.0005 substitutions/site/Myr) for the mitochondrial markers (Sweet and Johnson, 2015), and a divergence rate of 0.53% (0.00265 ± 0.0005 substitutions/site/Myr) for β-fibint 7 (Johnson and Clayton, 2000). Two independent runs of 20,000,000 generations were generated starting from different random points, corroborating the −lnL values and values of the effective sample sizes (ESSs) > 200 for all of the parameters in Tracer 1.7 (Rambaut et al., 2018). We used LogCombiner for combining the output files and TreeAnnotator (Drummond and Rambaut, 2015) for obtaining a maximum credibility tree (MCCT) with a burn-in of 25% and a posterior probability limit of 0.5%.
2.4
Phylogeography
As a measure of genetic differentiation between populations, we calculated genetic fixation from the FST estimator, and interpreted its values based on Hartl and Clark (1997). From this index, the rate of gene flow between populations or Nm can be estimated, the higher this value the higher the gene flow (Nei, 1987). Genetic distances were measured using the Nei corrected distance parameter DXY (Nei, 1987). Populations were identified according to the phylogeny and haplotype networks. Historical demography analyses for each clade were calculated using the R2 test, which has been shown to have a superior statistical power in small sample sizes (Ramos-Onsins and Rozas, 2002). The significance of the R2 test was estimated through 1000 coalescent simulations for each clade. All phylogeographic tests and genetic diversity estimates were obtained in DnaSP v.6.12.01 (Rozas et al., 2017). Given that ND2 and COI have different sample size numbers, we built different haplotype networks in Network 10.1 using a Median Joining Network algorithm (Bandelt et al., 1999). We assessed genetic variation in L. verreauxi via an Analysis of Molecular Variance (AMOVA) in Arlequin 3.5.2.2 (Excoffier et al., 2005). Significance was estimated via two independent analyses with 16,000 permutations based on pairwise differences: one included all populations, while another included only populations present in the IT, to test whether the levels of variation were consistent when only these populations are considered.
2.5
Ecological niche models
Since Pleistocene climatic changes have significantly influenced the evolution and biogeography of lowland taxa in the IT (Castillo-Chora et al., 2021; Sánchez-González et al., 2023), we estimated paleodistribution models for populations of L. verreauxi during the Middle Holocene (MH, ~6000 years ago), Last Glacial Maximum (LGM, ~21,000 years ago), and Last Interglacial (LIG, ~120,000–140,000 years ago), as well as for the Present. For model building, we used the data corresponding to occurrence records of L. verreauxi that we used in molecular analyses, supplemented with those downloaded from the Global Biodiversity Information Facility (GBIF; https://doi.org/10.15468/dl.f7wwag). Duplicated and erroneously georeferenced records were filtered and eliminated for analyses, and we performed a randomized spatial thinning with the R package “spThin” (Aiello-Lammens et al., 2015), which allowed us to build highly accurate models by reducing biases associated to data concentration in some geographic areas. Based on the genetic structure, we divided the total data into three datasets corresponding to each lineage present in the IT. Given that geographic boundaries for each lineage in this lowland taxon are not clear, we generated polygons from occurrence data corresponding to unequivocally assigned genetic samples; this procedure allowed us to reduce identification bias for records from biological collections identified below the species level (Appendix Fig. S1).
We downloaded environmental variables for the present and for past scenarios corresponding to the MH, LGM, and LIG from WorldClim 1.4 (Hijmans et al., 2005) with a resolution of 2.5 min. For the MH and the LGM models, we used layers from three different Atmosphere-Ocean General Circulation Models: CCSM4 (Community Climate System Model) (Gent et al., 2011), MIROC-ESM (Model for Interdisciplinary Research on Climate) (Hasumi and Emori, 2004), and MPI-ESM-P (Max Planck Institute) (Baehr et al., 2015); layers corresponding to the present and the LIG were downloaded from the CCSM4 model. We retained environmental variables with a greater contribution to the present distribution ranges, while variables with high collinearity were eliminated via the variance inflation factor (VIF), with a VIF >10 as implemented in the R package “corrplot” (Wei and Simko, 2021). We therefore selected nine bioclimatic variables: mean diurnal range (BIO2), isothermality (BIO3), max temperature of warmest month (BIO5), mean temperature of wettest quarter (BIO8), precipitation of wettest month (BIO13), precipitation of driest month (BIO14), precipitation seasonality (BIO15), precipitation of warmest quarter (BIO18) and precipitation of coldest quarter (BIO19).
We used the “Grinnell” v.0.0.21 R package (Machado-Stredel et al., 2021) to define an accessibility area (M) for each lineage (Appendix Fig. S2) based on estimates of the historical accessible area derived from simulations considering biological processes such as dispersal, colonization, and information on climate change over time. These M areas were further refined through comparisons with the terrestrial ecoregions (Dinerstein et al., 2017), and biogeographic provinces of the Neotropics layers (Morrone, 2014). All environmental models corresponding to the present and projections to the past were generated using the maximum entropy algorithm in Maxent v3.4.4 (Phillips et al., 2017) as implemented in “KUENM” v 1.1.9 package (Cobos et al., 2019), that is designed to strengthen the model calibration and creation process, via the generation of a set of candidate models to select one or several optimal models from certain parameterizations. Finally, to explore potential regions with low reliability for model transference to the past, we performed extrapolation risk analyses via the Mobility-oriented parity (MOP) test in “NicheToolbox” (Osorio-Olvera et al., 2020). The final models were edited in QGIS v3.18 (QGIS Development Team, 2022; http://www.qgis.org/)
2.6
Niche conservatism and range-breaking tests
Due to the absence of conspicuous geographic barriers across the lowland distributional range of L. verreauxi, we explored if genetic structure has been maintained through ecological niche differentiation in the different lineages detected in Mesoamerica. We compared the ecological niches through paired tests of niche conservatism between each lineage; additionally, we tested whether environmental differences across the IT may act as an ecological barrier separating these lineages. All these analyses were performed using “ENMTools” (Warren et al., 2021).
Niche conservatism models were generated with Maxent models (Phillips et al., 2017) to perform pairwise comparisons between the observed lineages. These analyses compare the overlapping of empirical models with a null distribution constructed from a set of pseudo-replicates using Monte Carlo methods. Empirical and pseudoreplicated values were calculated from measures of niche overlap (Schoener’s D) via 100 replicates to generate the null distribution across all tests. We test whether niches differ from what is expected by chance via a niche equivalence test (or identity test) and the symmetric background similarity test (Warren et al., 2021). Niche equivalence and similarity tests were conducted using model-based techniques test as implemented in ENMTools (Warren et al., 2021). To assess the role of environmental differences as a barrier for dispersal of the lineages across the IT, we used blob and ribbon range-breaking tests (Glor and Warren, 2011). Blob range-breaking test evaluates whether regions occupied by each lineage are more environmentally different than expected by chance; if so, the environmental differences would represent an abrupt biogeographical barrier. On the other hand, ribbon range-breaking test assesses whether unsuitable habitat divides areas of suitable habitat for each of the lineages which may function as a biogeographical boundary (Warren et al., 2021). We selected two areas that may functioning as a barrier for dispersal (the “ribbon”) in this lowland taxon: 1) the Sierra del Tolistoque, in the southern IT, and 2) the wetland region in the southern Gulf of Mexico. These two areas may serve as barriers due to apparent environmental differences on each side, maintaining lineage divergence despite confirmed overlap between the NA and NCA lineages in the southern IT. To optimize the ribbon range-breaking test, we generated a set of random points within the two selected areas due to the small number of occurrence records.
3.
Results
We obtained 48 sequences corresponding to ND2 (1041 base pairs, bp), 26 to COI (379 bp) and 16 to β-fibint 7 sequences (1136 bp). mtDNA topologies obtained from BI and ML analyses were highly similar. The tree showed two highly supported monophyletic groups for L. verreauxi (Fig. 1), separated by the IT. One group included the North American (NA) clade, distributed west of the IT, and a few from northwestern Chiapas (geographically east of the IT); the other main group, distributed east of the IT, is subdivided in two groups, each of which present further well-supported groupings. Samples from the Soconusco region in the Pacific slope in southeastern Mexico were grouped in a Northern Central American Clade (NCA), which is sister to the Yucatán Peninsula clade (YUC), while samples from Southern Central America and northern South America (SCA) formed a clade sister to west and central South America samples (SA).
Our mtDNA phylogeographic analyses showed two sympatric lineages in the southern IT region (Fig. 1). The NA lineage is distributed across the region in both the Pacific and Gulf slopes, while the NCA lineage is apparently restricted to the dry forests south of the IT in the Pacific slope. The nDNA tree (Appendix Fig. S3) did not show a clear geographic structure. Most of the samples of Leptotila species used as outgroup were included in a group sister to several L. verreauxi samples. The other group included most of samples distributed west of the IT (angelica and capitalis), two samples from NCA (bangsi), and an outgroup species, L. jamaicensis. These groups shared indels that coincide in position and length, with low nodal support (PP > 0.76 and BS > 70).
3.1
Divergence times
Our data suggest a divergence event during the Late Pliocene – Lower Pleistocene transition at approximately 2.5 Mya [1.94–3.43 Mya, 95% Height Posterior Density (HPD)] between the group west of the IT (the NA clade) and the group east of the IT. Further divergence events in the group east of the IT occurred between the clades SCA–SA and the clades NCA and YUC during the Lower Pleistocene around 1.89 Mya (1.39–2.49, HPD). The most recent divergences occurred between NCA and Yucatan at 0.9 Mya (0.62–1.34, HPD) and between SCA and SA at 0.92–1.78 Mya. A recent differentiation event between populations of the Tres Marías islands (L. v. capitalis) and rest of NA clade occurred around 100 kya (0.04–0.18, HDP).
3.2
Phylogeography
The haplotype network for the mDNA markers showed the same structure as the phylogeographic tree (Fig. 2). We found 25 haplotypes in ND2 and 8 haplotypes in COI. In both haplotype networks, we found a relatively large number of mutational steps separating two main clades roughly divided at the IT (25 in ND2 and 11 in COI), while NCA and Yucatán clades are separated by a smaller number of mutational steps (14 in ND2 and 5 in COI). We detected a star-like pattern in the NA clade; shapes are unclear for the other clades. The demographic R2 test showed population expansion for the NA clade (P < 0.05); however, R2 values for NCA and YUC were not statistically significant (P > 0.05 for both).
Figure
2.
Phylogenetics of Leptotila verreauxi. (A) mtDNA haplotype network. Black dots represent mutational steps, while circles indicate possible extinct or unsampled haplotypes. Circle size represent the frequency of haplotypes. (B) mtDNA tree. Numbers above nodes correspond to Bootstrap/posterior probability values; numbers below nodes indicate 95% HPD (High Posterior Density) ranges for divergence dates (in Mya). The dotted black box groups samples from the Tres Marías Islands. Red labels in the terminals correspond to sympatric samples from the Sierra del Tolistoque in the southern IT. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
The FST, Nm, and DXY indices showed a strong differentiation in most of the populations (Table 1). All the FST values are close to 1.0, signifying high genetic fixation and nearly complete haplotypic separation, in all lineages, as expected in genetically structured populations. Nm in the other hand, indicated less than one migrant every four generations. Genetic distances, as measured by the Nei corrected distance parameter DXY revealed high levels of genetic differentiation (>2%) for all paired comparisons, except for NCA-YUC with 0.00748 (0.7%).
Table
1.FST values (below the diagonal) and DXY values (above the diagonal) for the different L. verreauxi lineages. Values in the parentheses represent the number of migrants (Nm) per generation.
The AMOVA analysis (Table 2) including all populations showed that 95.9% of the genetic variation is found between populations, while only 4.06% is found within populations. The results are similar when the NA-NCA lineages are included.
Table
2.
Results of AMOVA of the L. verreauxi populations.
Distribution models for the present and projections to the past showed high performance in each of the lineages (Table 3), with models that are significantly better than expected by chance (AUC ratio >1) and p-values <0.001. For the NA clade, environmental variables with highest contribution were precipitation seasonality and precipitation of wettest month (BIO 15 = 28.1%, BIO 13 = 22.3%, respectively); while for NCA were isothermality and mean diurnal range (BIO 3 = 48.3%, BIO 2 = 30.4%, respectively). In the case of YUC, precipitation of wettest month and precipitation of warmest quarter (BIO 13 = 48.4%, BIO 18 = 8.9%, respectively).
Table
3.
Summary of model performance statistics per lineage.
The areas of high environmental suitability for the three Mesoamerican lineages showed clear geographic changes in comparison to present (Fig. 3). During the LIG, NCA and YUC showed a smaller geographic distribution in small and isolated habitat patches. NA showed an increase in the geographic extension of the environmental suitability in southeastern Mexico and northern Central America. During the LGM, NA experienced a reduction in suitable areas, being restricted to three main lowland areas in both slopes and an almost complete disappearance of suitability areas south of the IT. Regions such as YUC exhibited a decrease in environmental suitability, while NCA lineage showed an increase of the suitable conditions towards the MH. It is important to observe that all the three lineages experienced an increase in the environmental suitable areas; these areas have been apparently maintained towards the Present. According to environmental suitability maps, current contact zones across IT may have not been formed until recently (at least since LGM), as most of suitable areas appear as allopatric through time intervals for which ecological modeling was conducted.
Figure
3.
Ecological niche models of each lineage during the Present, Middle Holocene, Last Glacial Maximum and Last Interglacial. Dark areas in each map in the Present correspond to areas of climatic stability.
Because results of D were similar and I tend to express greater variability, for niche comparisons we focused on the Schoener’s D. The niche identity tests revealed that the niches of each of the three lineages are not identical (i.e, equivalent; empirical Schoener D values lower than the pseudoreplicated values; P < 0.01; Table 4; Appendix Fig. S4). On the other hand, the background similarity tests did not reject the null hypothesis for the NA-NCA comparison, while in the other two comparisons, the null hypothesis was rejected, indicating niche divergence (i.e., the ecological niches of these clades are less similar than expected based on their background environments).
Table
4.
Schoener’s D values of the niche comparison tests (Identity, Similarity and Blob rangebreak).
Lineages
Identity test
Similarity test
Blob rangebreak
NA–NCA
0.297**
0.271ns
0.282ns
NA–YUC
0.169**
0.148*
0.159ns
YUC–NCA
0.090**
0.103*
0.099ns
The significance of the test is shown ns, non-significant; P < 0.05; **P < 0.01.
Range-breaking blob tests did not show empirical values with a significant deviation from the null hypothesis (P > 0.05; Table 4; Appendix Fig. S4), suggesting that no biogeographical barrier is associated with environmental differences across the IT ecological gradient. However, the ribbon range-breaking test for the NA and NCA lineages suggested the presence of a region with less suitable conditions for the NCA lineage in the contact zone (D value, P < 0.01 for the comparison of the NCA lineage vs. the ribbon; Table 5; Appendix Fig. S5). Similarly, the comparison between the NA and YUC lineages suggests the presence of a region of less suitable conditions for the YUC lineage (D value, P < 0.05 for YUC vs. the ribbon; Table 5; Appendix Fig. S6).
Table
5.
Schoener’s D values from the ribbon-rangebreak tests.
L1
L2
L1 vs. L2
L1 vs. Ribbon
L2 vs. Ribbon
Outside vs Ribbon
NA
NCA
0.297ns
0.184ns
0.049**
0.154**
NA
YUC
0.086**
0.253ns
0.263*
0.326ns
The significance of the test is shown ns, non-significant; *, P < 0.05; **, P < 0.01.
Our study illustrates both emergent patterns found in different taxa across the study region, as well as discrepancies between phenotypically-based and molecular-based taxonomies in L. verreauxi.
4.1
Phylogeography
Our phylogenetic tree (Fig. 2) agrees with divergence events described in other studies along the Mexican Pacific slope and the IT in which previously unexpected evolutionary patterns are being revealed (Arcangeli et al., 2018; Castillo-Chora et al., 2021; Sánchez-González et al., 2021, 2023; Llanes-Quevedo et al., 2022). On the other hand, our analyses also show that traditional phenotypically-based taxonomy (Fig. 1) may either fail to correspond with evolutionary patterns as differentiation may lie further east from the IT, perhaps up to the Usumacinta River delta, as in other taxa (Howell and Webb, 1995; Bradley et al., 2008; Ramírez-Díaz et al., 2023); or partially support them, as differentiation between NA and NCA corresponds to the ecological gradient, but for which sympatry in the southern IT (southern Oaxaca) between these same lineages may be (theoretically) unexpected.
Our phylogeographic tree showed further differentiation events within each group divided by the IT (Fig. 1). In the group west of the IT, within the NA clade, all individuals from Tres Marías islands (currently in subspecies capitalis) were grouped, showing a genetic structure consistent with simple paraphyly (Omland et al., 2006), in which lineage sorting seems to be incomplete. Simple paraphyly in this case may be the result of a small number of founders that likely dispersed to the islands (180–40 kya), which is consistent with islands emergence (McCloy et al., 1988) and with dispersal times in other insular bird species (Ortiz-Ramírez et al., 2018). In the group east of IT, genetic differentiation between sister clades YUC and NCA is relatively low (0.7%), also suggesting a recent divergence event. Genetic divergence between east and west Mesoamerican clades in L. verreauxi is very deep, underlying the role of the IT as a significant barrier for differentiation in lowland taxa. Values above traditional species thresholds for species recognition (Seutin et al., 1995; Hebert et al., 2004; Fuchs et al., 2021) were found between NA–NCA (3.8%), as well as between NA and YUC (3.6%). Our analyses also found deep differentiation in the SCA and the SA despite having a small number of samples.
4.2
Biogeography
Previous research in mid-sized New World doves (Leptotila, Geotrygon and Zenaida) suggested that Leptotila dispersed to North America from South America after the closure of the Isthmus of Panama about 3.4 Mya (1.9–5.8 Mya). According to Johnson and Weckstein (2011), this dispersal event allowed for both, the evolution of a northern clade (including all L. verreauxi taxa), and its radiation following a stepping-stone fashion back into South America. However, our results (Fig. 2) suggest two main differentiation events: the first involves North America and Central America + South America, while in the second one the earliest differentiation occurred in South America, thus suggesting that the Central American taxa (NCA and YUC) are the most recent diverged lineages. Our differentiation scenario differs from the stepping-stone biogeographic scenario suggested in Johnson and Weckstein (2011).
Based on our results, lineage differentiation in Mesoamerican L. verreauxi occurred during the Late Pliocene-Pleistocene (Fig. 2) at about 2.6 Mya (1.9–3.4 Mya) with the split of lineages west (NA) and east of the IT (NCA-YUC). This divergence event is consistent with climatic changes that led to vegetation changes across the region, likely severing connections between formerly continuous ecosystems and therefore, allowing persistence of suitable conditions within small patches (Fig. 3, see also Castillo-Chora et al., 2021), which reduced gene flow between lowland populations (e.g., Sánchez-González et al., 2021, 2023). Both The star-like pattern in the haplotype network and R2 statistic for the NA clade support the biogeographic pattern described above, suggesting recent population growth probably due to the expansion of suitable conditions across the region. Demographic patterns in the other clades are unclear from our data.
Divergence date for the west and east groups is also consistent with proposed geographic processes in the region, such as formation of a seaway crossing the IT at approximately 2–3 Mya during the late Pliocene (Barrier et al., 1998), which has been frequently invoked as the main promoter of vicariant events in the region for the differentiation in unrelated lowland taxa, such as mammals (Whitmore and Stewart, 1965), geckos (Butler et al., 2023), and montane birds (Ornelas et al., 2013). The seaway hypothesis, however, has been challenged (Wyatt Durham et al., 1955; Mulcahy et al., 2006), and alternative hypothesis have been proposed, such as the ecological niche divergence (see Sánchez-González et al., 2023). Although the hypothesis related to a seaway across the IT may explain west-east divergence pattern, further differentiation events across the north to south ecological gradient have been described in other vertebrates, such as reptiles (Parkinson et al., 2000; Zarza et al., 2008; Gray et al., 2019; Butler et al., 2023).
4.3
Secondary contact
Our analyses revealed a secondary contact zone for the NA and NCA lineages south of the IT (Fig. 1), in which genetically divergent individuals are syntopic in the seasonally dry forests foothills of the Pacific slope of the Sierra del Tolistoque, Oaxaca. Sympatric lineages with a deep genetic divergence, such as these in L. verreauxi, have been rarely reported in birds (Hogner et al., 2012; Benham and Cheviron, 2019; Pons et al., 2019), but also in other groups (Tominaga et al., 2009; Xiao et al., 2012).
Different causes may result in the establishment of a secondary contact zone (Benham and Cheviron, 2019). Sympatric overlap of distinct lineages is typically interpreted as the result of formerly isolated divergent populations that resumed secondary contact (Hogner et al., 2012; Benham and Cheviron, 2019). Pleistocene climatic fluctuations may have promoted the divergence of the west and east groups across the IT, which may have been restricted to allopatric stable climatic areas (Fig. 3). After the disappearance of the barrier (either geographical or ecological), populations are expected to expand to ecologically suitable areas. However, environmental differences in the landscape and biological interactions (e.g., competitive exclusion) between taxa with similar niches may slow down the recolonization process (Anderson et al., 2002), acting as a reinforcement that may help in maintaining differentiation over time.
4.4
Ecological niche differentiation
Estimating ecological niches using intraspecific lineages or subspecies as units can represent their niche more accurately, without ignoring local adaptive responses as would occur when estimating models at the species level (Smith et al., 2019). Our comparisons showed that the ecological niches of the three L. verreauxi lineages are non-equivalent. However, our analyses detected niche similarity in the NA-NCA comparison, but niches were not similar for the clades found to the east of the IT. In the case of NA–NCA, both lineages are occupying different environments and climatic conditions, but ecological niches are similar, and therefore may evolve slowly (Peterson et al., 1999).
The similarity in the ecological niches of the NA and NCA may explain, in part, their sympatric distribution in the southern IT. For the clades eastern of the IT, comparisons in NA-YUC and YUC-NCA indicated that niches are not equivalent nor similar, suggesting niche divergence.
According to range-breaking tests, environmental differences across the north to south ecological gradient in the IT do not represent an abrupt environmental barrier for dispersal of the populations of L. verreauxi between habitats. Thus, our results suggest that the secondary contact zone in southern IT may have less suitable conditions for dispersal of NCA, than for NA. Thus, NCA populations are seemingly less tolerant to dry forest climatic conditions on Pacific slope of the IT, which likely prevents its dispersal further west. NA populations, however, seem to be tolerant to both humid and dry conditions, suggesting that other factors, such as the role of biological interactions at local scales (e.g., competition) between both lineages should be considered in maintaining this limited secondary contact zone. A further possibility is related to the seasonal differences in the precipitation between the dry forests of the southern IT and the mesic forests of the Soconusco region in coastal Chiapas (Binford, 1989): the seasonal increase of humidity would allow the dispersal of the mesic-adapted NCA further north to be in contact with NA in the southern IT. We suggest that this secondary contact zone, has been maintained through limited ecological suitability for NCA to the west of the IT, while to the east, divergence may have been maintained via ecological interactions or seasonal movements. We lack data on lineages north of the IT, where NA and YUC may overlap due to the absence of conspicuous barriers; a shallow genetic divergence however between these two lineages, suggests a corridor across the Usumacinta River drainage.
The Isthmus of Tehuantepec is an area long-recognized as a geographic limit for several North American bird families (Navarro-Sigüenza et al., 2014), as well as a promoter of speciation and ecological differentiation, particularly for montane bird taxa, which have different species divided by the Tehuantepec Valley (Barber and Klicka, 2010). However, its role as a barrier for lowland bird taxa is still poorly understood (see Binford, 1989). Patterns of divergence in lowland populations in the Isthmus region should be viewed with caution due to a dynamic biotic history, which due to changes in ecological conditions has likely promoted ecological niche differentiation, as suggested in several bird taxa (Castillo-Chora et al. 2021; Sánchez-González et al, 2021, 2023; Llanes-Quevedo et al., 2022).
Declaration of competing interest
The authors declare no conflict of interest.
Ethics statement
Specimens and tissues used in this project were collected under permits granted by the Secretaría del Medio Ambiente y Recursos Naturales (SEMARNAT) of the Mexican Government.
CRediT authorship contribution statement
Orlando J. Espinosa-Chávez: Writing – original draft, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. Adolfo G. Navarro-Sigüenza: Writing – review & editing, Supervision, Data curation. Hernando Rodríguez-Correa: Writing – review & editing, Supervision. Luis A. Sánchez-González: Writing – review & editing, Writing – original draft, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.
Acknowledgments
This paper is part of the fulfillments of OJE-C master’s degree project in the Posgrado en Ciencias Biológicas, UNAM, who also had a CGEP Scholarship for UNAM postgraduate studies. Many thanks to the financial support for this study by PAPIIT- UNAM grant (IN222817) to LAS-G, and the François Vuilleumier Fund for Neotropical Bird Research from the Neotropical Ornithological Society (NOS) awarded to OJE-C. Special tanks to Uri O. García-Vázquez and David A. Prieto-Torres for comments to early drafts. We appreciate comments by two anonymous reviewers, which greatly improved our paper.
Arroyo FA, Pawlowska TE, Choat JH, Clements KD, Angert ER. Recombination contributes to population diversification in the polyploid intestinal symbiont Epulopiscium sp. type B. ISME J. 2019;13: 1084-97.
Barnes A, Kaur A, Augenbraun M. An unusual presentation of prostatic abscess due to Actinomyces turicensis and Peptostreptococcus. Cureus. 2020;12: e8665.
Berry D. The emerging view of Firmicutes as key fibre degraders in the human gut. Environ Microbiol. 2016;18: 2081-3.
BirdLife International. Gyps himalayensis (amended version of 2017 assessment). The IUCN red list of threatened species 2017: e. T22695215A118594518. 2017. . Accessed 10 Sept 2020.
Blumstein DT, Rangchi TN, Briggs T, de Andrade FS, Natterson-Horowitz B. A systematic review of carrion eaters' adaptations to avoid sickness. J Wildl Dis. 2017;53: 577-81.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QⅡME allows analysis of high-throughput community sequencing data. Nat Method. 2010;7: 335-6.
Chung O, Jin S, Cho YS, Lim J, Kim H, Jho S, et al. The first whole genome and transcriptome of the cinereous vulture reveals adaptation in the gastric and immune defense systems and possible convergent evolution between the Old and New World vultures. Genome Biol. 2015;16: 215.
Colston TJ, Jackson CR. Microbiome evolution along divergent branches of the vertebrate tree of life: what is known and unknown. Mol Ecol. 2016;25: 3776-800.
David LA, Maurice CF, Carmody RN, Gootenberg DB, Button JE, Wolfe BE, et al. Diet rapidly and reproducibly alters the human gut microbiome. Nature. 2014;505: 559-63.
de Paepe K, Kerckhof FM, Verspreet J, Courtin CM, Van de Wiele T. Inter-individual differences determine the outcome of wheat bran colonization by the human gut microbiome. Environ Microbiol. 2017;19: 3251-67.
Diakou A, Norte AC, de Carvalho IL, Núncio S, Nováková M, Kautman M, et al. Ticks and tick-borne pathogens in wild birds in Greece. Parasitol Res. 2016;115: 2011-6.
Edgar RC. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Method. 2013;10: 996-8.
Gao H, Jiang Q, Ji H, Ning J, Li C, Zheng H. Type 1 diabetes induces cognitive dysfunction in rats associated with alterations of the gut microbiome and metabolomes in serum and hippocampus. BBA-Mol Basis Dis. 2019;1865: 165541.
Garg S, Kaul SC, Wadhwa R. Cucurbitacin B and cancer intervention: chemistry, biology and mechanisms. Int J Oncol. 2018;52: 19-37.
Ghonimy A, Zhang DM, Farouk MH, Wang Q. The impact of carnitine on dietary fiber and gut bacteria metabolism and their mutual interaction in monogastrics. Int J Mol Sci. 2018;19: 1008.
Grigor'eva IN. Gallstone disease, obesity and the Firmicutes/Bacteroidetes ratio as a possible biomarker of gut dysbiosis. J Pers Med. 2021;11: 13.
Hird SM. Evolutionary biology needs wild microbiomes. Front Microbiol. 2017;8: 725.
Houston DC, Cooper JE. The digestive tract of the whiteback griffon vulture and its role in disease transmission among wild ungulates. J Wildl Dis. 1975;11: 306-13.
Huang P, Zhang Y, Xiao K, Jiang F, Wang H, Tang D, et al. The chicken gut metagenome and the modulatory effects of plant-derived benzylisoquinoline alkaloids. Microbiome. 2018;6: 211.
Huang G, Wang X, Hu Y, Wu Q, Nie Y, Dong J, et al. Diet drives convergent evolution of gut microbiomes in bamboo-eating species. Sci China Life Sci. 2021;64: 88-95.
Jarvis ED, Mirarab S, Aberer AJ, Li B, Houde P, Li C, et al. Whole-genome analyses resolve early branches in the tree of life of modern birds. Science. 2014;346: 1320-31.
Kau AL, Ahern PP, Griffin NW, Goodman AL, Gordon JI. Human nutrition, the gut microbiome and the immune system. Nature. 2011;474: 327-36.
Keenan SW, Engel AS, Elsey RM. The alligator gut microbiome and implications for archosaur symbioses. Sci Rep. 2013;3: 2877.
Kelly D, Yang L, Pei Z. Gut microbiota, Fusobacteria, and colorectal cancer. Diseases. 2018;6: 109.
Kocijan I, Prukner-Radovčić E, Beck R, Galov A, Marinculić A, Sušić G. Microflora and internal parasites of the digestive tract of Eurasian griffon vultures (Gyps fulvus) in Croatia. Eur J Wildl Res. 2009;55: 71-4.
Kogut MH. Issues and consequences of using nutrition to modulate the avian immune response. J Appl Poultry Res. 2017;26: 605-12.
Kohl KD. Diversity and function of the avian gut microbiota. J Comp Physiol B. 2012;182: 591-602.
Kohl KD, Connelly JW, Dearing MD, Forbey JS. Microbial detoxification in the gut of a specialist avian herbivore, the Greater Sage-Grouse. FEMS Microbiol Lett. 2016;363: fnw144.
Kumar S, Stecher G, Li M, Knyaz C, Tamura K. MEGA X: molecular evolutionary genetics analysis across computing platforms. Mol Biol Evol. 2018;35: 1547-9.
Lagier JC, Armougom F, Million M, Hugon P, Pagnier I, Robert C, et al. Microbial culturomics: paradigm shift in the human gut microbiome study. Clin Microbiol Infect. 2012;18: 1185-93.
Langille MGI, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes JA, et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol. 2013;31: 814-21.
Lee WY. Avian gut microbiota and behavioral studies. Korean J Ornithol. 2015;22: 1-11.
Lee WJ, Hase K. Gut microbiota-generated metabolites in animal health and disease. Nat Chem Biol. 2014;10: 416-24.
Levin D, Raab N, Pinto Y, Rothschild D, Zanir G, Godneva A, et al. Diversity and functional landscapes in the microbiota of animals in the wild. Science. 2021;372: e5352.
Ley RE, Hamady M, Lozupone C, Turnbaugh PJ, Ramey RR, Bircher JS, et al. Evolution of mammals and their gut microbes. Science. 2008;320: 1647-51.
Lu X, Ke D, Zeng X, Gong G, Ci R. Status, ecology, and conservation of the Himalayan griffon Gyps himalayensis (Aves, Accipitridae) in the Tibetan plateau. Ambio. 2009;38: 166-73.
Margalida A, Oliva-Vidal P. The shadow of diclofenac hangs over European vultures. Nat Ecol Evol. 2017;1: 1050.
Mateos-Hernández L, Crespo E, de la Fuente J, de la Lastra JMP. Identification of key molecules involved in the protection of vultures against pathogens and toxins. In: Radis-Baptista G, editor. An integrated view of the molecular recognition and toxinology: from analytical procedures to biomedical applications. Rijeka: InTech; 2013. p. 241-65.
Meadows JA, Wargo MJ. Carnitine in bacterial physiology and metabolism. Microbiology. 2015;161: 1161-74.
Meng X, Lu S, Lai XH, Wang Y, Wen Y, Jin D, et al. Actinomyces liubingyangii sp. nov. isolated from the vulture Gypaetus barbatus. Int J Syst Evol Microbiol. 2017a;67: 1873-9.
Meng X, Lu S, Wang Y, Lai XH, Wen Y, Jin D, et al. Actinomyces vulturis sp. nov., isolated from Gyps himalayensis. Int J Syst Evol Microbiol. 2017b;67: 1720-6.
Meng X, Lu S, Yang J, Jin D, Wang X, Bai X, et al. Metataxonomics reveal vultures as a reservoir for Clostridium perfringens. Emerg Microbes Infect. 2017c;6: e9.
Meng X, Lai XH, Lu S, Liu S, Chen C, Zhou D, et al. Actinomyces tangfeifanii sp. nov., isolated from the vulture Aegypius monachus. Int J Syst Evol Microbiol. 2018;68: 3701-6.
Miyake S, Ngugi DK, Stingl U. Phylogenetic diversity, distribution, and cophylogeny of giant bacteria (Epulopiscium) with their surgeonfish hosts in the Red Sea. Front Microbiol. 2016;7: 285.
Moleón M, Sánchez-Zapata JA, Margalida A, Carrete M, Owen-Smith N, Donázar JA. Humans and scavengers: the evolution of interactions and ecosystem services. Bioscience. 2014;64: 394-403.
Ogada D, Shaw P, Beyers RL, Buij R, Murn C, Thiollay JM, et al. Another continental vulture crisis: Africa's vultures collapsing toward extinction. Conserv Lett. 2016;9: 89-97.
Padayachee T, Nzuza N, Chen W, Nelson DR, Syed K. Impact of lifestyle on cytochrome P450 monooxygenase repertoire is clearly evident in the bacterial phylum Firmicutes. Sci Rep. 2020;10: 13982.
Pinevich AV, Andronov EE, Pershina EV, Pinevich AA, Dmitrieva HY. Testing culture purity in prokaryotes: criteria and challenges. Antonie Van Leeuwenhoek. 2018;111: 1509-21.
Plaza PI, Blanco G, Madariaga MJ, Boeri E, Teijeiro ML, Bianco G, et al. Scavenger birds exploiting rubbish dumps: pathogens at the gates. Transbound Emerg Dis. 2019;66: 873-81.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41: D590-6.
Ramirez C, Coronado J, Silva A, Romero J. Cetobacterium is a major component of the microbiome of giant Amazonian fish (Arapaima gigas) in Ecuador. Animals. 2018;8: 189.
Roggenbuck M, Bærholm Schnell I, Blom N, Bælum J, Bertelsen MF, Sicheritz-Pontén T, et al. The microbiome of New World vultures. Nat Commun. 2014;5: 5498.
Rubinstein MR, Wang X, Liu W, Hao Y, Cai G, Han YW. Fusobacterium nucleatum promotes colorectal carcinogenesis by modulating E-cadherin/β-catenin signaling via its FadA adhesin. Cell Host Microbe. 2013;14: 195-206.
Simon K, Verwoolde MB, Zhang J, Smidt H, de Vries RG, Kemp B, et al. Long-term effects of early life microbiota disturbance on adaptive immunity in laying hens. Poult Sci. 2016;95: 1543-54.
Smiddy M, Sleator RD, Patterson MF, Hill C, Kelly AL. Role for compatible solutes glycine betaine and L-carnitine in listerial barotolerance. Appl Environ Microbiol. 2004;70: 7555-7.
Trevelline BK, Fontaine SS, Hartup BK, Kohl KD. Conservation biology needs a microbial renaissance: a call for the consideration of host-associated microbiota in wildlife management practices. P Roy Soc B-Biol Sci. 2019;286: 20182448.
Waite DW, Taylor MW. Characterizing the avian gut microbiota: membership, driving influences, and potential function. Front Microbiol. 2014;5: 223.
Waite DW, Taylor MW. Exploring the avian gut microbiota: current trends and future directions. Front Microbiol. 2015;6: 673.
Wei F, Wu Q, Hu Y, Huang G, Nie Y, Yan L. Conservation metagenomics: a new branch of conservation biology. Sci China Life Sci. 2019;62: 168-78.
Wu GD, Chen J, Hoffmann C, Bittinger K, Chen YY, Keilbaugh SA, et al. Linking long-term dietary patterns with gut microbial enterotypes. Science. 2011;334: 105-8.
Yao R, Xu L, Lu G, Zhu L. Evaluation of the function of wild animal gut microbiomes using next-generation sequencing and bioinformatics and its relevance to animal conservation. Evol Bioinform Online. 2019;15: 1176934319848438.
Zepeda Mendoza ML, Roggenbuck M, Manzano-Vargas K, Hansen LH, Brunak S, Gilbert MTP, et al. Protective role of the vulture facial skin and gut microbiomes aid adaptation to scavenging. Acta Vet Scand. 2018;60: 61.
Zhai Q, Yin R, Yu L, Wang G, Tian F, Yu R, et al. Screening of lactic acid bacteria with potential protective effects against cadmium toxicity. Food Control. 2015;54: 23-30.
Zhang G, Li C, Li Q, Li B, Larkin DM, Lee C, et al. Comparative genomics reveals insights into avian genome evolution and adaptation. Science. 2014;346: 1311-20.
Zhou C, Wang G, Yu H, Geng Y, Wu W, Tu H, et al. Genome-wide analysis reveals the genomic features of the turkey vulture (Cathartes aura) as a scavenger. Mol Genet Genom. 2019;294: 679-92.
Table
1.FST values (below the diagonal) and DXY values (above the diagonal) for the different L. verreauxi lineages. Values in the parentheses represent the number of migrants (Nm) per generation.
DownLoad:
Email Alerts
RSS Feeds