Genomics and morphometrics reveal the adaptive evolution of pikas

Rui-Xiang Tang, Jiao Wang, Yi-Fei Li, Cheng-Ran Zhou, Guan-Liang Meng, Feng-Jun Li, Yue Lan,Megan Price, Lars Podsiadlowski, Yan Yu, Xu-Ming Wang, Yin-Xun Liu,, Bi-Song Yue, Shan-Lin Liu,Zhen-Xin Fan,7,＊, Shao-Ying Liu,＊

1 Key Laboratory of Bioresources and Ecoenvironment (Ministry of Education), College of Life Sciences, Sichuan University, Chengdu,

Sichuan 610065, China

2 Sichuan Academy of Forestry, Chengdu, Sichuan 610081, China

3 Key Laboratory of Birth Defects and Related Diseases of Women and Children of MOE, Department of Pediatrics, West China Second University Hospital, Sichuan University, Chengdu, Sichuan 610041, China

4 BGI-Shenzhen, Shenzhen, Guangdong 518083, China

5 Zoological Research Museum Alexander Koenig, Bonn D-53113, Germany

6 Department of Entomology, China Agriculture University, Beijing 100193, China

7 Sichuan Key Laboratory of Conservation Biology on Endangered Wildlife, College of Life Sciences, Sichuan University, Chengdu,

Sichuan 610065, China

ABSTRACT Pikas (Lagomorpha: Ochotonidae) are small mouselike lagomorphs. To investigate their adaptation to different ecological environments during their dispersal from the Qinghai-Xizang (Tibet) Plateau(QTP), we collected 226 pikas and measured 20 morphological characteristics and recorded habitat information. We also sequenced the genome of 81 specimens, representing 27 putative pika species.The genome-wide tree based on 4 090 coding genes identified five subgenera, i.e., Alienauroa, Conothoa,Lagotona, Ochotona, and Pika, consistent with morphometric data. Morphologically, Alienauroa and Ochotona had similar traits, including smaller size and earlier divergence time compared to other pikas.Consistently, the habitats of Alienauroa and Ochotona differed from those of the remaining subgenera. Phylogenetic signal analysis detected 83 genes significantly related to morphological characteristics, including several visual and hearingrelated genes. Analysis of shared amino acid substitutions and positively selected genes (PSGs) in Alienauroa and Ochotona identified two genes, i.e.,mitochondrial function-related TSFM (p.Q155E) and low-light visual sensitivity-related PROM1 (p.H419Y).Functional experiments demonstrated that TSFM-155E significantly enhanced mitochondrial function compared to TSFM-155Q in other pikas, and PROM1-419Y decreased the modeling of dynamic intracellular chloride efflux upon calcium uptake.Alienauroa and Ochotona individuals mostly inhabit different environments (e.g., subtropical forests) than other pikas, suggesting that a shift from the larger ancestral type and changes in sensory acuity and energy enhancement may have been required in their new environments. This study increases our understanding of the evolutionary history of pikas.

Keywords: Pikas; Genomics; Morphometrics;Adaptive evolution; Sensory and energy functions

lNTRODUCTlON

Pikas (genusOchotona; Lagomorpha: Ochotonidae) are small mouse-like lagomorphs (Wang et al., 2020). Once widely distributed across Eurasia, Africa, and North America, most extant pikas are found in China, especially in the Hengduan Mountains, with a few species inhabiting North America,Central Asia, and the Russian Far-East (Lissovsky et al.,2019; Liu et al., 2017; Wang et al., 2020). Pikas are larger than mice but smaller than rabbits, and exhibit size differences among species. Extant pikas inhabit alpine areas, tundra,steppes, and other cold environments (Wang et al., 2020;Wilson and Smith, 2015) and play important ecological roles in stabilizing plant communities and promoting plant growth due to grazing (Smith At, 2018; Wang et al., 2020; Wilson and Smith, 2015).

The taxonomy and phylogenetic relationships within the genusOchotonahave been explored using different mitochondrial genes, nuclear genes, and exome data (Ge et al., 2013; Wilson DE, 2005; Koju et al., 2017; Lanier and Olson, 2009; Lissovsky, 2014; Liu et al., 2017; Melo-Ferreira et al., 2015; Niu Yd, 2004; Yu et al., 2020). Liu et al. (2017)proposed five subgenera (Conothoa,Ochotona,Pika,Alienauroa, andLagotona) based on cytochromeb(cytb)sequences, while Wang et al. (2020) groupedLagotonaintoPikaand supported only four subgenera (Conothoa,Ochotona,Pika, andAlienauroa) based on whole-genome coding sequences. The genetic basis of high-altitude adaptations in pikas has also been investigated (Henry and Russello, 2013; Rankin et al., 2017; Wang et al., 2020), but with an emphasis on adaptations to hypoxia, ultraviolet radiation, and cold tolerance. However, pikas have occupied various ecological environments during their dispersal from the Qinghai-Xizang (Tibet) Plateau (QTP) to other parts of Eurasia and North America (Wang et al., 2020). Indeed, the major clades/subgenera of extant pikas are classified into four ecotypes according to their distribution and ecological environment. TheConothoamountain group is mainly distributed in the Himalayas and surrounding habitats dominated by talus; theAlienauroaforest group is primarily found along the eastern edge of the QTP in humid forest areas dominated by moss and leaf litter; thePikanorthern group is distributed at high latitudes with rock and steppe habitats; and theOchotonashrub-steppe group is mainly distributed in the northern areas of the QTP dominated by shrub-steppe (Wang et al., 2020).

