The red panda (Ailurus fulgens), an endangered Himalaya-endemic mammal, has been classified as two subspecies or even two species – the Himalayan red panda (A. fulgens) and the Chinese red panda (Ailurus styani) – based on differences in morphology and biogeography. However, this classification has remained controversial largely due to lack of genetic evidence, directly impairing scientific conservation management. Data from 65 whole genomes, 49 Y-chromosomes, and 49 mitochondrial genomes provide the first comprehensive genetic evidence for species divergence in red pandas, demonstrating substantial inter-species genetic divergence for all three markers and correcting species-distribution boundaries. Combined with morphological evidence, these data thus clearly define two phylogenetic species in red pandas. We also demonstrate different demographic trajectories in the two species: A. styani has experienced two population bottlenecks and one large population expansion over time, whereas A. fulgens has experienced three bottlenecks and one very small expansion, resulting in very low genetic diversity, high linkage disequilibrium, and high genetic load.

( A and C ) The Chinese red panda. ( B and D ) The Himalayan red panda. (A and B) The face coat color of the Chinese red panda is redder with less white on it than that of the Himalayan red panda. (C and D) The tail rings of the Chinese red panda are more distinct than those of the Himalayan red panda, with the dark rings being more dark red and the pale rings being more whitish. Photo credit: (A) Yunfang Xiu, Straits (Fuzhou) Giant Panda Research and Exchange Center, China; does not require permission. (B) Arjun Thapa, Institute of Zoology, Chinese Academy of Sciences. (C) Yibo Hu, Institute of Zoology, Chinese Academy of Sciences. (D) Chiranjibi Prasad Pokheral, Central Zoo, Jawalkhel, Lalitpur, Nepal; does not require permission.

The red panda (Ailurus fulgens), an endangered Himalaya-endemic mammal, was once widely distributed across Eurasia but is now restricted at the southeastern and southern edges of the Qinghai-Tibetan Plateau within an altitude range of 2200 to 4800 m ( 3 ). On the basis of differences in morphology (e.g., skull morphology, coat color, and tail ring) and geographic distribution ( Fig. 1 and table S1), red pandas are classified into two subspecies, the Himalayan subspecies (A. f. fulgens Cuvier, 1825) and the Chinese subspecies (A. f. styani Thomas, 1902) ( 4 , 5 ). Morphologically, the Chinese subspecies has much larger zygomatic breadth, the greatest skull length, stronger frontal convexity, more distinct tail rings, and redder face coat color with less white on it ( Fig. 1 ) ( 5 , 6 ). On the basis of these morphological differences, C. Groves even proposed that the two subspecies should be updated as two distinct species: the Himalayan red panda (A. fulgens) and the Chinese red panda (A. styani) ( 6 ). The Nujiang River is considered the geographic boundary between the two species ( 7 ). The Himalayan red panda is distributed in Nepal, Bhutan, northern India, northern Myanmar, and Tibet and western Yunnan Province of China, while the Chinese red panda inhabits Yunnan and Sichuan provinces of China. The subspecies or species classification has remained controversial largely due to the lack of genetic evidence, and their distribution boundary may also be inaccurate because of the morphological similarity of red pandas on both sides of the Nujiang River ( 6 , 8 , 9 ). For instance, the skull size and morphology of red pandas from southeastern Tibet were more similar to those of the Chinese red panda than the Himalayan red panda ( 6 ). Although previous studies attempted to use mitochondrial DNA or microsatellite markers to explore this problem, the very small sample size from the Himalayan red panda and the limited ability of the molecular markers resulted in failure to resolve the species delimitation ( 10 – 12 ). Next-generation sequencing technology not only provides whole-genome data but also enables the identification of Y chromosome sequences in nonmodel animals, which were difficult to obtain previously ( 13 , 14 ). Thus, it is now feasible to use whole genomes, Y chromosomes, and mitochondrial genomes to comprehensively delimit species, subspecies, and populations. Here, with sufficient sampling of the Himalayan red panda, we performed whole-genome resequencing, Y chromosome single-nucleotide polymorphism (SNP) genotyping, and mitochondrial genome assembly of wild red pandas covering most of the distribution ranges of the two species, aiming to clarify species differentiation, population divergence, demographic history, and the impacts of population bottlenecks on genetic evolutionary potential.

