Adaptive Genetic Variation in Black‐and‐White Snub‐Nosed Monkeys (Rhinopithecus bieti): Low Diversity and the Role of Balancing Selection

TL;DR

This study examines the low genetic diversity in black-and-white snub-nosed monkeys and its implications for conservation.

Contribution

The study provides new insights into adaptive genetic variation and balancing selection in an endangered primate species.

Findings

Low heterozygosity and polymorphism were observed in MHC class I exons of Rhinopithecus bieti.

Positive selection and trans-species polymorphism were detected, indicating balancing selection.

Low genetic diversity may reduce the species' ability to adapt to environmental changes.

Abstract

Genetic variation provides the raw material for natural selection, enabling species to maintain adaptive potential, respond to environmental changes, and resist pathogens. Reduced genetic diversity can severely compromise long‐term viability, particularly in small, isolated populations prone to inbreeding, genetic drift, and restricted gene flow—a vicious cycle known as the “extinction vortex”. Assessing genetic diversity in threatened species is therefore critical for effective conservation strategies. The black‐and‐white snub‐nosed monkey (
Rhinopithecus bieti
) is an Endangered primate that has experienced significant population decline and habitat fragmentation, raising concerns about its genetic diversity. We utilized major histocompatibility complex (MHC) class I genes, whose encoded proteins recognize antigens central to immune responses, to assess the adaptive genetic diversity of a semi‐provisioned subpopulation of this species. Species‐specific multi‐locus primers targeting exons 2 and 3 of MHC class I genes were designed using published
R. bieti
whole‐genome sequences. Amplicon‐based next‐generation sequencing was employed to genotype these exons in the studied subpopulation inhabiting Baima Snow Mountain National Nature Reserve, Yunnan, China. A total of 16 MHC class I sequences (7 exon 2 sequences and 9 exon 3 sequences) were identified from 47 individuals and assigned to 5 loci. Exon 2 exhibited low heterozygosity (H
e = 0.349) and moderate polymorphism (PIC = 0.281), whereas exon 3 showed extremely low heterozygosity (H
e = 0.147) and low polymorphism (PIC = 0.131). In addition, positive selection signatures were detected in both exons, and phylogenetic analyses indicated trans‐species evolutionary patterns in class I loci. These results underscore the role of balancing selection in maintaining adaptive genetic variation. However, low genetic diversity is likely to have diminished the studied subpopulation's capacity to adapt to environmental change, thereby undermining its long‐term viability. This study emphasizes the urgent need to assess adaptive genetic diversity across all
R. bieti
populations in order to develop targeted management strategies. The data generated in this studied subpopulation provide the baseline for comparison.

We investigated the adaptive genetic diversity within a population of Rhinopithecus bieti using amplicon sequencing. From 47 individuals, we identified 16 MHC class I sequences, revealing a low level of adaptive genetic variation. Although evidence for balancing selection, including positive selection and trans‐species polymorphism, was observed, these selective forces did not ameliorate the observed low genetic diversity. Our findings underscore the urgent need to assess adaptive genetic diversity across the species to inform conservation management strategies.

Figures

Haplotype composition of MHC class I exons 2 (a) and 3 (b) in Rhinopithecus bieti.