As pika species inhabit different environments, and environmental stress was a major driving force during pika evolution (Parsons, 2005; Kristensen et al., 2020), we hypothesized that certain functions and pathways may be under different selection in different pikas. To study the evolutionary history of pikas, we measured 20 external,cranial, and dental characteristics, and recorded habitat information of 226 adult pikas collected from 107 sampling sites in China. We also sequenced the genomes of 81 specimens representing 27 putative species of pika and reconstructed a genome-wide tree based on 4 090 coding genes. Based on extensive genomic, morphological, and habitat datasets, we: (1) assessed the morphological variation and evolution of different pikas; (2) explored the adaptations of different pikas and the relationship between molecular genetics and phenotypic characteristics; and (3) performed functional experiments based on cellular mutations.

MATERlALS AND METHODS

Ethics statement

All samples were obtained following the Guidelines of the American Society of Mammalogists (ASM guidelines) (Sikes et al., 2011) and the laws and regulations of China for the protection of terrestrial wild animals (State Council Decree 1992). The BGI Institute of Review Board of Bioethics and Biosafety (BGI-IRB) (No. FT17005) approved all collection and research protocols. Voucher specimens were deposited in the Sichuan Academy of Forestry, Chengdu, China.

Sample and data collection

In this study, we selected 226 pika specimens collected from 107 sampling sites in China over the past 30 years (33 samples were provided by the Northwest Institute of Plateau Biology, Chinese Academy of Sciences, Xining, China).Specimens were captured by randomly placing cage-type small-animal traps at the sampling sites. After capture, the pikas were anesthetized with ether on the spot and measured for external features. Of the captured individuals, 226 adult specimens were transported to the laboratory for processing,and the rest were releasedinsitu. Upon return to the laboratory, the 226 adult specimens were anesthetized and euthanized by cervical dislocation. The pikas were then dissected, and tissue samples were extracted and stored in 95% analytical alcohol. All specimens were adults with intact skulls for morphometric analyses. Ecological conditions were recorded at the time of field sampling (e.g., elevation, habitat,landscape type). Species identification using morphological characteristics was based on the Taxonomic Index of theOchotonain China (Feng, 1985), Taxonomic list and distribution of mammal species and subspecies in China(Wang, 2003), and previous articles (Liu et al., 2017; Wang,1988; Gong et al., 2000), as well as drawings and morphological descriptions in the literature (Feng, 1985;Hoffmann Rs, 2005; Smith, 2009; Sokolov et al., 2009; Wilson et al., 1993). We also sequenced the whole genomes of 81 individuals, although some samples were sequenced at low coverage.

High-throughput sequencing

Total genomic DNA was extracted from muscle and liver tissue from the 81 samples using a Gentra Puregene Tissue Kit (Qiagen, USA) according to the manufacturer’s protocols.We constructed genomic libraries (paired-end libraries with insert sizes of ～300 bp) following BGISEQ 500 (Huang et al.,2017) or Illumina Hiseq 2500 (USA) (Zhang, 2021) protocols(details in Supplementary Table S1). We then sequenced each library and generated ＞45 Gb of data (10-30X) for one representative specimen of each morphologically distinct species, and ～5 Gb of data (1-3X) were generated for the remaining specimens of each species. In the high-coverage datasets, low-quality reads were excluded if one or more of the following criteria were met: (1) N-content ＞10%; (2)adapter-contaminated reads with adapter sequences overlapping reads by ＞50%; (3) read length below Q10 ＞20%.For the low-coverage datasets, the low-quality read filtering criteria and sequencing platform information are listed in Supplementary Table S1. In total, almost 2 Tb of 150-bp paired-end data were obtained.

Measurement of morphological characteristics

We measured 20 phenotypic characteristics in the 226 adults with intact skulls to construct the morphological dataset.Individuals were identified as adults by inspection of their teeth(see detailed criteria in Supplementary Methods). We measured external, cranial, and dental characteristics for all specimens. External measurements were recorded upon capture in the field to an accuracy of 0.5 mm. Cranial and dental characteristics were measured in the laboratory using an electronic vernier caliper to an accuracy of 0.01 mm.External measurements included head and body length (HBL),ear length (from notch to top of ear, EL), and hind foot length(HFL) excluding claws. Cranial measurements included skull greatest length (SGL), skull basal length (SBL), occipital condyle to nasal bone length (OCNL), zygomatic breadth (ZB),minimal interorbital width (minimum distance across frontal bone between orbits, IOW), skull height (vertical distance from ventral surface of bullae to top of cranium, SH), auditory bulla length (ABL), eyepit (eye socket) length (maximum inner diameter of long axis of eyepit, EPL), eyepit breadth(maximum inner diameter of short axis of eyepit, EPB), nasal bone length (NBL), and nasal bone breadth (NBB). Dental measurements included upper molar occlusal surface length(UOSL) and lower molar occlusal surface length (LOSL). We also recorded four discrete morphological features, i.e.,presence of congenital tragus (CT); presence of oval foramen(OF); number of incisive and palatal foramen (IF and PF); and shape of skull (Skull). In addition, we calculated ratios between features, including ratio of ear length to head and body length (EHR) and ratios of auditory bulla length, eyepit length, and eyepit breadth to skull greatest length (ASR,ESLR, ESBR).