The delimitation of species, subspecies, and population is fundamental for insights into the biology and evolution of species and effective conservation management. Traditionally, species, subspecies, or population delimitation is based on reproductive isolation, geographic isolation, and/or morphological differences and does not consider the role of gene flow. The misclassification of basal taxa will result in erroneous or misleading conclusions about the species’ evolutionary history and adaptive mechanisms, and potentially inappropriate conservation management decisions for threatened species ( 1 , 2 ).

RESULTS AND DISCUSSION

Genomic evidence of two phylogenetic species in red pandas We performed whole-genome resequencing for 65 wild red pandas, with an average of 98.7% genome coverage and 13.9-fold sequencing depth for each individual based on the red panda reference genome (belonging to the Chinese red panda) of 2.34 Gb (15). Using the SNP-calling strategy of the Genome Analysis Toolkit (GATK), we identified a total of 4,932,036 SNPs for further analysis (table S4). On the basis of the whole-genome SNPs, the phylogenetic tree, principal components analysis (PCA), and ADMIXTURE results revealed substantial genetic divergence between the two species, providing the first genomic evidence of species differentiation (Fig. 2, B to D). The middle Himalaya population (MH) belonging to the Himalayan red panda was first divergent from the populations of the Chinese red panda (Fig. 2, B and D). Furthermore, four distinct genetic populations were identified: MH (n = 18), eastern Himalaya-Gaoligong (EH-GLG, n = 3 and 13, respectively), Xiaoxiangling-Liangshan (XXL-LS, n = 12 and 8, respectively), and Qionglai (QL, n = 10) (Fig. 2, B to D; fig. S1; and table S5). The individual SLL1 is the only sampled red panda from the Saluli Mountains (SLL), and its genetic assignment implied gene flow between the SLL population and its adjacent XXL and GLG populations (Fig. 2C). Because of the very small sample size, SLL1 was excluded in any population-level analyses. Traditionally, MH, EH, and the GLG individuals at the western side of the Nujiang River were classified as the Himalayan red panda, while the GLG individuals at the eastern side of Nujiang River, XXL, LS, and QL belonged to the Chinese red panda (7). Our results did not support the Nujiang River as the species distribution boundary because the EH and part of the GLG population at the western side of the Nujiang River clustered into a genetic population with other GLG individuals at the eastern side (Fig. 2, B to D). This EH-GLG genetic clustering was supported by morphological evidence that the morphology of red panda skulls from southeastern Tibet (namely, the EH population in this study) was more similar to that of the Chinese red panda than the Himalayan red panda (6). In addition, two individuals from Myanmar (GLG5 and GLG6) also clustered within the EH-GLG genetic cluster, suggesting that the Myanmar population belongs to the Chinese red panda. Thus, we infer that the Yalu Zangbu River, the largest geographic barrier to dispersal between the two species, may be the potential boundary for species distribution (Fig. 2A), although additional samples need to be collected from Bhutan and India to verify this inference. Fig. 2 Population genetic structure based on whole genomes, Y chromosome SNPs, and mitochondrial genomes of red pandas. (A) The geographic distribution of wild red panda samples under the background of habitat suitability. Red, QL population; purple, XXL-LS population; blue, SLL population; pink, EH-GLG; dark red, MH. (B) Maximum likelihood phylogenetic tree based on whole-genome SNPs, with the ferret as the outgroup. The values on the tree nodes indicate the bootstrap support of ≥50%. (C) ADMIXTURE result based on whole-genome SNPs with K = 2 to 7. (D) PCA result based on whole-genome SNPs. (E) Network map based on eight Y chromosome SNP haplotypes. (F) Network map based on 41 mitochondrial genome haplotypes. Within the Chinese red panda, we further found population genetic differentiation. EH-GLG first diverged with XXL-LS-QL and then QL separated from XXL-LS (Fig. 2, B and C). Notably, we did not detect genetic substructure within EH-GLG spanning the famous Three Parallel Rivers (Nujiang River, Lancangjiang, and Jinshajiang), suggesting that the three large rivers did not hinder the gene flow of red pandas. This result is consistent with data from microsatellite markers (12). Our Y chromosome SNP and mitochondrial genome results also supported the substantial divergence between the two species (Fig. 2, E and F; figs. S2 and S3; and tables S6 to S8). The haplotype networks and phylogenetic trees of both eight Y chromosome SNP (Y-SNP) haplotypes from 49 male individuals and 41 mitochondrial genome haplotypes from 49 individuals showed that the MH haplotypes (Himalayan red panda) clustered together and separated from the haplotypes of the Chinese red panda, highlighting the notable genetic divergence between the two species. In summary, regardless of the whole-genome SNPs, Y-SNPs, or mitochondrial genomes, notable genetic differentiation was found between the two species. Our comprehensive investigations reveal two evolutionarily significant units in red pandas. Under the phylogenetic species concept (16), it is reasonable to propose two species: the Himalayan red panda (A. fulgens) and the Chinese red panda (A. styani). This phylogenetic species classification was supported by their morphological differences (6). The Y chromosome SNP and mitochondrial genome results revealed a female-biased gene flow pattern in red pandas (Fig. 2, E and F). Within the Chinese red panda, we observed different phylogeographic patterns between the mitochondrial genome and Y chromosome. The distribution of mitochondrial haplotypes was mixed and was not associated with the geographic sources of the individuals. By contrast, the distribution of Y-SNP haplotypes demonstrated an obvious phylogeographic structure: The haplotypes of EH-GLG were separated from those of XXL-LS-QL, and no shared Y-SNP haplotypes were found. These contrasting phylogeographic patterns reflected a female-mediated historical gene flow, implying female-biased dispersal and male-biased philopatry in red pandas. This dispersal pattern differs from the male-biased dispersal found in most mammals (17) but is similar to that of another bamboo-eating mammal, the giant panda (18, 19).