Multiple sequence alignments of the amino acid sequences deduced from the MHC class I sequences of Rhinopithecus bieti. Dots (.) denote amino acids identical to the first sequence; dashes (−) indicate deletions relative to the first sequence. Asterisks (*) mark inferred antigen‐binding sites (ABSs), plus symbols (+) mark inferred T‐cell receptor (TCR) binding sites (both inferred from HLA equivalents; Maccari et al. 2018). Hashtag symbols (#) denote positively selected sites detected via PAML v4.7 (Yang 2007).

Phylogenetic trees of Rhinopithecus bieti MHC class I exon 2 (a: Bayesian; c: maximum likelihood) and exon 3 (b: Bayesian; d: maximum likelihood) with homologous sequences from closely related species. Branch support values (Bayesian posterior probabilities > 55; UFBoot2 support > 50) are shown at relevant nodes.

Tables

TABLE 1. Sequence diversity of MHC class I genes in Rhinopithecus bieti .

Exon	Nucleotide			Amino acid
	k	D	Pi	k′	D′	Pi′
exon 2	35.000 ± 3.944	0.150 ± 0.068	0.133 ± 0.017	21.762 ± 8.988	0.321 ± 0.144	0.247 ± 0.102
exon 3	37.917 ± 3.975	0.164 ± 0.108	0.139 ± 0.028	21.000 ± 13.119	0.317 ± 0.213	0.233 ± 0.146
p	0.164	0.183	0.183	0.193	0.305	0.226

TABLE 2. Genetic diversity of MHC class I loci in Rhinopithecus bieti .

Exon	Locus	N	A e	A e ′	H o	H e	H e′	PIC
exon 2	loc1 (01/04)	47	1.367	1.357	0.319	0.271	0.263	0.232
	loc2 (02/03)	46	1.941	1.927	0.652	0.490	0.481	0.367
	loc3 (05/06)	47	1.394	1.376	0.340	0.286	0.273	0.243
	Mean		1.567	1.553	0.437	0.349	0.339	0.281
exon 3	loc1 (01/04)	47	1.160	1.167	0.149	0.139	0.143	0.128
	loc2 (02/03)	47	1.136	1.143	0.085	0.121	0.125	0.112
	loc3 (07/08)	47	1.066	1.065	0.064	0.062	0.061	0.060
	loc4 (05/06)	10	1.342	1.364	0.300	0.268	0.267	0.222
	Mean		1.176	1.185	0.150	0.147	0.149	0.131

TABLE 3. Rate of non‐synonymous substitutions (d N) and synonymous substitutions (d S) of MHC class I sequences in Rhinopithecus bieti .

Exon	Sites	N	d N	d S	ω	Z	p
exon 2	ABSs	24	0.339 ± 0.088	0.116 ± 0.058	2.922	2.315	0.011
	Non‐ABSs	64	0.109 ± 0.022	0.103 ± 0.032	1.058	0.333	0.370
	All	88	0.163 ± 0.026	0.106 ± 0.027	1.538	1.652	0.051
exon 3	ABSs	17	0.208 ± 0.074	0.193 ± 0.118	1.078	0.113	0.455
	Non‐ABSs	75	0.142 ± 0.024	0.195 ± 0.048	0.728	1.093	0.138
	All	92	0.153 ± 0.023	0.200 ± 0.043	0.765	1.065	0.144

TABLE 4. Results of likelihood‐ratio tests for codon evolution in MHC class I sequences of Rhinopithecus bieti .

Exon	Models compared	df	Test statistic	Significance (p)
exon 2	M2a vs. M1a	2	16.460	< 0.001
	M3 vs. M0	4	43.134	< 0.001
	M8 vs. M7	2	16.586	< 0.001
exon 3	M2a vs. M1a	2	10.866	< 0.001
	M3 vs. M0	4	53.126	< 0.001
	M8 vs. M7	2	15.194	< 0.001

TABLE 5. Codon evolution model and sites of positive selection of MHC class I genes of Rhinopithecus bieti (maximum likelihood analysis).

Exon	Model	k	Log‐likelihood	Estimates parameters	Positively selected sites
exon 2	M0(one ratio)	1	−766.098	ω = 0.681	None
	M1a(neutral)	2	−753.714	P 0 = 0.506 (P 1 = 0.493) ω 0 = 0.050, ω 1 = 1.000	Not allowed
	M2a(selection)	4	−745.484	P 0 = 0.411, P 1 = 0.574 (P 2 = 0.013) ω 0 = 0.000, ω 1 = 1.000, ω 2 = 38.980	23A, 69A, 72T **
	M3(discrete)	5	−744.531	P 0 = 0.586, P 1 = 0.401 (P 2 = 0.012) ω 0 = 0.148, ω 1 = 1.882, ω 2 = 57.757	4M *, 5R*, 8Y, 10A **, 11V, 17A, 20Y, 21F, 22F**, 23A **, 25G, 34R, 40E, 42P, 47R, 48A, 49P, 51V, 55G, 61Q *, 62E , 64Q **, 65N **, 66M *, 67K, 68T, 69A **, 70T, 71Q, 72T **, 73Y **, 74R, 75G, 76N, 78R, 79T, 81L*, 82R
	M7(beta)	2	−753.802	p = 0.024, q = 0.021	Not allowed
	M8(beta & ω)	4	−745.509	P 0 = 0.986 (P 1 = 0.013) p = 0.007, q = 0.005, ω = 39.612	23A, 65N , 69A, 72T **
exon 3	M0(one ratio)	1	−913.652	ω = 0.399	None
	M1a(neutral)	2	−893.261	P 0 = 0.699 (P 1 = 0.300) ω 0 = 0.109, ω 1 = 1.000	Not allowed
	M2a(selection)	4	−887.828	P 0 = 0.722, P 1 = 0.183 (P 2 = 0.094) ω 0 = 0.163, ω 1 = 1.000, ω 2 = 7.066	4I *, 6E, 59A, 60R **, 61A *, 64Q, 65H, 72T
	M3(discrete)	5	−887.089	P 0 = 0.267, P 1 = 0.614 (P 2 = 0.118) ω 0 = 0.000, ω 1 = 0.377, ω 2 = 5.850	4I **, 6E *, 23E, 58A , 59A, 60R **, 61A **, 64Q, 65H, 72T, 79R, 91A
	M7(beta)	2	−894.940	p = 0.261, q = 0.402	Not allowed
	M8(beta & ω)	4	−887.343	P 0 = 0.889 (P 1 = 0.110) p = 0.773, q = 1.940, ω = 6.236	4I *, 6E, 58A , 59A, 60R **, 61A *, 64Q, 65H, 72T *

Full text

Introduction

1

Genetic diversity is a crucial indicator of a population's ability to adapt to environmental challenges and selective pressures associated with climate change and emerging pathogens (Reed and Frankham 2003). Species with low genetic variation encounter inbreeding risks, increasing the likelihood of expressing recessive deleterious alleles, eroding heterozygote advantages, and lowering overall population fitness (Hedrick and Garcia‐Dorado 2016; Frankham 2015; Kardos et al. 2021; Willi et al. 2022). Consequently, increased genetic diversity enhances both survival and evolutionary potential (Keller and Waller 2002). In addition, the conversion of natural forests into anthropogenic landscapes has given rise to patterns of habitat fragmentation resulting in small isolated local populations characterized by restricted gene flow, inbreeding, and genetic drift (Frankham et al. 2017). This self‐reinforcing cycle further diminishes genetic diversity, accelerates population decline, intensifies isolation, and ultimately traps populations in an extinction vortex (Nabutanyi and Wittmann 2021). Therefore, a more complete understanding of the genetic diversity of threatened (Vulnerable, Endangered, and Critically Endangered) species is essential for designing effective conservation strategies that promote population persistence and recovery (Frankham et al. 2002; Kohn et al. 2006; Allendorf et al. 2010).

The black‐and‐white snub‐nosed monkey ( Rhinopithecus bieti ) is an Endangered species of colobine primate (colobines are characterized by an enlarged, low acid, multichambered stomach and high microbacterial diversity that facilitates the breakdown of difficult‐to‐digest complex carbohydrates and plant toxins; Liu et al. 2022) endemic to China (Long et al. 2020). This species inhabits coniferous, mixed coniferous‐broadleaf, and broadleaf forests at elevations ranging from 3000 to 4700 m, the highest elevational distribution for a non‐human primate (Liu et al. 2007, 2009; Zhou et al. 2014). In this regard, black‐and‐white snub‐nosed monkeys face significant challenges associated with thermoregulation during long, cold winters and hypoxia at these high elevations (Zhou et al. 2014). Moreover, the distribution of black‐and‐white snub‐nosed monkeys is restricted to a narrow corridor between the Lancang (Mekong) and Jinsha (Yangtze) Rivers in the Yunling Mountains, where populations are fragmented across anthropogenic landscapes that include farmland, roads, and human settlements (Xia et al. 2022; Cui et al. 2008; Kuang et al. 2020). Microsatellite analyses have shown that habitat fragmentation has significantly limited gene flow across the species range, resulting in five genetically distinct population clusters (Liu et al. 2009). While the species' population size has increased from < 2000 in 1994 to approximately 3800 in 2021 in response to conservation efforts (Yu et al. 2016; Liu et al. 2008; Zhao et al. 2019; Li et al. 2024), assessments reveal that genome‐wide heterozygosity (0.034%) in R. bieti is the lowest among the five extant Rhinopithecus species ( R. bieti , R. roxellana , R. brelichi , R. strykeri , and R. avunculus ) (Kuang et al. 2020). In this species, a restricted distribution, small population size, and prolonged geographic isolation appear to have accelerated the loss of genetic diversity and inbreeding depression, posing significant threats to the species' survival (Ralls et al. 2018).

Accurate genetic diversity assessments face sampling challenges, particularly for populations that exploit large home ranges and inhabit mountainous terrain. In this regard, non‐invasive sampling methods (e.g., fecal or hair collection) minimize species disturbance, but typically yield degraded DNA prone to contamination. This can limit genomic analyses. Consequently, studies using these samples primarily rely on neutral markers (e.g., SSRs [simple sequence repeats], mitochondrial D‐loop) for diversity and structure assessments (Andrews et al. 2016; De Barba et al. 2010; Carroll et al. 2018). These methods, however, cannot detect environmentally driven adaptive variation (Neefjes et al. 2011). In this regard, the major histocompatibility complex (MHC) gene family serves as a preferred molecular marker for comprehensively assessing adaptive genetic diversity in wildlife populations for which hair and fecal samples are most easily obtained.

MHC genes encode proteins essential for pathogen recognition and antigen presentation in vertebrates. Class I and II molecules exhibit distinct functions: Class I presents endogenous antigens (e.g., viral proteins) to initiate cytotoxic responses, while class II presents exogenous antigens (e.g., bacterial proteins) to trigger humoral immunity (Neefjes et al. 2011; Trowsdale and Knight 2013; Hulpke and Tampé 2013). Their antigen‐binding domains, which in class I are encoded by exons 2–3 and in class II are encoded by exon 2, exhibit extreme polymorphism, reflecting adaptations associated with pathogen‐mediated coevolutionary arms races (Sutton et al. 2011; Herdegen et al. 2014). For adaptive loci such as MHC, natural selection, particularly pathogen‐mediated balancing selection (e.g., negative frequency‐dependent selection or heterozygote advantage), typically targets these specific antigen‐binding domains and tends to preserve adaptive variation. This results in a level of polymorphism that is generally higher than that of neutral loci, whose polymorphisms are predominantly shaped by demographic history and genetic drift. Therefore, comparing the diversity patterns between adaptive loci (e.g., MHC) and neutral loci (e.g., genome‐wide single nucleotide polymorphisms or microsatellites) identify the strength and mode of natural selection acting on specific genes (Niskanen et al. 2014; Herdegen‐Radwan et al. 2021). Notably, in mammals, MHC class I genes have experienced frequent intra‐species gene duplication (resulting in paralogous genes), and evolve at a faster rate than MHC class II genes, which are primarily characterized as orthologous genes (Takahashi et al. 2000; Piontkivska and Nei 2003; Kuduk et al. 2012). A relatively rapid evolutionary rate, drives MHC class I genes to accumulate substantial adaptive genetic variation. This results in marked interspecific divergence, thereby limiting cross‐species research that uses MHC class I genes (Sommer 2005; Radwan et al. 2020).

A recent study that explored adaptive genetic diversity in R. bieti using five MHC class II loci reported relatively high heterozygosity and moderate polymorphism (PIC = 0.383, H e = 0.500; Yan et al. 2024). However, the genetic variation and maintenance mechanisms of MHC class I genes in this species remain uncharacterized. In order to elucidate the genetic diversity of MHC class I genes in a semi‐provisioned subpopulation of black‐and‐white snub‐nosed monkeys, we utilized whole‐genome sequencing data (ASM169854v2; Yu et al. 2016) and designed species‐specific multi‐locus amplification primers targeting MHC class I exons 2 and 3. These exons encode essential functional domains that bind peptides to T cells (Wieczorek et al. 2017). Subsequently, amplicon sequencing was employed to explore the genotypes of these exons.

The specific objectives of this study are to (i) assess the genetic heterozygosity and polymorphism of MHC class I genes in a semi‐provisioned subpopulation of R. bieti ; and (ii) identify different selective agents (including positive selection and trans‐species evolution) that maintain MHC class I diversity in this subpopulation of this Endangered primate species.

Materials and Methods

2

Study Site and Sampling

2.1

This study focuses on an isolated population of R. bieti inhabiting the Xiangguqing area within the Baima Snow Mountain National Nature Reserve, Yunnan Province, China. With an estimated population size of 430–480 individuals, the Xiangguqing population represents one of the largest among the 24 known natural populations of this species (Li and Kuang 2025). Its habitat exhibits a pronounced altitudinal gradient and a set of diverse vegetation types, making it ecologically representative of the species' distributional range (Li et al. 2023). To facilitate scientific monitoring, the Xiangguqing population was split by reserve managers into a smaller semi‐provisioned subpopulation and a larger non‐provisioned subpopulation (Zhu et al. 2016). The home ranges of the two bands overlap, with individuals occasionally transferring between them.

The samples for this study were collected between September and December 2020 from the semi‐provisioned subpopulation, which receives supplementary food throughout the year (approximately 19.6% of their diet; Li et al. 2023), and forages independently for the remainder of their diet. During the study period, the semi‐provisioned subpopulation contained 57–62 individuals, representing approximately 11.88%–14.42% of the entire Xiangguqing population. Long‐term monitoring of this subpopulation has enabled researchers to achieve precise individual identification and kinship assessment, thereby facilitating targeted sampling at the individual level.

Non‐invasive methods were employed to collect fecal samples. A total of 60 fecal samples from 47 individuals were collected, which included 10 one‐male units (OMUs) and one all‐male unit (AMU) within the subpopulation. The sampled individuals included 27 adults, 9 subadults, and 11 juveniles. All fecal samples were collected within 15 min of being voided. Disposable gloves and masks were worn during sample collection to minimize the risk of contamination. The fecal samples were preserved in DETs solution (20% DMSO, 0.25 M EDTA, 100 mM Tris–HCl, saturated NaCl, pH 8.0) at a ratio of 1:3 (sample: buffer). All samples were subsequently stored at −80°C for long‐term preservation.

In this research, no individuals were captured and no biological tissues were collected. During sampling, individuals were attracted through food provisioning to enable tracking, identification, and fecal sample collection. No harm was caused to any individual during the process. The research protocol was approved by the Ethical Committee of Northwest University, Xi'an, China (Resolution No.: AWC‐20210402M).

Molecular Experiments

2.2

DNA Extraction

2.2.1

Genomic DNA was extracted from fecal samples using the QIAamp DNA Stool Mini Kit (QIAGEN, Germany), according to the manufacturer's protocol.

Primer Development

2.2.2

To comprehensively characterize MHC class I gene sequences, a multi‐primer amplification strategy was employed. Reference sequences for MHC class I exons 2 and 3 from rhesus macaques ( Macaca mulatta ) (GenBank accession numbers: NW005093093, NW021160161) were aligned against the whole‐genome sequence of R. bieti using Genome BLAST. We prioritized the earlier genome assembly (ASM169854v2; Yu et al. 2016) over the more recent one (PGDP_RhiBie; Wu et al. 2023) for two primary reasons: first, it has undergone manual review, annotation, and optimization by RefSeq, thus achieving superior annotation quality—a feature absent in the new assembly; second, it exhibits better assembly continuity (103,881 scaffolds vs. 377,151 scaffolds). This combination of high‐quality annotation and robust continuity provided a more reliable foundation for our analysis. We identified 18 putative MHC class I sequences. Pseudogenes containing frameshift mutations or premature stop codons were excluded, leaving four candidate functional MHC class I sequences (FG1–FG4). Species‐specific primers (30 pairs) targeting exons 2 and 3 of these four sequences were designed with DNASTAR software (https://www.dnastar.com/), ensuring selective amplification of functional MHC class I loci while blocking pseudogenes. Initial PCR amplification and Sanger sequencing were employed to eliminate primers with poor amplification efficiency, non‐specific products, or unexpected fragment lengths. Seven primer pairs were retained for population‐level amplification (Table S1).

To meet the sequencing technical requirements (Illumina NovaSeq 6000, 250 bp paired‐end reads) and ensure an appropriate amplicon length (recommended: 300–500 bp), a maximum optimization threshold of 450 bp was established. Nested primers were designed for three primer pairs (Ex2_FG3, Ex3_FG1, Ex3_FG2) with amplicons approaching or exceeding this threshold (Table S1). Figure S1 schematically illustrates the primer‐binding regions designed on the four template sequences.

For sample multiplexing during amplicon sequencing, barcode sequences were appended to the four specific primers (Ex2_FG1, Ex2_FG4, Ex3_FG3, and Ex3_FG4) and three nested primers, creating fusion primers. Barcode incorporation primers were synthesized by Sangon Biotech (Shanghai) Co. Ltd.

PCR Amplification

2.2.3

The seven selected primer pairs were employed in combination to amplify MHC class I exons 2 and 3 from R. bieti fecal DNA samples. For primer pairs requiring single‐round amplification (Ex2_FG1, Ex2_FG4, Ex3_FG3, Ex3_FG4), PCR reactions were prepared in 25‐μL volumes containing: 12.5 μL 2× HotStarTaq Plus Master Mix (QIAGEN, Germany), 0.5 μL each of forward and reverse primers (10 μmol/L), 2 μL DNA template, 0.25 μL bovine serum albumin (BSA; TAKARA, Japan), 2.5 μL 10× CoralLoad Concentrate (QIAGEN, Germany), and 6.75 μL nuclease‐free water. Thermal cycling conditions included an initial denaturation at 95°C for 5 min, followed by 35 cycles of 95°C for 30 s, primer‐specific annealing (individual T a in Table S1) for 30 s, and extension at 72°C for 40 s, with a final extension at 72°C for 10 min.

As for primer pairs requiring nested PCR (Ex2_FG3, Ex3_FG1, Ex3_FG2), the first‐round reactions were prepared in 12.5‐μL volumes comprising: 6.25 μL 2× HotStarTaq Plus Master Mix (QIAGEN, Germany), 0.25 μL each of outer primers (10 μmol/L), 1 μL DNA template, 0.125 μL bovine serum albumin (BSA; TAKARA, Japan), 1.25 μL 10× CoralLoad Concentrate (QIAGEN, Germany), and 3.375 μL nuclease‐free water. The thermal cycling conditions were identical to those used for single‐round amplification. The nested PCR reactions (25 μL) contained: 12.5 μL 2× Taq Master Mix (Runde, China), 0.5 μL each of inner primers (10 μmol/L), 2.5 μL of the diluted primary PCR product, and 9 μL nuclease‐free water. The thermal profile was modified as follows: initial denaturation at 95°C for 5 min, followed by 30 cycles of 95°C for 30 s, primer‐specific annealing (see Table S1) for 30 s, and extension at 72°C for 30 s, with a final extension at 72°C for 10 min.

Amplicon Sequencing

2.2.4

All amplified products were purified using the UE PCR Clean‐Up Kit (UElandy, China), quantified using the Qubit High‐Sensitivity Kit, diluted to 10 ng/μL, pooled, and sequenced as 250 bp paired‐end reads on the Illumina NovaSeq 6000 platform (Novogene, Beijing, China).

Data Analysis

2.3

MHC Genotyping

2.3.1

Raw FASTQ files were demultiplexed using sample‐specific barcode sequences with exact barcode matching through the atlas‐utils demux tool (v0.0.1; https://github.com/jameslz/biostack‐suits). Quality filtering was conducted using Trimmomatic (v0.39) through the following steps: (i) Removal of bases at the 5′ and 3′ ends with quality scores < 2; (ii) Truncation of reads using a 4‐base sliding window when the average quality score within the window fell below 15; (iii) Discarding of reads shorter than 36 bp after truncation. Paired‐end reads were merged using the USEARCH fastq_mergepairs command (v11.2.64; http://drive5.com) with default parameters (Edgar and Flyvbjerg 2015). Only high‐quality merged reads were retained, and unmerged reads were discarded. Merged reads were aligned to primer sequences, and those with > 2 nucleotide mismatches were excluded. Primer regions were trimmed from matched reads. Non‐redundant unique reads were extracted, and their abundance was quantified using atlas‐utils unique (v0.0.1; https://github.com/jameslz/biostack‐suits). Potential PCR artifacts and chimeras were removed using the methods described by Santos et al. (2017) and Sommer et al. (2013). To ensure the reliability of identified variants, each DNA sample was subjected to at least two independent PCR amplifications, which were performed on separate days using independently prepared reaction mixtures. Only those MHC class I sequences consistently detected across all PCR replicates were retained as R. bieti alleles for subsequent analysis. The isolated sequences were then subjected to BLAST alignment against the R. bieti genome to evaluate their homology with the MHC class I sequences of this species.

R. bieti MHC class I sequences were assigned to specific loci using MHC‐TYPER v1.0 (Huang et al. 2019). This assignment was achieved by leveraging two types of genetic information—allele co‐occurrence patterns and genotype frequencies—and integrating error penalty mechanisms, likelihood and Bayesian information criterion (BIC) calculations, and a simulated annealing algorithm to identify the optimal allele configuration across duplicated loci. Linkage disequilibrium (LD) among loci and haplotype construction were analyzed using SHEsis (http://analysis.bio‐x.cn/myAnalysis.php; Shi and He 2005), with analyses based on D′ values.

Genetic Diversity

2.3.2

Nucleotide sequences were translated into amino acid sequences using MEGA v7.0 (Kumar et al. 2016). Antigen‐binding sites (ABSs) and T‐cell receptor (TCR) binding sites in MHC class I sequences were inferred based on homology with human HLA class I sequences (Maccari et al. 2018), followed by multiple sequence alignment of the amino acid sequences. Nucleotide diversity (Pi) and the average number of nucleotide differences (k) for R. bieti MHC class I sequences were calculated using DnaSP v5.10.1 (Librado and Rozas 2009). Sequence divergence was quantified by calculating average nucleotide genetic distance (D) under the Tamura‐Nei model (Tamura and Nei 1993) in MEGA v7.0. To evaluate the variation characteristics of amino acid sequences, the average amino acid genetic distance (D′) was calculated using the JTT matrix (Jones et al. 1992) in MEGA v7.0. The mean proportion of differing sites across all sequence pairs was used to characterize amino acid diversity (Pi′) with the p‐distance model (Nei and Kumar 2000), and the average number of amino acid differences (k′) among all sequence pairs was calculated using the ‘Number of differences’ model.

For each MHC class I locus (exons 2 and 3), the number of alleles (N), polymorphism information content (PIC), observed heterozygosity (H o), and expected heterozygosity (H e) were computed using Cervus v3.0.7 (Kalinowski et al. 2007). The effective number of alleles (A e) was calculated using GenAlEx 6.5 (Peakall and Smouse 2012). Additionally, to mitigate potential kinship bias among the 47 individuals that might confound population genetic analyses, we conducted an unbiased estimation of H e and A e using the N e (effective population size) estimator proposed by Nomura (2008). This approach identifies putative non‐relatives based on the pairwise kinship coefficient θ^xx′l¯, which is estimated via Weir's (1996) method. A 90% threshold was set for the selection of putative non‐relatives, corresponding to the number of social units (10 OMUs) in the sampled subpopulation.

Selective Pressure Analysis

2.3.3

The ω value (ratio of non‐synonymous to synonymous substitution rates, d N /d S) was calculated for antigen‐binding sites (ABSs), non‐antigen‐binding sites (non‐ABSs), and all amino acid sites in exons 2 and 3 using the Nei‐Gojobori method with a Jukes‐Cantor correction (Nei and Gojobori 1986) in MEGA v7.0. One thousand bootstrap resampling replicates were used to estimate standard errors. To detect amino acid residues under positive selection and estimate ω values, the CODEML program in PAML v4.7 (Yang 2007) was used. Six codon‐substitution models were tested under hypotheses of positive selection, neutral evolution, and purifying selection. The M0 (single‐ratio model) assumes a uniform ω across all sites; M1a (nearly neutral model) allows sites with ω ≤ 1 (neutral + purifying selection); M2a (positive selection model) extends M1a by adding a class of sites with ω > 1; M3 (discrete model) permits two or more freely estimated ω classes; M7 (β‐distribution model) fits ω values to a β‐distribution (0 ≤ ω ≤ 1); and M8 (β + ω model) extends M7 by adding a class of sites with ω > 1. Likelihood‐ratio tests (LRTs) were performed to compare models allowing for different selection intensity among sites (M2a vs. M1a, M3 vs. M0, and M8 vs. M7). Tests to evaluate the significance of positive selection were conducted under a β‐distribution framework.

Phylogenetic Analysis

2.3.4

Homologous sequences of MHC class I exons 2 and 3 from a total of 7 primate species were obtained from NCBI (https://www.ncbi.nlm.nih.gov; Table S2), supplemented with 58 unpublished MHC Class I sequences from R. roxellana , a congener of the black‐and‐white snub‐nosed monkey. To root the phylogenetic tree, MHC class I allele Tube‐W22 (GenBank accession: JQ353673.1)—derived from the northern tree shrew ( Tupaia belangeri ), the closest extant sister group of primates—served as the outgroup. All sequences were aligned using MAFFT v7.505 (Katoh and Standley 2013), and unreliable alignment regions containing substantial gaps or low‐conservation sites were removed using trimAl v1.5.rev0 (Capella‐Gutiérrez et al. 2009). The best‐fitting nucleotide substitution models for phylogenetic reconstruction were determined separately for maximum likelihood and Bayesian analyses using ModelFinder (Kalyaanamoorthy et al. 2017). The maximum likelihood tree was constructed in IQ‐TREE v3.0.1 (Wong et al. 2025) under the best‐fit model, with node support assessed by 1000 UFBoot2 replicates (Minh et al. 2013) and 1000 SH‐aLRT tests. Bayesian inference was performed in MrBayes v3.2.7 (Ronquist et al. 2012) using four Markov chains (temp = 0.2) run for 1000,000 generations with sampling every 1000 generations and 25% burn‐in discarded. Node support was quantified by posterior probabilities. Final trees were visualized and annotated in iTOL v6 (Letunic and Bork 2021).

Results

3

MHC Allele Assignment

3.1

Genotyping was successfully achieved for 47 individuals, which represents 75.8% to 82.5% of the entire semi‐provisioned subpopulation. Amplicon sequencing identified 16 MHC class I sequences, including 7 exon 2 sequences and 9 exon 3 sequences. For exon 2, the number of alleles per individual ranged from 4 to 7 (mean = 5.3). Exon 3 exhibited a similar range of alleles per individual (4 to 7 alleles; mean = 4.6). BLAST alignment was performed against the R. bieti genome to verify the homology of the isolated sequences with R. bieti MHC class I sequences. All obtained sequences were confirmed to be derived from R. bieti , with their sequence similarity to published R. bieti genomic fragments (TaxID: 61621) ranging from 96% to 100% (Table S3).

Using MHC‐TYPER v1.0, the 16 MHC class I sequences were assigned to five loci based on allele‐locus mapping and the inherent exon 2–3 linkage within MHC loci: loc1 (Rhbi‐E201/04*; Rhbi‐E301/04*), loc2 (Rhbi‐E202/03*; Rhbi‐E302/03*), loc3 (Rhbi‐E205/06*; Rhbi‐E307/08*), loc4 (Rhbi‐E305/06*), and loc5 (Rhbi‐E207*; Rhbi‐E309*). Haplotype reconstruction via SHEsis identified 10 exon 2 haplotypes and 11 exon 3 haplotypes (Figure 1).

Haplotype composition of MHC class I exons 2 (a) and 3 (b) in Rhinopithecus bieti.

Genetic Variation

3.2

The amino acid sequence alignment of the 16 R. bieti MHC class I sequences is presented in Figure 2. Sequence diversity analysis of exons 2 and 3 using DnaSP v5.10.1 revealed that, at the nucleotide level, the mean diversity parameters for exon 2 (k = 35.000 ± 3.944, D = 0.150 ± 0.068, Pi = 0.133 ± 0.017; Table 1) were numerically lower than those for exon 3 (k = 37.917 ± 3.975, D = 0.164 ± 0.108, Pi = 0.139 ± 0.028; Table 1). However, Mann–Whitney U tests showed that these differences were statistically non‐significant (Table 1). In contrast, at the amino acid level, the mean diversity parameters for exon 2 (k′ = 21.762 ± 8.988, D′ = 0.321 ± 0.144, Pi′ = 0.247 ± 0.102; Table 1) were numerically higher than those for exon 3 (k′ = 21.000 ± 13.119, D′ = 0.317 ± 0.213, Pi′ = 0.233 ± 0.146; Table 1). These differences were also non‐significant (Table 1).

Multiple sequence alignments of the amino acid sequences deduced from the MHC class I sequences of Rhinopithecus bieti. Dots (.) denote amino acids identical to the first sequence; dashes (−) indicate deletions relative to the first sequence. Asterisks () mark inferred antigen‐binding sites (ABSs), plus symbols (+) mark inferred T‐cell receptor (TCR) binding sites (both inferred from HLA equivalents; Maccari et al. 2018). Hashtag symbols (#) denote positively selected sites detected via PAML v4.7 (Yang 2007).*

TABLE 1: Sequence diversity of MHC class I genes in Rhinopithecus bieti .

When calculating genetic diversity parameters, monomorphic loci—defined as those in which all individuals within the sampled subpopulation carry an identical allele—exhibited no genetic variability and thus contributed no informative data for a diversity assessment. Their inclusion in the analysis may even interfere with the accurate interpretation of results. Therefore, the monomorphic locus loc5 (Rhbi‐E207*; Rhbi‐E309*) was excluded from the genetic diversity calculations. Analyses conducted using Cervus v3.0.7 showed that exon 2 of the studied subpopulation exhibited low heterozygosity and moderate polymorphism (H e = 0.349, PIC = 0.281; Table 2). In contrast, exon 3 showed extremely low heterozygosity and low polymorphism (H e = 0.147, PIC = 0.131; Table 2). Furthermore, unbiased estimates of genetic diversity parameters, after accounting for kinship effects, indicated that both exon 2 and exon 3 in the studied subpopulation exhibited low expected heterozygosity (exon 2: H e ′ = 0.339; exon 3: H e ′ = 0.149). In addition, the genetic diversity estimates for both exons were consistent between analyses with and without kinship correction.

TABLE 2: Genetic diversity of MHC class I loci in Rhinopithecus bieti .

Positive Selection

3.3

The selection parameter ω (d N /d S) was calculated using MEGA v7.0 for the amino acid sites in exons 2 and 3 of MHC class I, including ABSs, non‐ABSs, and all sites (Table 3). For ABSs in exon 2, ω significantly exceeded 1 (ω = 2.922, p = 0.011), indicating that there was strong positive selection at these sites in R. bieti 's MHC class I sequences. And, although non‐ABSs and all amino acid sites in exon 2, along with ABSs in exon 3, displayed ω > 1, these values were non‐significant. Non‐ABSs and all amino acid sites in exon 3 showed values of ω < 1 and were non‐significant.

TABLE 3: Rate of non‐synonymous substitutions (d N) and synonymous substitutions (d S) of MHC class I sequences in Rhinopithecus bieti .

Selection pressure analysis using PAML v4.7 (CODEML program) revealed that under the Akaike Information Criterion (AIC), the M2a (positive selection model), M3 (discrete model), and M8 (β + ω model) better fit the evolutionary patterns of MHC gene evolution than the other three models (Table 4). Site‐specific tests identified amino acid residues under positive selection with posterior probabilities exceeding 95% (Table 5). For exon 2, M2a detected 3 positively selected residues (including 72T under strong selective pressure), M3 identified 38 residues (13 significant, including 72T), and M8 confirmed 4 residues (notably 72T). For exon 3, M2a identified 8 residues (e.g., 4I, 60R, 61A), M3 detected 12 residues (e.g., 4I, 6E, 60R, 61A), and M8 confirmed 9 residues (e.g., 4I, 60R, 61A, 72T). Most positively selected residues were associated with antigen binding or T‐cell receptor interaction.

TABLE 4: Results of likelihood‐ratio tests for codon evolution in MHC class I sequences of Rhinopithecus bieti .

TABLE 5: Codon evolution model and sites of positive selection of MHC class I genes of Rhinopithecus bieti (maximum likelihood analysis).

Trans‐Species Evolution

3.4

Bayesian and Maximum Likelihood (ML) trees were constructed to investigate the phylogenetic relationships of MHC class I genes among R. bieti and other primates (Figure 3). Trees constructed using the two methods exhibited highly congruent topologies. Phylogenetic analysis revealed a discordance between MHC class I gene relationships and species relationships, with alleles from different species intermixing in the same clusters—consistent with trans‐species polymorphism. No monophyletic clade consisting exclusively of R. bieti alleles was observed. However, specific R. bieti sequences formed well‐supported branches with those of R. roxellana , such as Rhbi‐E204* and Rhbi‐E201* clustering with Rhro‐E204* (Figure 3a,c), and Rhbi‐E302* grouping with Rhro‐E303* (Figure 3b,d). These results suggest high sequence homology between these alleles in these two closely related species.

Phylogenetic trees of Rhinopithecus bieti MHC class I exon 2 (a: Bayesian; c: maximum likelihood) and exon 3 (b: Bayesian; d: maximum likelihood) with homologous sequences from closely related species. Branch support values (Bayesian posterior probabilities > 55; UFBoot2 support > 50) are shown at relevant nodes.

Discussion

4

This study provides the first assessment of MHC class I diversity in the Endangered black‐and‐white snub‐nosed monkey ( Rhinopithecus bieti ), focusing on a semi‐provisioned subpopulation. Our findings indicate a low level of adaptive genetic diversity in this subpopulation, a pattern likely associated with strong genetic drift that is characteristic of small, isolated populations. While signatures of positive selection and trans‐species polymorphism—both lines of evidence for balancing selection—were detected in this study, these selective forces did not ameliorate the low level of genetic diversity in this subpopulation. Instead, genetic drift appears to have had dominant effect, overwhelming the capacity of balancing selection to sustain genetic diversity, thereby leading to the current low level of adaptive genetic variation. By focusing on this Endangered R. bieti subpopulation, this research provides an informative case study for elucidating the genetic effects of habitat fragmentation on population genetic variation, and establishes a robust empirical foundation for evaluating the adaptive genetic status of this species.

Sequence Variation Characteristics

4.1

The present study revealed that, compared with exon 3, the diversity parameters of exon 2 were numerically lower at the nucleotide level but higher at the amino acid level. Although this direct comparison lacked statistical significance (Table 1), the observed trend indicated that while exon 2 is relatively conserved at the nucleotide level, its variation is more likely to induce amino acid substitutions. This result is consistent with the findings of positive selection analysis—the antigen‐binding site (ABSs) region of exon 2 was subject to significant positive selection (ω = 2.922, p = 0.011; Table 3). Specifically, nucleotide variations in exon 2 tend to be nonsynonymous mutations (which directly alter the amino acid sequence), and the accumulation of nonsynonymous mutations is a typical signature of positive selection—whose primary purpose is to enhance amino acid diversity at ABSs in order to adapt to diverse pathogen environments. In contrast, exon 3 exhibited the characteristic of “higher nucleotide diversity but lower amino acid diversity” relative to exon 2, reflecting a tendency to retain the original amino acid sequence and avoid compromising the structural stability of the encoded protein. This characteristic is also consistent with selection analysis results: both ABSs and non‐ABSs were under neutral selection (ABSs: ω = 1.078, p = 0.455; non‐ABSs: ω = 0.728, p = 0.138; Table 3).

To further explore evolutionary differences of MHC class I genes in R. bieti , we compared its diversity parameters with those of its closely related species R. roxellana (golden snub‐nosed monkey) (Dong et al., unpublished data). The results showed that the nucleotide diversity parameters (k, D, Pi) and amino acid diversity parameters (k′, D′, Pi′) of exon 2 were generally comparable between R. bieti and R. roxellana (k: 35.000 ± 3.944 vs. 33.830 ± 3.472; D: 0.150 ± 0.068 vs. 0.142 ± 0.039; Pi: 0.133 ± 0.017 vs. 0.128 ± 0.005; k′: 21.762 ± 8.988 vs. 20.855 ± 4.954; D′: 0.321 ± 0.144 vs. 0.294 ± 0.081; Pi′: 0.247 ± 0.102 vs. 0.237 ± 0.056). In the case of exon3, however, both types of diversity parameters were higher in R. bieti than in R. roxellana (k: 37.917 ± 3.975 vs. 28.082 ± 3.171; D: 0.164 ± 0.108 vs. 0.115 ± 0.042; Pi: 0.139 ± 0.028 vs. 0.105 ± 0.006; k′: 21.000 ± 13.119 vs. 14.280 ± 5.251; D′: 0.317 ± 0.213 vs. 0.199 ± 0.078; Pi′: 0.233 ± 0.146 vs. 0.160 ± 0.059). These results suggest that despite a lower number of alleles (16 vs. 58), R. bieti has retained a considerable level of sequence variation, which may provide a genetic basis for responding to potential pathogen‐mediated selective pressures.

Genetic Diversity

4.2

The number of individuals successfully genotyped at exon 3 loc4 (N = 10) was significantly lower than that at other loci. Due to the limited effective sample size, the genetic diversity parameters for this locus should be interpreted with extreme caution. We speculate that the main reasons for this observation include: (1) the primer‐binding region of this locus exhibits high sequence variability, which may reduce primer binding efficiency to template DNA in some individuals, thereby preventing effective PCR amplification across all surveyed subjects. And (2), the distribution of MHC genes is inherently heterogeneous among individuals. Within the same species, individuals often differ in their number and repertoire of MHC genes they carry, a phenomenon often referred to as “MHC gene copy number variation” or “presence/absence polymorphism”. For example, in humans, certain HLA‐DRB genes (e.g., DRB3, DRB4, DRB5) are present only in a subset of individuals (Robbins et al. 1997). Therefore, the low genotyping success rate observed for loc4 in this study likely reflects that this locus represents a low‐frequency or highly variable MHC gene within the studied subpopulation.

The genome‐wide heterozygosity of Rhinopithecus bieti is extremely low (H e = 0.034%; Kuang et al. 2020), which is consistent with strong genetic drift effects resulting from a small and fragmented population structure, reflecting a significant erosion of the species' genomic genetic diversity under drift. In the present study, the heterozygosity of MHC class I genes in this semi‐provisioned subpopulation (H e: exon 2 = 0.349; exon 3 = 0.147; Table 2) was relatively low, but it did not exhibit the same pattern of extremely low diversity as observed at the genome‐wide level. From an evolutionary perspective, in populations dominated by strong genetic drift, the diversity of functional genes typically decreases synchronously with the genomic background, in the absence of positive selection. However, the finding that MHC genes did not follow the genome‐wide trend toward extremely low diversity in this study indicates that, given the critical role of MHC genes in immune responses, pathogen‐mediated selection pressure has mitigated the effects of genetic drift to a certain extent, thereby preserving key polymorphisms at these adaptive loci. It should be noted that MHC heterozygosity was calculated based on functionally critical coding regions (exons 2 and 3), whereas genome‐wide heterozygosity encompasses a large number of neutral or nearly neutral loci. These two metrics differ inherently in terms of their sources of variation (mutation rates), evolutionary drivers (selective forces), and computational scope, making direct quantitative comparisons inappropriate. Nevertheless, contrasting their diversity patterns under the pressure of genetic drift remains biologically meaningful, providing evidence that pathogen‐mediated balancing selection maintains adaptive variation.

A recent study conducted of this same subpopulation demonstrated that neutral microsatellites exhibited moderate heterozygosity and high polymorphism (SSRs: H e = 0.565, PIC = 0.520; Yan et al. 2024), whereas adaptive MHC class I genes displayed low levels of diversity (with both H e and PIC values). This pattern appears inconsistent with the general expectation—that under pathogen‐mediated selective pressure, adaptive loci typically maintain higher polymorphism than found in neutral loci. It should be noted that the comparison between adaptive and neutral loci in the present study focuses on differences in diversity patterns rather than discrepancies in absolute values. This comparison facilitates the elucidation of evolutionary forces acting on functional genes. Beyond the differences in evolutionary drivers, the observed divergence in diversity patterns between the two types of loci may also be attributable to methodological biases in microsatellite genotyping. Specifically, in the microsatellite diversity analysis (Yan et al. 2024), 8 out of the initial 17 loci were excluded due to amplification failure or the generation of unidentifiable, noisy capillary electrophoresis peaks, with only 9 loci retained for the final analysis. This filtering process likely excluded loci with lower polymorphism, thereby leading to an overestimation of the actual microsatellite diversity in the subpopulation. Future studies will need to validate these findings through expanded sampling, cross‐population comparisons, genome‐wide resequencing of MHC recombination patterns, and unbiased microsatellite screening.

We also found that MHC class I alleles (16 in this study) exceeded class II alleles (11 across DQA1, DQB1, DPB, and DRB loci; Yan et al. 2024) in number within the same subpopulation, aligning with vertebrate‐wide trends. For instance, the IPD‐IMGT/HLA database (April 2025) lists 28,683 HLA class I alleles versus 12,877 class II alleles in humans. Similarly, in a study of olive baboons ( Papio anubis ; n = 22), class I A and B loci harbored 16 and 32 alleles, respectively, compared to 9 DQA and 27 DRB class II alleles (Petersen et al. 2022). When combining individuals from multiple geographic populations, both the Humboldt penguin ( Spheniscus humboldti ) (n = 1 00) and the Magellanic penguin ( Spheniscus magellanicus ) (n = 75) showed higher class I allele counts (236) relative to class II (126) (Sallaberry‐Pincheira et al. 2016). This disparity likely arises from differential selection pressures acting on class I and class II genes. MHC class I molecules primarily present endogenous antigens (e.g., viral peptides from infected cells), requiring rapid adaptation to counteract the high mutation rates characteristic of viruses, particularly RNA viruses (Sanjuán and Domingo‐Calap 2016). Consequently, class I genes maintain elevated polymorphisms to ensure broad‐spectrum recognition of epitopes. In contrast, MHC class II molecules target structurally conserved exogenous antigens (e.g., bacterial proteins) with slower evolutionary rates, leading to MHC class II loci experiencing relatively weaker selective pressures and retaining fewer alleles (van de Weijer et al. 2015; Minias et al. 2019; Piertney and Oliver 2006). Unexpectedly, MHC class I diversity parameters (H e and PIC) in the current study were lower than those of class II loci (H e = 0.500, PIC = 0.383; Yan et al. 2024). This result may stem from past or recent selection, geographic bottlenecks, or genetic drift acting on this subpopulation. Additionally, primer amplification biases may have impeded the detection of class I alleles, thereby systematically underestimating MHC class I diversity.

Compared to the golden snub‐nosed monkey ( Rhinopithecus roxellana ), which has an estimated remaining population of some 22,500 individuals, R. bieti exhibits a smaller population size (~3800 individuals) and lower genetic diversity. Both R. bieti and R. roxellana are endemic to China and are Endangered (EN) (Long et al. 2020; Long and Richardson 2021). These lineages appear to have split from a common ancestor some 1.4 Mya (Kuang et al. 2023), with R. roxellana radiating north and R. bieti entering the Himalayan region. Genome‐wide analyses reveal that R. bieti has the lowest heterozygosity (H e = 0.034%) among the four other extant Rhinopithecus species. For example, heterozygosity in R. bieti is lower than that of R. roxellana (H e = 0.043%), R. avunculus (H e = 0.042%), R. brelichi (H e = 0.069%), and R. strykeri (H e = 0.036%) (Kuang et al. 2020). In terms of adaptive genetic diversity, Dong et al. investigated MHC class I diversity in an R. roxellana population (n = 113) from the Qinling Mountains (Dong et al., unpublished data). These authors identified 58 MHC class I sequences (30 from exon 2 and 28 from exon 3), with individuals carrying 5–15 alleles—far exceeding the 16 MHC class I sequences (4–7 alleles per individual) observed in the current study of this subpopulation of R. bieti . This disparity may result from multiple factors. Population history indicates that R. roxellana experienced two genetic bottlenecks (mid‐Pleistocene ~2.00 Mya and late Pleistocene ~0.10–0.40 Mya) and two expansions (~1.00 Mya and 0.05–0.07 Mya), with the current population distributed across three regions: Minshan‐Qionglai, Qinling, and Shennongjia (Li et al. 2002; Kuang et al. 2023). In contrast, R. bieti populations have declined dramatically (Zhou et al. 2014), with approximately 3800 individuals distributed in five genetically isolated regions confined to a narrow corridor between the Lancang and Jinsha Rivers in the Yunling Mountains (Zhao et al. 2019). The population expansion phases in R. roxellana likely provided opportunities for the accumulation of genetic variation, while smaller population size and geographic isolation appear to have exacerbated genetic drift in R. bieti . Pathogen‐mediated selection pressures may also have played a role. Studies have demonstrated that pathogen diversity and abundance generally increase as climates shift from colder to warmer regions (Cohen et al. 2020; Rohr and Cohen 2020; Altizer et al. 2013). R. roxellana inhabits a range of forested habitats across an altitudinal range of between 1400 and 2800 m. Accounting for the approximate 0.6°C temperature decrease per 100‐m rise in elevation (Barry and Chorley 2009), the mean annual temperature in R. roxellana habitats (10.2°C at 1600 m in the Qinling Mountains; Hou 2018) is slightly higher than in R. bieti habitats (7.5°C at 3200 m in Baima Snow Mountain; Li et al. 2023). In addition, in the northern part of R. bieti 's habitat on Baima Snow Mountain, extreme minimum temperatures can reach −25.4°C, and sub‐zero temperatures (< 0°C) can last for 140 to 200 days annually, while the frost‐free period typically ranges from 80 to 90 days (Li et al. 1982). In contrast, the Qinling Mountains—home to the northernmost populations of R. roxellana —experience an extreme minimum temperature of −19°C (Lu and Li 2006), with the frost‐free period lasting roughly 150 days (Hou 2018). Thus, the results of our research suggest that, based on the observed thermal differences between the habitats of the two species, the evolutionary adaptations of R. bieti to prolonged exposure to low‐temperature environments may potentially reduce pathogen‐mediated selective pressures acting on its MHC diversity. This inference requires further verification, specifically by conducting a correlation analysis among MHC diversity, long‐term climate data, and pathogen distribution data. When combined with the demographic history of R. bieti , the climatic characteristics of its habitat may interact synergistically with the species‐specific traits of small population size and high habitat fragmentation, thereby contributing to the relatively limited MHC variation.

In addition to species' evolutionary history, factors such as population and sample size, sample quality, and genotyping techniques significantly influence allele detection rates. Studies show that increasing sample size and improving sample quality (e.g., using blood samples instead of lower‐quality fecal samples) can enhance allele detection efficiency. Expanding sampling across multiple populations also improves detection rates, as environmental and pathogen‐mediated selective pressures vary between populations of the same species, leading to genetic differentiation (Fumagalli 2013). Furthermore, advancements in genotyping sensitivity increase the likelihood of identifying rare alleles (Salk et al. 2018). In practice, variability in sample types (blood, feces, hair), sampling scale, geographic distribution, and genotyping methods may introduce biases. Therefore, a comprehensive evaluation of these variables is essential when interpreting genetic diversity results.

To date, among snub‐nosed monkey species, detailed MHC class I sequence data are only available for R. bieti (this study) and R. roxellana (Dong et al., unpublished data). Conducting a systematic comparative analysis of MHC diversity across all species of the genus Rhinopithecus (including the critically Endangered R. brelichi , R. strykeri , and R. avunculus ) represents a core priority for future research. This analytical framework will help clarify the roles of phylogeny, demographic history (e.g., bottlenecks, population size), and ecological factors (e.g., pathogen pressure across different habitats) in shaping the polymorphism of adaptive immune genes. Employing standardized, high‐throughput genotyping methods in follow‐up studies is essential to minimize technical biases and enable robust cross‐species comparisons.

Balancing Selection

4.3

Two lines of evidence support the role of balancing selection in shaping MHC diversity in R. bieti . First, our analyses revealed significant positive selection on the antigen‐binding sites (ABSs) of the MHC class I exon 2 (ω = 2.922, p = 0.011; Table 3), as well as on multiple amino acid residues in exons 2 and 3, which are associated with antigen binding and T‐cell receptor interactions (Table 5). These findings suggest that pathogen‐mediated selection pressures promote adaptive genetic variation at functionally critical MHC loci, thereby maintaining polymorphisms (Piertney and Oliver 2006; Zhang et al. 2018; Dearborn et al. 2022). Second, trans‐species polymorphism, which refers to the retention of ancestral alleles across descendant species, provides strong evidence in support of balancing selection (Hedrick 1998; Bernatchez and Landry 2003; Slade et al. 2019). Our phylogenetic analyses demonstrated that MHC class I clades do not align with species boundaries. Instead, they form mixed clusters (Figure 3), indicating that these alleles diverged prior to speciation events and persisted through balancing selection. This trans‐species evolutionary pattern is widely reported in non‐human primates (Pechouskova et al. 2015) and other mammals (Gong et al. 2023), and fish (Geng et al. 2021), further supporting the role of balancing selection in maintaining MHC diversity. Notably, the sequences of exons 2 and 3 from R. bieti formed distinct branches with homologous sequences from R. roxellana (e.g., Rhbi‐E302* with Rhro‐E303*; Figure 3b,d), suggesting recent shared ancestry rather than long‐term independent evolution.

However, despite the detection of signals indicative of historical balancing selection acting on the MHC genes of this species, the MHC diversity of the studied subpopulation remains at a low level. This seemingly contradictory pattern is most likely attributable to the demographic history and neutral evolutionary processes experienced by this subpopulation. In small, fragmented populations, genetic drift can rapidly shape the level of genetic variation (Robertson 1962; Strand et al. 2012). Generally speaking, the efficacy of balancing selection in maintaining diversity is contingent upon population size. In large populations, pathogen‐mediated balancing selection plays a dominant role in shaping genetic variation at MHC loci, rendering these loci more variable than neutral ones (Takahata 1990). In contrast, as population size markedly declines, the effects of genetic drift intensify significantly (Kimura 1983), ultimately making drift the predominant evolutionary force even at loci under the influence of balancing selection. The duration at which the legacy of past balancing selection persists as “a long‐term effective population size” (Takahata 1990) depends on the severity of population contraction, as well as the intensity and nature of historical selection pressures. For instance, studies on black grouse ( Tetrao tetrix ) have shown that balancing selection is the key process influencing MHC loci in large and medium‐sized populations; whereas when populations decline to an extremely small size and become geographically isolated, genetic drift emerges as the dominant factor (Strand et al. 2012).

Given R. bieti 's small population size and distribution across highly fragmented habitats, its effective population size and limited pattern of gene flow are expected to continue to decline. Under such circumstances, strong genetic drift is likely to become the dominant evolutionary force, gradually eroding the MHC polymorphism maintained by historical balancing selection. Therefore, the pattern of diversity observed in this study may reflect the following dynamic process: allelic polymorphism shaped by historical balancing selection (manifested as positive selection and trans‐species polymorphism) has failed to maintain high diversity levels due to intense genetic drift during recent population bottlenecks and isolation events. In conclusion, the current MHC diversity pattern of this studied subpopulation represents the legacy of long‐term balancing selection and strong recent genetic drift, with the latter likely being dominant in the context of its current small population size.

Our study has revealed that the adaptive genetic diversity (MHC class I) of this small and isolated semi‐provisioned primate subpopulation is at a relatively low level. This indicates an elevated risk from environmental change due to reduced adaptive potential, which may compromise its long‐term viability. Therefore, we propose introducing individuals into the subpopulation to implement genetic rescue. However, given that such gene flow may also increase the risk of outbreeding depression, sustained monitoring and careful management remain essential. In future research, it is urgent to integrate genome‐wide data of both neutral and adaptive genetic variation to assess genetic variation across other populations throughout R. bieti 's distributional range. This will not only facilitate an update of its overall genetic profile, but also enable the quantification of adaptive differentiation among distinct populations and the delineation of conservation units (CUs)—particularly those requiring urgent human intervention. In doing so, it will provide a more robust scientific basis for both short‐term and long‐term conservation planning for this Endangered primate species.

Conclusions

5

This study investigated the characteristics of adaptive genetic variation in a semi‐provisioned subpopulation of the Endangered black‐and‐white snub‐nosed monkey ( Rhinopithecus bieti ), using MHC class I genes as markers. The results revealed that the studied subpopulation exhibits low levels of heterozygosity and polymorphism at both exon 2 and exon 3. Evidence of positive selection and trans‐species evolution suggests that balancing selection has historically maintained MHC polymorphism. However, driven by demographic history and neutral evolutionary processes, the level of adaptive genetic diversity in this subpopulation remains low. This implies a decline in their capacity for pathogen resistance and environmental adaptation, thereby posing a severe threat to the long‐term viability of the targeted subpopulation and the broader species.

To support the conservation of R. bieti , we proposed that future research and management practices consider genetic rescue through scientifically regulated gene flow. On the one hand, large‐scale and long‐term genetic monitoring should be carried out, with sampling expanded to systematically collect adaptive genetic diversity data from other populations. In addition, genome‐wide genetic variation data, covering both neutral and adaptive loci, should be integrated to scientifically delineate conservation units (CUs) for this species. Moreover, ecological corridors among fragmented habitats should be restored or maintained to enhance habitat connectivity across the species' distributional range. This will help mitigate the risk of population isolation, facilitate natural gene flow, reduce the degree of inbreeding, and elevate the level of genetic variation.

Author Contributions

Fei Long: data curation (lead), formal analysis (equal), investigation (equal), validation (equal), visualization (lead), writing – original draft (lead), writing – review and editing (equal). Mei‐Jing Ying: data curation (equal), formal analysis (lead), investigation (lead), visualization (equal). Tian‐Qi Shi: formal analysis (equal), investigation (equal), validation (equal). Paul A. Garber: writing – review and editing (equal). Ru‐Liang Pan: writing – review and editing (equal). Jian‐Dong Lai: investigation (equal), resources (equal). Zhi‐Pang Huang: funding acquisition (equal), resources (lead). Bao‐Guo Li: conceptualization (equal), funding acquisition (equal), resources (equal), supervision (equal). Pei Zhang: conceptualization (lead), funding acquisition (lead), methodology (lead), project administration (lead), resources (equal), supervision (lead), writing – review and editing (lead).

Ethics Statement

This study involves non‐human primates ( Rhinopithecus bieti ) and does not involve human subjects. This research was approved by the Ethical Committee of Northwest University, Xi'an, China (Resolution No.: AWC‐20210402M), in compliance with the institutional guidelines for the care and use of non‐human primates.

Consent

The authors have nothing to report.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Appendix S1: ece373070‐sup‐0001‐AppendixS1.zip.