We analyzed morphometric variation among the 226 adult specimens using principal component analysis (PCA) in R v.4.0.2 (Null et al., 2011). We applied the Kaiser-Meyer-Olkin(KMO) and Bartlett’s tests in the R Psych package v2.1.9(Revelle, 2020) to check PCA fitness. We performed Permutational Multivariate Analysis of Variance(PERMANOVA) in the R Vegan package v2.5.7 (Oksanen et al., 2019) to examine statistical differences at the species and subgenus levels. Heatmaps were drawn using the R package ComplexHeatmap v2.10.0 (Gu et al., 2016). Box plots were drawn using the R package ggplot2 v3.3.5, with the significance threshold determined using thet-test.

Molecular data pre-processing and construction of ortholog datasets

We downloaded the reference genome and full gene dataset of the North American Pika (Ochotonaprinceps, vOchPri2.0)from Ensemble (Hubbard et al., 2002) to obtain the orthologous genes of our sequenced samples. We then obtained the lineage-specific single-copy ortholog list of the Euarchontoglires group from the BUSCO database (v10.1)(Sim?o et al., 2015) and extracted the corresponding orthologous genes from the full gene dataset ofO.princeps.To avoid mapping uncertainty that may be caused by paralogs, genes with high similarity within the gene set were removed using BLASTn (v2.6.0+, e value＜1e-5) (Altschul et al., 1990). Subsequently, we obtained aO.princepssinglecopy orthologous gene dataset with 5 684 complete genes.We aligned the whole-genome sequencing (WGS) data of each sequenced individual to the ortholog dataset using BWAMEM v0.7.17 (Li, 2013) with default parameters and obtained the corresponding genes using the consensus calling function in BCFtools v.1.10.2 (Danecek et al., 2021). High-quality coding sequences (CDS) were extracted from the mapping results and subjected to pipeline analysis, as detailed in the Supplementary Methods.

We aligned the corresponding protein sequences using MAFFT v.7.453 (Katoh and Standley, 2013) and conducted CDS alignment using PAL2NAL v14 (Suyama et al., 2006)based on the protein alignments to obtain an individual-level nuclear gene CDS dataset. Representative sample sequences of each morphologically distinct species with the highest sequencing depth were selected to construct a species-level nuclear gene CDS dataset. In addition, four species with published genomes and protein gene sequences (Oryctolagus cuniculus,Peromyscusmaniculatusbairdii,Musmusculus,andRattusnorvegicus) were designated as the outgroups.Their sequence data were downloaded from the NCBI database and combined with the corresponding CDS dataset of pikas. Species-level CDS was aligned using PRANK v.140603 (Loytynoja, 2014) based on protein alignments.

Phylogenetic analyses

For the combined individual-level nuclear gene dataset, we used ModelTest-NG v0.1.6 (Darriba et al., 2020) to identify the best model for each gene based on corrected Akaike Information Criterion (AICc) and the recommended parameters for tree construction. We then inferred the phylogenetic tree using RAxML v8.2.12 (Stamatakis, 2014)with the parameter -m GTRGAMMAI and 100 bootstrap replicates. The best-scoring maximum-likelihood (ML) trees for individual genes and their bootstrap replicate trees were used as input files for ASTRAL-III v.5.15.1 (Zhang et al., 2018) to generate a nuclear gene phylogenetic tree under the multilocus bootstrapping model. We also ran ASTRAL-III without a multi-locus bootstrapping model. The stability of the tree topology was confirmed using default parameters with the dataset containing only the top-scoring ML trees and the dataset containing the processed top-scoring ML trees (with branches showing bootstrap support ＜10% removed), as recommended by ASTRAL. The final species tree was achieved using ASTRAL-III based on the multispecies coalescent model. Due to the limitations of ASTRAL when estimating branch length (i.e., branch lengths are only estimated for internal not terminal branches, are in coalescent units, and are prone to underestimation due to statistical noise in gene tree estimation), direct use of ASTRAL results was not conducive for subsequent analysis. Consequently, branch lengths of the final phylogenetic tree were re-estimated using substitutions per site with ExaML v.3.0.22 (Kozlov et al.,2015). We used the ASTRAL bootstrap replicate trees and determined bootstrap support in the phylogenetic tree using RAxML with the “-f b” parameter.

Divergence time inference