Species/population demographic and divergence history The pairwise sequentially Markovian coalescent (PSMC) analysis results showed that the demographic history of red panda could be traced back to approximately 3 million years (Ma) ago, and the two red panda species experienced obviously different demographic histories (Fig. 3A). The Chinese red panda from EH-GLG, XXL-LS, and QL experienced similar demographic trajectories: two population bottlenecks and one large population expansion. This species suffered from an obvious population decline approximately 0.8 Ma ago, which coincided with the occurrence of the Naynayxungla Glaciation (0.78 to 0.5 Ma ago). The population decline resulted in the first bottleneck approximately 0.3 Ma ago, mostly likely caused by the Penultimate Glaciation (0.3 to 0.13 Ma ago) (20). After the glaciations, the populations started to expand and reached a climax approximately 50 thousand years (ka) ago. Then, the arrival of the last glaciations again resulted in rapid population decline, and the second bottleneck occurred during the Last Glacial Maximum (~20 ka ago) (20). Fig. 3 Demographic history, divergence, and admixture of two red panda species and their populations. (A) PSMC analysis revealed different demographic histories of the two species, with a generation time (g) of 6 years and a mutation rate (μ) of 7.9 × 10−9 per site per generation. The time axis is logarithmic transformed. (B) Fastsimcoal2 simulation reconstructed the divergence, admixture, and demographic history of red panda species and populations. The time axis is logarithmic transformed, and the number of migrants per year between two adjacent populations is shown beside each arrow. (C) TreeMix analysis detected significant gene flow from the EH-GLG to XXL-LS populations. s.e., standard error. The Himalayan red panda from MH underwent a demographic history differing from that of the Chinese red panda: three population bottlenecks and one small expansion (Fig. 3A). The difference began with the first population bottleneck approximately 0.25 Ma ago. In contrast to the subsequent population recovery of the Chinese red panda, the Himalayan red panda continued to decrease and then went through a second bottleneck approximately 90 ka ago. Afterward, the population started to increase very slowly, but soon the population again decreased due to the last glaciations. The PSMC results showed that even at the climax of population growth (~50 ka ago), the effective population size of the Himalayan red panda was only approximately 35% that of the Chinese red panda. In addition, the Bayesian skyline plot (BSP) analyses based on mitochondrial genomes indicated that both species experienced recent population declines most likely caused by the Last Glacial Maximum, supporting the PSMC results (fig. S4). The different demographic trajectories may result from geographic and climate differences. The Chinese red panda was mainly distributed in the Hengduan Mountains rather than the platform or adjacent edges of the Qinghai-Tibetan Plateau and thus might have suffered less impact of the Pleistocene glaciations. The interglacial warm climate and the vast region of the Hengduan Mountains might have facilitated the rapid population expansion of the Chinese red panda (3). By contrast, the Himalayan red panda lived in the adjacent southern edge of the Qinghai-Tibetan Plateau and might have suffered severe impact of the Pleistocene glaciations. Even during the interglacial period, the geographic proximity to glaciers and limited potential habitat might have restricted this species’ population recovery (21). In Holocene, the climate might have less impact on red panda populations (21), while increasing human activities became the main factor driving recent red panda population declines, which have been detected by microsatellite marker-based Bayesian population simulations (12). We further uncovered the species/population divergence history using Fastsimcoal2 simulation. On the basis of the comparison of alternative population divergence models, we determined the best-support divergence/demography model (Fig. 3B, fig. S5, and table S9). The divergence between the Himalayan (MH) and Chinese red pandas (EH-GLG, XXL-LS, and QL) occurred 0.22 Ma ago, coincident with the first population bottleneck of the two species caused by the Penultimate Glaciation. Next, EH-GLG and XXL-LS-QL diverged 0.104 Ma ago. The divergence may have resulted from the widely unsuitable habitat located in the Daxueshan and SLL Mountains (21). Last, XXL-LS and QL diverged 26 ka ago, which was most likely caused by the Last Glacial Maximum. After the population divergence, MH, EH-GLG, and QL suffered from population decline, whereas XXL-LS experienced population growth. Asymmetrical gene flow was detected between adjacent divergent populations (Fig. 3B). After the early divergence between the two species, more gene flow occurred from the Chinese red panda to the Himalayan red panda. Regardless of historical or current gene flow, EH-GLG seemed to be the source population of gene flow with more gene flow into other adjacent populations, among the four genetic populations (Fig. 3B). This implies that EH-GLG might be the historical dispersal source of red pandas. TreeMix analysis also detected significant gene flow from EH-GLG to XXL-LS (Fig. 3C and fig. S6), consistent with the Fastsimcoal2 result.

