Species limits and hybridization in Andean leaf-eared mice (Phyllotis)

Abstract

Leaf-eared mice (genus Phyllotis) are among the most widespread and abundant small mammals in the Andean Altiplano, but species boundaries and distributional limits are often poorly delineated due to sparse survey data from remote mountains and high-elevation deserts. Here we report a combined analysis of mitochondrial DNA variation and whole-genome sequence (WGS) variation in Phyllotis mice to delimit species boundaries, to assess the timescale of diversification of the group, and to examine evidence for interspecific hybridization. Estimates of divergence dates suggest that most diversification of Phyllotis occurred during the past 3 million years. Consistent with the Pleistocene Aridification hypothesis, our results suggest that diversification of Phyllotis largely coincided with climatically induced environmental changes in the mid- to late Pleistocene. Contrary to the Montane Uplift hypothesis, most diversification in the group occurred well after the major phase of uplift of the Central Andean Plateau. Species delimitation analyses revealed surprising patterns of cryptic diversity within several nominal forms, suggesting the presence of much undescribed alpha diversity in the genus. Results of genomic analyses revealed evidence of ongoing hybridization between the sister species Phyllotis limatus and P. vaccarum and suggest that the contemporary zone of range overlap between the two species represents an active hybrid zone.

Keywords: Altiplano; Andes; Phyllotini; Sigmodontinae; geographic range limits; introgression; species delimitation.

Figure 1.

Distribution limits of Phyllotis species and geographic sampling coverage in the Central Andes and adjoining lowlands. A) Ranges of Phyllotis mice in the P. darwini species group, based on patterns of morphological and DNA marker variation (Jayat et al., 2021; Ojeda et al., 2021; Steppan and Ramírez, 2015; Storz et al., 2024). B) Distribution of 169 sampling localities, representing sites of origin for 448 Phyllotis specimens used in the survey of cytb and WGS variation.

Figure 2.

Calibrated maximum clade credibility tree showing Bayesian estimates of phylogenetic relationships and divergence times within the genus Phyllotis. Estimates of the 95% Highest Posterior Distributions interval for the divergence times are shown for main clades. Node support is shown only for those cases in which Bayesian posterior probability values were <1. Specimens in the clade labeled ‘P. vaccarum*’ carry cytb haplotypes that group with haplotypes of P. limatus, even though whole-genome sequence data confirmed their identity as P. vaccarum (Storz et al., 2024).

Figure 3.

Maximum clade credibility depicting the delimitation schemes inferred from GMYC (red bars) and PTP (blue bars). Gaps in the vertical bars denote units delimited by each method, and asterisks denote splits with support values >0.75. Continuous gray bars denote current taxonomic designations for nominal species. Terminal labels depict the haplotype classes of sequences that were retained to construct the non-redundant matrix of cytb haplotypes. Specimens in the clade labeled ‘P. vaccarum*’ carry cytb haplotypes that group with haplotypes of P. limatus, even though whole-genome sequence data confirmed their identity as P. vaccarum (Storz et al., 2024).

Figure 4.

Maximum likelihood tree estimated from coding sequence of complete mitochondrial genomes for a set of 11 nominal Phyllotis species. Numbers adjacent to internal nodes denote ultrafast bootstrap support values for each clade. Within the taxon currently recognized as P. darwini, the species delimitation analysis identified two highly distinct subdivisions (see Fig. 3). Representatives of both internal subdivisions form distinct clades in the mitogenome tree, which we labeled ‘P. darwini south’ and ‘P. darwini’ north.

Figure 5.

Genomic variation among species of Phyllotis based on 137 samples representing 11 nominal species. A) Genomic principal component analysis (PCA) of genome-wide variation (PC1 vs PC2). Two distinct clusters of nominal P. darwini specimens, ‘darwini South’ and ‘darwini North’, are distinguished along the PC1 axis. B) Plot of PC1 vs PC3 separates P. limatus and P. vaccarum along the PC3 axis, and reveals a single specimen, UACH9099 (designated P. limatus based on mtDNA haplotype), which has a PC3 score intermediate between the two species. C) Map of collecting localities and distribution limits of P. limatus and P. vaccarum. UACH9099 comes from a site located in a narrow zone of range overlap between the two species in northern Chile. The map also shows the distribution of mice that are identified as P. vaccarum on the basis of whole-genome sequence data, but which carry mtDNA haplotypes that are more closely related to those of P. limatus (denoted as ‘P. vaccarum*’ in the inset tree diagram). D) Structure plot showing clear distinction between P. limatus and P. vaccarum (n=20 and 51, respectively). The putative hybrid specimen, UACH9099, was assigned almost exactly equal ancestry proportions from the two species.

Figure 6.

Windowed PCA of a P. vaccarum x P. limatus hybrid. PC1 was computed in overlapping 1 Mbp windows along the genome for a subset of 50 P. vaccarum (green), 20 P. limatus (blue), and the putative hybrid, UACH9099 (red). Mean PC1 values for each species are shown as white lines and the mean value between both species’ averages is shown as a grey line. UACH9099 features a mosaic genome, with its local ancestry alternating between P vaccarum, P. limatus, or a point intermediate between the two species. (A) Windowed PCA of chromosomes 1-19. (B) High resolution visualization of PC 1 along chromosome 1.

Figure 7.

Revised distribution limits of species in the Phyllotis darwini species group based on mtDNA and WGS data. Filled circles denote collection localities that helped define geographic range limits.