Divergence times of the species tree were estimated based on nuclear genes using MCMCTree in the PAML v.4.9h package(Yang, 2007) with approximate likelihood determined using the“REV” (GTR, model=7) model. All ambiguous codon sites were removed from the “species-level nuclear gene set” and only the third codon sites of them were extracted to generate a new supergene. TrimAl v1.4.15 (Capella-Gutierrez et al.,2009) was used to refine the alignment. The final species tree only retained topological information and was extracted as the input tree for subsequent dating analysis. Fossil calibration points were sourced from previous research (Meredith et al.,2011) and the Paleobiology Database (https://paleobiodb.org/), setting priors of divergence time for the tree nodes. The root was set to 79.5 million years ago (Ma). TheRodentia-Lagomorphasplit,Leporidae-Ochotonidaesplit, andMus-Rattussplit were constrained at 71.5-94.1, 47.4-56.9, and 11-13 Ma, respectively. In addition, the oldest extant pika fossils (0.116-0.160 Ma) were used to calibrate the occurrence ofOchotonidae. The calibration constraints were specified with soft boundaries using 2.5% tail probabilities above and below their limits using the built-in function in MCMCTree. The prior for the mean substitution rate was estimated using BaseML. The independent rate model(clock=2 in MCMCTree) was used to specify the rate priors for internal nodes. The overall substitution rate (rgene_gamma)was set to G (1, 6.3085) and the rate drift parameter(sigma2_gamma) was set to G (1, 4.5). Each run was executed with 5 000 000 iterations as burn-in and sampling every 200 iterations until 100 000 samples were collected.Finally, MCMC analyses were run thrice with different random seed numbers to check convergence, with similar results found.

Phylogenetic comparative analyses

To investigate the relationship between molecular genetics and phenotypic characteristics in extant pikas, we measured phylogenetic signals, i.e., degree to which related species are phenotypically similar, using the Blomberg’s K statistic(Blomberg et al., 2003) implemented in RASP v4 (Yu et al.,2020). We calculated the average values of the continuous traits for each species, then extracted single-gene phylogenetic trees with complete representative species from the previously obtained best-scoring ML trees. We estimated the phylogenetic signal for each trait using different gene phylogenetic trees and the final species tree.

We selected the best-fit model by evaluating the AICc values of different models using the R phytools v.0.7-90 package (Revell, 2012). Among the models, the morphological characteristics with continuous data were more suitable for the Brownian motion model, while the discrete data were more suitable for an equal-rates (ER) model. According to the fitness test results, we performed fast estimation of the ML ancestral states for the continuous traits following the Brownian motion model and computed 95% confidence intervals (CIs) for these estimates using the fastAnc function in R phytools v.0.7-90 package (Revell, 2012). We evaluated the evolutionary traits of the discrete characteristics using the ace function with an ER model in the R ape v.5.4-1 package(Paradis and Schliep, 2019).

Evolutionary analyses

Based on morphological clustering, species inAlienauroaandOchotonaexhibited more similar morphological characteristics than other species, butO.curzoniaeandO.dauuricashowed significant deviations in morphological PCA and inhabited different environments compared to other species withinOchotona(Supplementary Table S2). A separate tree was generated by pruning these two species from the final species tree to avoid outlier effects on evolutionary analysis. Each species node in these two subgenera and their ancestral branches were then labeled as foreground clades.Subsequently, all gene sequences of the remaining 25 species were extracted from the species-level nuclear gene CDS dataset.

Selection analyses were performed using the CODEML program in the PAML package (Yang, 2007). We used the branch-site model to identify positively selected genes (PSGs)in the foreground clades. For each gene, we compared model A (“model 2, NSsites 2, fix_omega 0, omega 1; some sites in the foreground branch are under positive selection) against model A null (“model 2, NSsites 2, fix_omega 1, omega 1”; all sites evolved neutrally) and used theChi-square test to evaluate the significance of the compared likelihood ratio tests(LRTs). Genes showing significant differences (P≤0.05) were considered as PSGs. In accordance with Yu et al. (2016), we searched for common amino acid substitutions in the foreground and background branches of the PSGs. We expanded the sample size to incorporate all individuals to exclude variation in individual species.

Functional gene annotations were downloaded from OrthoDB v10.1 (https://www.orthodb.org/) including Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and InterPro domains. Gene symbols were converted from the reference genome ofO.princeps. Enrichment analyses were performed using g:Profiler (https://biit.cs.ut.ee/gprofiler/gost) (Raudvere et al.,2019), with a Benjamini-Hochberg false discovery rate (FDR≤0.05) significance threshold and the background gene set derived from the Ensembl annotation ofO.princeps(version:e104_eg51_p15_3922dba). Protein-protein interaction networks of PSGs were constructed using the STRING(Szklarczyk et al., 2017) protein interaction database.Cytoscape v3.8.2 (Shannon et al., 2003) was used to visualize the protein-protein interaction network.

lntracellular chloride ion measurement

We tested the functional effects of the substitution PROM1-H419Y following previous research (Hori et al., 2019). Mouse embryonic fibroblasts (MEFs) were maintained in 1640 medium (Gibco, Life Technologies, USA) with 10% fetal bovine serum (Gibco, Life Technologies, USA) supplemented with penicillin-streptomycin (Gibco, Life Technologies, USA).The chloride-sensitive fluorescent indicator MQAE (MCE,USA) was used according to the manufacturer’s instructions.Briefly, the MEFs were incubated with 5 mmol/L MQAE in Krebs-HEPES buffer (20 mmol/L HEPES-NaOH, 128 mmol/L NaCl, 2.5 mmol/L KCl, 2.7 mmol/L CaCl2, 1 mmol/L MgSO4,16 mmol/L glucose) for 1 h at 37 °C. Calcium ionophore was added at a final concentration of 15 μmol/L (A23187, MCE)and temporal changes in chloride ions were detected using a FV3000 confocal microscope (Olympus, Japan) at 1 min intervals.

Cell Mito Stress Test

We performed cell mitochondrial functional analysis to investigate the functional changes in the substitution TSFMQ155E. A Seahorse XF Cell Mito Stress Test Kit (Agilent) was used to evaluate the oxygen consumption rate (OCR). Briefly,10 000 human umbilical vein endothelial cells (HUVECs) were seeded into an 8-well Agilent Seahorse XF Cell Culture Microplate, separately. Then, 1.5 μmol/L oligomycin, 2.0 μmol/L FCCP, and 0.5 μmol/L rotenone-antimycin A were dissolved in an assay medium and loaded into the sensor cartridge. Details on plasmid construction are provided in the Supplementary Methods.

RESULTS

Taxon sampling and datasets

We measured 20 morphological characteristics, including external, cranial, and dental characteristics, of 226 pika specimens collected over the past 30 years (see Methods,Supplementary Table S2). We also obtained environmental information on the collected specimens, including elevation,habitat, landscape type, and microhabitat (Supplementary Table S2). The specimens were assigned to five subgenera and 27 species (nominal or tentatively named) collected from 12 provinces and autonomous regions of China. Many specimens were topotypes (i.e., specimens collected at locality of original type) (Figure 1A; Supplementary Table S1).The complete analysis pipeline is shown in Supplementary Figure S1.

We sequenced the genome of 81 individuals, resulting in 2.09 Tb of clean bases. Four samples were removed due to low library quality and thus were not used for subsequent analyses (Supplementary Table S1). We obtained 5 684 fulllength single-copy orthologous genes spanning 872 931 477 nucleotides (coding sequence) of the 77 specimens fromO.princeps, which were used for nuclear gene dataset construction. After removal of low confidence genes and potential paralogous genes, we retained 4 090 coding genes(“nuclear gene CDS dataset”), with average, maximum, and minimum coverages of 83.94%, 93.67%, and 63.40% for all samples, respectively (Supplementary Table S3).

Molecular phylogeny and morphological data

We reconstructed the phylogeny of extant pikas using nuclear genes. All pika samples were clearly separated from theOryctolaguscuniculus(Leporidae) outgroup and constituted a monophyletic group in the nuclear gene phylogenetic tree.Furthermore, based on the phylogenetic tree (Supplementary Figure S2), the genusOchotonawas recovered as five independent evolutionary branches, representing the five subgeneraAlienauroa,Conothoa,Lagotona,Ochotona, andPika, consistent with the morphological results (Figure 1A, B).The inter-relationship within the genus was (Conothoa,(Alienauroa, ((Lagotona,Pika),Ochotona))), withConothoabeing the most basal clade.

We used the same method as Wang et al. (2020), but with a larger dataset and more loci for estimating the divergence time. We computed divergence times using the species tree based on nuclear genes and a Bayesian relaxed clock(MCMCTree) (Yang, 2007). From our estimation, extant pikas originated 12.05 Ma in the middle Miocene (Figure 1C),consistent with estimates from previous studies (11.62-14.65 Ma) (Ge et al., 2012; Koju et al., 2017; Lanier and Olson,2009; Wang et al., 2020). Extant pikas rapidly diverged in their early speciation from 12.05 Ma to 10.65 Ma.Conothoawas the earliest clade to diverge (12.05 Ma), and the last common ancestor ofConothoaexisted in the late Miocene (7.54 Ma).Alienauroadiverged 11.23 Ma in the Miocene, and the last common ancestor appeared in the late Neogene (4.29 Ma).Finally,Ochotonadiverged 10.65 Ma, with the last common ancestor at 3.39 Ma, andLagotonaandPikadiverged 8.43 Ma. Compared with previous studies, the larger dataset and loci used for estimation here may result in later divergence times at both the subgenus and species level.

After obtaining the phylogenetic relationship, we next analyzed the morphological data. Initial observations of external, cranial, and dental characteristics revealed 27 distinct patterns (Supplementary Table S2), each representing a putative species of pika. This included 26 described species(Liu et al., 2017; Wang et al., 2020) and one previously evaluated but unnamed new species (Wang et al., 2020). We selected one species from each subgenus and noted its skull and external ear characteristics to distinguish unique features(Figure 1B). For example,O.huanglongensis(Alienauroa) had a flatter skull shape than the other subgenera and a distinctive triangular protrusion (congenital tragus) on its ear not found in any other subgenus. The incisive and palatal foramina were separate inO.mantchurica(Pika) but connected in the other subgenera. Additionally,O.gloveri(Conothoa) had an oval foramen, which had disappeared in the other subgenera(Figure 1B).

Figure 1 Geographical distribution, morphological characteristics, and species tree with divergence times of pikas

Taxonomic identifications were corroborated by PCA(Figure 2A), which indicated that the five subgenera were distinct from each other, although some samples were highly similar in several morphological characteristics. At the subgenus level,Conothoa(blue) andPika(orange) were most similar andAlienauroa(green) andOchotona(purple) were most similar, with their centers closer than to any other subgenera, despite not being evolutionarily close (Figure 1C).The PERMANOVA test results showed that the morphological characteristics of the five subgenera differed significantly(P＜0.001), with significant differences between each subgenus pair (P＜0.014) (Supplementary Table S4).Approximately two-thirds of the studied species differed significantly in the 20 morphological characteristics, indicating that a large proportion of species can be identified using morphometrics alone (see details in Supplementary Results).

Consistency in habitat changes and morphological variations

We also recorded collection site information for the 226 individuals, covering 107 sites ranging from 537-4 690 m in elevation (Supplementary Table S2). Through the combination of our results and other records (Wang et al., 2020), we identified thatAlienauroaandOchotonamainly inhabited lowlatitude subtropical forest, with a five-layer vertical structure(megaphanerophyte, dunga-runga, high shrub, low shrub, and herbaceous). Beyond that,O.thomasiandO.nubrica(bothOchotona) inhabited high-altitude dense scrub with an abundant herbaceous layer, whereasO.dauuricaandO.curzoniae(bothOchotona) inhabited prairie or shrub habitat with less vegetation. The other three subgenera (Pika,Lagotona, andConothoa) resided in environments with simple vegetation structures, such as scree, arid valley, prairie, and desert scrub. Several species (e.g.,O.roylii,O.gaoligongensis,O.nigritia,O.mantchurica, andO.coreana)lived along the edge of forest, with a simple two-layer vertical structure (megaphanerophyte and shrub (or herbaceous)).Overall, the habitats ofAlienauroaandOchotonawere markedly different from those of the subgeneraPika,Lagotona, andConothoa(Supplementary Table S2).

To test whether the different ecological environments had different effects on the evolution of morphological characteristics in pikas, we plotted a heatmap using the 20 morphological characteristics. Overall, the morphological characteristics ofAlienauroaandOchotonawere smaller than those of the other three subgenera (Conothoa,Lagotona, andPika), whileO.dauuricaandO.curzoniaewere larger in size than other species in the subgenusOchotona(Figure 2B). Box plots showed that most differences in size were significant between the (AlienauroaandOchotona) and (Conothoa,Lagotona, andPika) groups (Supplementary Figure S3),including EPL, EPB, EL, and ABL (Figure 2C-F). In total, 18 of the 20 morphological characteristics showed significant differences between the above two groups, except for EHR and ASR (Supplementary Figure S3).

We subsequently reconstructed the ancestral states of 24 morphological characteristics along phylogenetic lineages of extant pikas to further investigate their phenotypic evolution(Supplementary Figure S4 and Table S5). Notably, onlyAlienauroaspecies evolved a congenital tragus (CT).Furthermore, consistent with heatmap analysis, after the early divergence ofAlienauroaandOchotona, both evolved similar and smaller morphological characteristics compared with the other pikas and their ancestors (Supplementary and Table S2;fossil evidence in Supplementary Results). For example,AlienauroaandOchotonaindividuals had shorter HBL(Figure 2G), EPL (Figure 2H), and ABL (Figure 2I). In contrast,Conothoa,Lagotona, andPikatraits remained as large as the ancestral types.

Phylogenetic signals and positive selection

To determine the effect of size trends in different phenotypic traits on genetics, we first measured phylogenetic signals. In total, 83 genes significantly related to at least one of the 16 continuous morphological characteristics were identified(Supplementary Table S6; Figure 3A). Among them, EPL had the greatest number of significantly related genes (40),including several genes related to visual ability, such asRPE65,NEIL3, andCROCC. In addition, EL had the second greatest number of significantly related genes (29), including several hearing-related genes, such asJAG1, andLARS2.We also detected other hearing-related genes that were significantly related to non-EL characteristics, such asLRP10in ABL.

Based on phylogenetic signal, morphological variation, and habitat data, we hypothesized that different environments enabled certain genes, especially sensory-related genes, to undergo the same amino acid substitutions, resulting in similar morphological and functional changes inAlienauroaandOchotonaand divergence from other pika subgenera.Therefore, we identified shared amino acid substitutions between (AlienauroaandOchotona) and (Conothoa,Lagotona, andPika). To limit species variation, we included all 73 individuals, after excluding fourO.curzoniaeandO.dauuricaindividuals because they showed significant deviation in morphological PCA (Figure 2A), smaller morphological characteristics (Figure 2B), and simpler environments (Supplementary Table S2) compared to otherAlienauroaandOchotonaspecies. Overall, we identified 408 common amino acid substitutions in 345 genes in allAlienauroaandOchotonapikas.

We identified a total of 181 PSGs in theAlienauroaandOchotonalineages, excludingO.curzoniaeandO.dauurica(Supplementary Table S7). These PSGs were enriched in functional categories related to visual ability, such as photoreceptor cell cilium (GO:0097733), non-motile cilium(GO:0097730), and photoreceptor outer segment(GO:0001750). In addition, abundant human phenotype ontologies (HP) related to visual ability, such as abnormal visual electrophysiology (HP:0030453), abnormal eye physiology (HP:0012373), and retinal dystrophy(HP:0000556), were significantly enriched in the PSGs(Figure 3B; Supplementary Table S8). We performed proteinprotein interaction analysis using all 181 PSGs and only illustrated networks with more than two genes (Figure 3C),which showed that the PSGs were highly correlated to each other. Among the above PSGs,RPE65andNEIL3were also detected in the phylogenetic signal. Other PSGs, such asPROM1,IQCB1,GUCA1B, andRP1, were also related to lowlight visual sensitivity. In addition to vision-related genes,mitochondrial function-related genes were also identified as PSGs, includingTSFMandTUFM.

Functional validation of substitutions in TSFM and PROM1

Of the common substitutions identified, 65 occurred in 53 genes also identified as PSGs (Supplementary Table S7).Combining amino acid mutations, positively selected genes,and the above-mentioned genes related to vision and mitochondrial function, we screened two loci, i.e., low-light visual sensitivity genePROM1(p.H419Y) and mitochondrial function-related geneTSFM(p.Q155E), for subsequent functional verification (Figure 3D). All sampledAlienauroaandOchotonaspecimens had the Prom1-419Y and TSFM-155E substitutions, whereas all other pikas had Prom1-419H and TSFM-155Q (Figure 3D).TSFMis mainly located in the mitochondria and is associated with mitochondrial products and function (Scala et al., 2019). The effects of TSFM-155Q and TSFM-155E on mitochondrial function in HUVECs showed that TSFM-155E significantly improved basal respiration but did not change proton leak (Figure 4A, B;Supplementary Tables S9, S10). Additionally, the maximal respiratory rate was markedly higher in HUVECs that transfected TSFM-155E plasmids, indicating enhancement in maximal electron transport chain activity.

Figure 2 Morphological characteristics of pikas

Figure 3 Phylogenetic signal analysis and shared amino acid substitutions

PROM1drives chloride ion efflux upon intracellular calcium ion uptake (Hori et al., 2019). Here, we investigated whether PROM1-H419Y affectsPROM1function. First, we overexpressed Prom1-419Y or Prom1-419H plasmids in MEFs, then used the chloride-sensitive fluorescent indicator MQAE to detect temporal changes in intracellular chloride ion levels upon calcium uptake. Significant chloride efflux was observed in the MEFs that exogenously expressed Prom1-419H when intracellular calcium uptake was provoked by the calcium ionophore A23187. Taken together, Prom1-419H promoted the modeling of dynamic intracellular chloride efflux upon calcium uptake compared to Prom1-419Y (Figure 4C, D;Supplementary Table S11).

DlSCUSSlON

Figure 4 Functional validation of identified variants on TSFM and PROM1

We collected 226 pika specimens, covering most of their distribution in China (Lissovsky et al., 2019; Liu et al., 2017;Wang et al., 2020), and generated a comprehensive dataset to study their environmental adaptation. Results showed that morphological characteristics can be used to identify pika species, but species with similar morphology need to be further identified using genomic sequences. Furthermore,genomic and morphological analyses supported the classification of five subgenera (Alienauroa,Conothoa,Lagotona,Ochotona, andPika) (Liu et al., 2017). TheAlienauroaandOchotonaspecies were morphologically smaller than the other pika species and their ancestors for nearly all 20 morphological characteristics examined (Figure 2B-I; Supplementary Figures S3, S4). This is consistent with the fossil records for extinct pikas, such asOchotonoides complicidensandOchotonalagreli, which indicate that they were morphologically larger than extant pikas (Institute of Vertebrate Paleontology, 1960) (see Supplementary Results).The smaller body size of extantAlienauroaandOchotonaspecies may be associated with the dispersal of pikas into different habitats with different environmental conditions.Recent study suggested that extant pikas may have originated from high-elevation areas of the QTP, with later dispersal to other parts of Eurasia and North America (Wang et al., 2020).During dispersal, different species of pika would have occupied different ecological environments, which may have influenced certain morphological characteristics. ForAlienauroaspecies, which inhabit mossy forest areas (Wang et al., 2020) (Supplementary Table S2), we identified microhabitats such as talus, moss, shrubs, grass, stone gaps,and bamboo forest. ForOchotonaspecies, which are typical forest and scrub-grassland-steppe dwellers (Wang et al.,2020) (Supplementary Table S2), we identified several microhabitats, including fallen logs, bushes, bamboo forest,talus, and scrub-grassland (Supplementary Table S2).Pika,Lagotona, andConothoaspecies, which belong to the northern and mountain groups (Wang et al., 2020), mainly inhabit scree, arid valley, shrub, prairie, and desert habitats,with much less vegetation than forests. Consistently,according to our results,AlienauroaandOchotonaspecies mostly inhabited subtropical forest habitats, whereas the subgeneraPika,Lagotona, andConothoamainly inhabited scree, arid valley, prairie, and desert scrub habitats(Supplementary Table S2). Dispersal into different environments may require different sensory and energy functions to avoid predators and find food, thus implying natural selection pressure for sensory and energy functions inAlienauroaandOchotonaspecies to adapt to their different environments.

Indeed, the phylogenetic signals suggested that the morphological evolution of different pikas focused on sensory functions, such as vision and hearing (Figure 3A). The identified PSGs inAlienauroaandOchotonaspecies(Supplementary Table S7) were mainly enriched in GO categories and HP ontologies related to visual ability(Supplementary Table S8). Several PSGs were closely related to low-light visual sensitivity, such asPROM1,RPE65,NEIL3,IQCB1,GUCA1B, andRP1. ThePROM1(prominin 1) gene,also known asCD133andAC133, encodes a pentaspan transmembrane glycoprotein, which is concentrated at the base of the photoreceptor outer segment and acts as a key regulator of disk morphogenesis during early retinal development (Fujinami et al., 2020; Maw et al., 2000).Variants ofPROM1can cause diseases such as retinal degeneration, cone-rod dystrophies, and retinitis pigmentosa(Carss et al., 2017; Cehajic-Kapetanovic et al., 2019; Song et al., 2011).RPE65is a protein-coding gene, and its protein is a component of the vitamin A visual cycle of the retina, which provides the 11-cis retinal chromophore of the photoreceptor opsin visual pigments (Jin et al., 2005). Mutations inRPE65can cause various diseases, including Leber congenital amaurosis and early-onset retinal dystrophy, which can lead to poor vision in dim light (Jo et al., 2019; Kumaran et al., 2020;Motta et al., 2020). TheNEIL3(Nei-like DNA glycosylase 3)gene, which belongs to a class of DNA glycosylases homologous to the bacterial Fpg/Nei family, is necessary for normal retinal lamination and retinal neuronal differentiation by rescuing the retinal homeobox knockdown phenotype (Pan et al., 2018). Other identified PSGs, such asIQCB1,GUCA1B,andRP1, are also related to low-light visual sensitivity (Otto et al., 2005; Sato et al., 2005; Silva et al., 2020).

Among the identified PSGs, onlyPROM1had a shared amino acid substitution (p.H419Y). All samples fromAlienauroaandOchotonacontained Prom1-419Y, whereas the other pikas contained Prom1-419H. Thus, we overexpressed the Prom1-419Y and Prom1-419H plasmids in MEFs. Results showed that PROM1-419Y inAlienauroaandOchotonainhibited the modeling of dynamic intracellular chloride efflux upon calcium uptake. Similarly, Prom1-knockout has a weakening effect on mice, resulting in lightdependent photoreceptor cell degeneration, which can be almost completely inhibited when Prom1-knockout mice are maintained in the dark (Dellett et al., 2015; Hori et al., 2019).As observed in the field, these pikas are diurnal animals(Smith At, 2018; Smith and Smith, 1986). This suggests that individuals with Prom1-knockout or reduced function, such as PROM1-419Y inAlienauroaandOchotonaspecies, may be better suited to living in subtropical forests, where additional layers of vegetation result in lower-light environments than found in the scree, arid valley, prairie, and desert scrub habitats ofPika,Lagotona, andConothoaspecies. We must note that vision is determined by many genes and loci, and even thePROM1gene has several mutations associated with retinal functions (Jespersgaard et al., 2019; Maw et al., 2000;Michaelides et al., 2010; Zhang et al., 2007).

In addition to vision-related genes, several mitochondrial function-related genes were also identified as PSGs, such asTSFMandTUFM. Mutations in theTSFMgene, which encodes a mitochondrial translation elongation factor enzyme,are associated with oxidative phosphorylation enzyme deficiency (Smeitink et al., 2006) and sensorineural hearing loss (Scala et al., 2019).TUFMencodes the protein Tu translation elongation factor involved in mitochondrial protein translation (Di Nottia et al., 2017). In addition,TSFMcatalyzes the exchange of guanine nucleotides on theTUFMprotein during elongation in mitochondrial protein translation (Wieden et al., 2002). Here,TSFMcontained one amino acid substitution (p.Q155E), i.e., TSFM-155E in theAlienauroaandOchotonasamples and TSFM-155Q in all remaining pika species. Functional analysis indicated that TSFM-155E inAlienauroaandOchotonaindividuals significantly enhanced mitochondrial function compared to TSFM-155Q in the other pikas (Figure 4A, B). Enhanced mitochondrial function may be related to the smaller body size ofAlienauroaandOchotonaindividuals. The mass-specific metabolic rate declines with increasing body size (Sukhotin et al., 2020). Larger organisms have lower metabolic turnover and energy demand per unit mass than smaller-bodied organisms across different animal taxa, different-sized organisms within a single species(Konarzewski and Ksiazek, 2013), and even within a single individual during ontogenesis (Maino and Kearney, 2014;Moses et al., 2008). We found thatAlienauroaandOchotonaindividuals were morphologically smaller than other pikas and their ancestors (Supplementary Figure S4 and Table S2). In addition, enhanced mitochondrial function can produce more adenosine triphosphate (ATP) to increase energy supply for activity (Gonzalez-Freire et al., 2015; Grevendonk et al., 2021;Hargreaves and Spriet, 2018), which may assist in the adaptation ofAlienauroaandOchotonaspecies to subtropical forest with a five-layer vertical structure through increased athletic ability (Supplementary Table S2). However, this mutation may also contribute to hypoxia and cold tolerance, as previous studies have found that different pikas exhibit adaptive evolution to hypoxia and cold (Rankin et al., 2017;Wang et al., 2020). Thus, further functional tests should be performed to confirm the effects of this mutation.

Our study utilized years of field collections to perform a comprehensive morphological and phenotypic study of 27 extant pika species in China. Analysis supported five subgenera within theOchotonagenus, i.e.,Alienauroa,Conothoa,Ochotona,Lagotona, andPika. Morphometrics were sufficiently robust to identify most pika species and should be used in conjunction with genomics to expand and strengthen pika taxonomy in the future. Our morphological and phenotypic data identified two distinctive subgenera, i.e.,AlienauroaandOchotona, which were smaller in size and retained more PSGs related to sensory acuity and energy enhancement compared to the other studied and ancestral pikas. TheAlienauroaandOchotonaindividuals also inhabited different environments (e.g., subtropical forest), suggesting that the shift from a larger ancestral type and changes in sensory acuity and energy may have been caused by dispersal into new environments or led to dispersal into new environments. Thus, this study provides a comprehensive analysis of pika taxonomy and increases our understanding of the evolution of environmental adaptation in pikas during dispersal. However, potentially important environmental factors should be investigated in future studies to better explore how the ecological environment has impacted the morphological evolution and adaptations of different pikas.

DATA AVAlLABlLlTY

The data that support our findings were published under NCBI BioProjectID PRJNA716776 and CNGB Nucleotide Sequence Archive (CNSA) project CNP0001439. The data were also deposited in GSA under accession No. PRJCA011005 and in Science Data Bank under DOI:10.57760/sciencedb.j00139.00019.

SUPPLEMENTARY DATA

Supplementary data to this article can be found online.

COMPETlNG lNTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRlBUTlONS

Conceptualization: Z.X.F. and S.Y.L. Sample collection:X.M.W., M.K.T., and S.Y.L. Bioinformatics analyses: R.X.T.,J.W., C.R.Z., F.J.L., L.Y., and Y.Y. Functional experiments:Y.F.L. Writing—original draft: R.X.T., J.W., C.R.Z., and Z.X.F.Writing —review & editing: M.P., L.P., B.S.Y., S.L.L., and S.Y.L. All authors read and approved the final version of the manuscript.

ACKNOWLEDGMENTS

We thank Rui Liao for assistance in collecting specimens in the field. We thank the Northwest Institute of Plateau Biology,Chinese Academy of Sciences, Xining, China for supporting our specimen inspection.