Genomic variation, linkage disequilibrium, and genetic load Whole-genome variation analysis revealed that EH-GLG had the highest genetic diversity (π = 6.994 × 10−4, θ w = 5.271 × 10−4), whereas the Himalayan red panda (MH) had the lowest genetic diversity (π = 3.523 × 10−4, θ w = 2.428 × 10−4) (Fig. 4A and table S10). Y-SNP and mitochondrial genomic variations also showed that the Himalayan red panda (MH) had the lowest genetic variations (Fig. 4A and table S10). Genome-wide linkage disequilibrium (LD) analysis demonstrated that the Himalayan red panda (MH) had higher level of LD and slower LD decay with a reduced R2 correlation coefficient becoming stable at a distance of approximately 100 kb, whereas the populations of the Chinese red panda exhibited rapid LD decay with a reduced R2 becoming stable at a distance of approximately 40 kb (Fig. 4B). The genomic variations and LD patterns imply different demographic histories of the two species and, in particular, reflect the genetic impacts of long-term population bottlenecks in the Himalayan red panda. Fig. 4 Genetic variation, genetic load, and natural selection of red pandas. (A) Genetic variations (nucleotide diversity) of different species and populations based on whole-genome SNPs, mitochondrial genomes, and Y chromosome SNPs. (B) LD of the four populations. (C) Ratios of homozygous derived deleterious or LoF variants to homozygous derived synonymous variants for different populations. The horizontal bars denote population means. (D) Distribution of θ π ratios (non-MH/MH) and Z(F ST ) values. Data points located to the left of the left vertical dashed lines and the right of the right vertical dashed lines (corresponding to the 5% left and right tails of the empirical θ π ratio distribution, respectively) and above the horizontal dashed line [the 5% right tail of the empirical Z(F ST ) distribution] were identified as selected regions for the MH (the Himalayan red panda, green points) and non-MH (the Chinese red panda, blue points) populations. We further analyzed the relationship between demographic history and genetic loads carried by different red panda populations, as deleterious variations should be removed more efficiently in larger populations (22, 23). We investigated the distributions of four types of variations [loss of function (LoF), deleterious, tolerated, and synonymous mutations] in protein-coding genes. We found that the ratios of homozygous derived deleterious or LoF variants to homozygous derived synonymous variants were higher in the Himalayan red panda (MH) than in the Chinese red panda; by contrast, the ratios of nonhomozygous derived deleterious or LoF variants to nonhomozygous derived synonymous variants were comparable between the two species (Fig. 4C). This genetic load pattern showed that the Himalayan red panda experiencing long-term population bottlenecks carried more homozygous LoF and deleterious mutations and thus suffers a higher risk of continuing population decline.