Assessing genetic admixture in fertility-associated genes of captive deer populations and its implications for reproductive management

List of abbreviations

BMP15	Bone Morphogenetic Protein 15
FOXL2	Forkhead Box L2
GDF9	Growth Differentiation Factor 9
MTNR1A	Melatonin Receptor 1A
MHCDQA1	Major Histocompatibility Complex Class II DQ Alpha 1
MHC	Major Histocompatibility Complex
DNA	Deoxyribonucleic Acid
EDTA	Ethylenediaminetetraacetic Acid
PCR	Polymerase Chain Reaction
dNTPs	Deoxynucleotide Triphosphates
MgCl₂	Magnesium Chloride
NCBI	National Center for Biotechnology Information
SNP	Single Nucleotide Polymorphism
PLINK	Name of a Genetics Software (not an acronym)
ANOVA	Analysis of Variance
Fst	Fixation Index
r	Pearson’s correlation coefficient
p	p-value (probability value)
UPM	Universiti Putra Malaysia
DVS	Department of Veterinary Services
AUP	Animal Use Protocol
CV	Cross-validation

1. Introduction

Wildlife farming constitutes an increasingly important component of global animal production systems, with deer farming occupying a prominent role due to its contributions to meat production, velvet antler harvesting, ecotourism, and conservation breeding initiatives [1]. The management of captive deer populations is increasingly relevant to the food industry in many parts of the world. Despite this growing interest, it is faced with challenges with respect to genetic control, reproductive efficiency, disease management, and behavioural compatibility. One very important issue that is often underestimated in multi-species captive systems is unintended hybridisation. This arises when closely related species or subspecies are maintained in proximity without strict breeding controls [2]. When such hybridisation occurs, it greatly affects population viability, particularly through its effects on fertility. The family Cervidae is made up of a diverse assemblage of species with distinct evolutionary histories. This is found in their ecological niche and reproductive strategies. Animals in this category include the spotted deer (Axis axis), the sambar deer (R. unicolor), and the Javan rusa (R. timorensis). Each of these species exhibits species-specific adaptations shaped by long-term evolutionary divergence [3].

When these species are managed together in captivity, opportunities for interspecific or intersubspecific mating increase, potentially leading to the breakdown of co-adapted genetic systems that underpin normal reproductive function [4]. Increasing evidence from captive and semi-managed cervid populations indicates that such hybridisation events are not uncommon and may contribute to reduced reproductive performance and long-term population instability [5].

Reproductive success in mammals is controlled by complex genetic networks involving genes that regulate ovarian development, gametogenesis, endocrine signaling, seasonal breeding responses, and immune–reproductive interactions. Several genes have been identified as particularly important in cervid and other mammalian reproductive biology. Bone morphogenetic protein 15 (BMP15) and GDF9 are oocyte-derived growth factors that regulate follicular development and ovulation. Disruptions in these factors have been shown to result in subfertility or infertility [6, 7]. FOXL2, on the other hand, is a transcription factor essential for ovarian differentiation and maintenance.  Ovarian failure as result of mutation in the FOXL 2 gene has been reported [8]. In addition, MTNR1A gene mediates seasonal breeding in cervids through its action on photoperiodism [9].

Seasonal breeding in cervids is strongly influenced by photoperiod, mediated largely through melatonin signaling pathways involving the MTNR1A gene [9]. In addition, genes within the major histocompatibility complex (MHC), such as MHCDQA1, influence reproductive success indirectly through immune competence and mate selection processes [10].

Hybridisation has the potential to interfere with these reproductive systems through multiple genetic mechanisms. One important outcome is outbreeding depression, which occurs when hybridisation disrupts co-adapted gene complexes that have evolved within species-specific genomic contexts [11]. Although hybrid offspring may sometimes exhibit heterosis for growth or survival traits, reproductive traits are particularly susceptible to negative epistatic interactions and genetic incompatibilities [12]. Such incompatibilities may arise through Dobzhansky-Muller interactions, in which alleles that function normally within their original genomic backgrounds interact adversely when combined in hybrids [13]. Hybridisation may also compromise species-specific adaptations related to reproductive timing, particularly in seasonally breeding species such as deer, where precise photoperiodic responses are critical for synchronising reproduction with environmental conditions [14]. Furthermore, altered epistatic relationships among fertility genes may disrupt signalling pathways essential for normal reproductive function [15].

In captive deer populations, the reproductive consequences of hybridisation may manifest as reduced conception rates, irregular oestrous cycles, increased embryonic loss, or complete reproductive failure [4]. These outcomes have significant economic implications for commercial deer farming operations and pose serious challenges for conservation breeding programmes aimed at maintaining genetically viable populations. Despite these risks, empirical studies directly linking genetic admixture within fertility-associated genes to reproductive performance in cervids remain limited. Existing research has often focused on genome-wide admixture without explicit attention to reproductive loci, or on individual genes without consideration of broader population structure [16].

Population genetic approaches offer powerful tools for examining the extent and consequences of hybridisation in managed populations. Model-based admixture analyses enable the estimation of individual ancestry proportions and the identification of genetic structure arising from historical or recent gene flow [17]. When applied to genes with direct functional relevance to reproduction, such analyses provide an opportunity to assess whether hybridisation has disproportionately affected genomic regions critical for fertility [18]. The ADMIXTURE algorithm estimates ancestry under a maximum-likelihood framework and allows the number of ancestral populations to be inferred through cross-validation, facilitating robust interpretation of population structure [17]. This approach has been successfully applied to investigate fitness-related traits in a range of wildlife species [19].

The present study addresses these gaps by analysing genetic admixture across five key fecundity-associated genes (BMP15, FOXL2, GDF9, MTNR1A, and MHCDQA1) in a captive deer population comprising Axis axis, R. unicolor, and R. timorensis. The population has experienced declining fertility, with interbreeding suspected as a contributing factor. By integrating sequence-based genetic data with admixture analysis and available reproductive records, this study aimed to (i) characterise the genetic structure and extent of admixture within and among the three species, and (ii) determine whether admixture patterns in fertility-associated genes are associated with reproductive performance. The findings are intended to inform evidence-based genetic management strategies for improving reproductive outcomes in captive cervid populations and to contribute to broader understanding of hybridisation effects on fitness-related traits in managed wildlife systems.

2. Materials and methods

2.1.  Study population and sample collection

The study was conducted on a captive deer population maintained at a commercial farming facility experiencing areported fertility decline. A total of 45 animals representing three cervid species were sampled: Axis axis (n = 11), R. timorensis (n = 19), and R. unicolor (n = 15). All individuals were permanently identified using ear tags, and available demographic information-including sex, age, and breeding history was recorded at the time of sampling. For female deer (n = 31), reproductive performance data, including conception rate, gestation length, and offspring survival, were obtained from farm breeding records covering the preceding three years.

Approximately 5 mL of whole blood was collected aseptically from the jugular vein of each animal into Ethylenediaminetetraacetic Acid (EDTA)-treated vacutainer tubes. Samples were transported to the laboratory under refrigerated conditions (4 °C) and processed within 24 hours of collection. All procedures involving animal handling and sample collection were approved by the Universiti Putra Malaysia Institutional Animal Care and Use Committee (Animal Use Protocol AUP-R047/2017) and were carried out in accordance with national guidelines for animal research.

2.2. DNA extraction and gene selection

Genomic DNA was extracted from whole blood samples using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. A NanoDrop 2000 device (Thermo Fisher Scientific, Waltham, MA, USA) was used to measure the concentration and purity of DNA spectrophotometrically. For further analysis, only samples with A260/A280 ratios between 1.8 and 2.0 were kept. Agarose gel electrophoresis was used to confirm DNA integrity.

Loci were identified based on their documented involvement in ovarian function, follicular development, seasonal breeding regulation, and immune–reproductive interactions in mammals. Five genes, the BMP15 gene, the FOXL2 gene, and the GDF9 gene, the MTNR1A gene, and the MHCDQA1 gene were selected for analysis.

2.3. Polymerase chain reaction and sequencing

Gene-specific primers targeting conserved coding regions of the five fecundity genes were designed using cervid reference sequences available in the GenBank database. Polymerase chain reaction (PCR) amplification was performed in 25 μL reaction volumes containing 50–100 ng of genomic DNA, 0.2 μM of each primer, 0.2 mM dNTPs, 1.5 mM MgCl₂, 1× PCR buffer, and 1 U of Taq DNA polymerase (Invitrogen, Carlsbad, CA, USA). PCR cycling conditions consisted of an initial denaturation at 95 °C for 5 minutes, followed by 35 cycles of denaturation at 95 °C for 30 seconds, primer-specific annealing temperatures for 30 seconds (Table 1), and extension at 72 °C for 1 minute. A final extension step was performed at 72 °C for 10 minutes.

Amplification products were visualised on 1.5% agarose gels stained with ethidium bromide. The QIAquick PCR Purification Kit (Qiagen) was used to purify the successful PCR products, and the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) was used for bidirectional sequencing. An ABI 3730xl DNA Analyser was used to perform the sequencing. Geneious Prime 2023 was used to build consensus sequences after raw chromatograms were manually examined (Auckland, New Zealand: Biomatters Ltd.).

Table 1. Selected primers for the amplification of protein coding sequences of BMP15, FOXL2, GDF9, MHCDQA1 and MTNR1A genes with their properties from Axis axis, R. timorensis and R. unicolor deer species reared at PTH Lenggong, Perak, Malaysia.

Genes	Primer sequence	Amplified region	Product size (bp)	Tm	GC content
BMP15		Exon 2	808
F	AGAGCCACTGTGGTTTACCG			59.97	55%
R	CTGCCGTTCGACGATTTCAC			59.91	55%
FOXL2		Exon1
F	TCCGGCTTCCTCAACAACTC		749	59.97	55%
R	GAAAGGGCAAAAAGGGCAGG			59.96	55%
GDF9		Exon1	777
F	TGATTCTCTGCCTTCTAGGGG			60.2	50%
R	AGGGAATGCCACCTGTGAAA			59.52	50%
MHC-DQA1		Exon 2	246
F	CTTGCTCCTCACTCCGACTCA			61.49	57.14%
R	GCTATGTTTCTCAGTGCACCC			59.26	52.38%
MTNR1A			856
F	GCCTGGCAGTTGCAGACCTG			60.5	55%
R	CATTTTTAAACGGAGTCCACC			60	45%

GC = Guanine and cytosine content, Tm = melting temperature of the primer.

2.4. Sequence processing and alignment

Sequence data were organised using a custom Python script (parse_fasta.py) to extract gene identity, species designation, and sample identifiers from sequence headers. Headers were standardised to the format >Gene|Species|SampleID, and individual FASTA files were generated for each gene. A corresponding metadata file linking samples to species was created to facilitate population genetic analyses.

Multiple sequence alignment for each gene was performed independently using MUSCLE v3.8.31 [20] with default parameters. Resulting alignments were inspected visually in Jalview v2.11.2 [21] to identify and correct potential alignment inconsistencies. Final alignments were exported in FASTA format for subsequent SNP detection.

2.5. SNP detection and quality filtering

A bespoke Python script (snp_to_plink.py) was used to extract variable sites from aligned sequences and transform them into formats that were compatible with population genetic software in order to identify single nucleotide polymorphisms. Polymorphic locations were defined as those with at least two nucleotide variations.

SNPs with more than 20% missing data across samples were excluded. In addition, loci with a minor allele frequency below 0.05 were removed to minimise the influence of rare variants and potential sequencing artefacts. Linkage disequilibrium pruning was conducted using PLINK with an r² threshold of 0.5, a sliding window of 50 SNPs, and a step size of 5 SNPs to reduce redundancy among closely linked markers.

Individuals with more than 20% missing genotype data were excluded from further analyses. After filtering, a total of 137 SNPs across the five genes were retained: BMP15 (26 SNPs), FOXL2 (32 SNPs), GDF9 (29 SNPs), MTNR1A (31 SNPs), and MHCDQA1 (19 SNPs).

2.6. Functional classification of SNPs into synonymous and non-synonymous mutations

Aligned sequences were color coded and read with “Color Align Conservation” on “Sequence Manipulation Suite (SMS) tools” to scan the entire length of the five gene segments for functional classification of the SNPs [22].

2.7. In silico protein-level analysis

A combination of two resources, EMBOSS Transeq [23] and ExPASy Translate [24], were used to translate all sequences into proteins. The translation was done using the following parameters: 1) Format of output: M, -, no spaces, compact 2) The forward strand of DNA. 3) Standard genetic code. 4) Trim: Sure, if necessary. Following translation, the proteins were aligned using MEGA X v10 and Sequence Manipulation Software (SMS) [22]. Following alignment, the patterns of variation were noted and utilised to categorise the proteins into groups of variants using "Protein Pattern Find" [25].

where p_i represents the proportion of ancestry attributed to the i^th ancestral component.

Differences in admixture metrics among species were evaluated using one-way analysis of variance (ANOVA). Associations between admixture measures and reproductive parameters (conception rate, gestation length, and offspring survival) were assessed using Pearson’s correlation coefficient. Partial correlation analyses controlling for age and parity were also performed. Conception rate was used because it was the most consistently recorded reproductive parameter in the farm records, while gestation length and offspring survival were also analyzed when available.

Partial correlation analysis was also performed; age and parity were controlled to account for potential confounding factors. Statistical analyses were conducted in R v4.2.2 [26], using the 'stats' package. fixation index (Fst) estimates were used to quantify the genetic differentiation among species using the hierfstat package in R [27].

2.10. Data visualization

Ancestry proportions were visualised using custom Python scripts (plot_admixture.py) based on the matplotlib library [28]. Individual ancestry profiles were displayed as stacked bar plots and were grouped by species and ordered to highlight admixture patterns. Cross-validation error values across K values were plotted to illustrate model selection. Scatter plots depicting relationships between admixture complexity and reproductive traits were generated with linear regression overlays. All figures were prepared at 300 dpi resolution in accordance with journal submission guidelines.

3. Results

3.1. Sequence data and SNP identification and analysis

High-quality sequence data were successfully generated for all the five fecundity-associated genes following alignment using MUSCLE. A total of 171 variable nucleotide positions were detected across the combined dataset. The number of polymorphic sites differed among loci, with FOXL2 and MTNR1A having the highest levels of sequence variability, while MHCDQA1 showed comparatively fewer polymorphisms. Quality control filters were applied to exclude SNPs with more than 20% missing data and minor allele frequencies below 0.05. A total of 137 SNPs were retained for population genetic analysis. The retained dataset provided sufficient resolution for admixture inference while maintaining data reliability.

3.2. Functional classification of SNPs into synonymous and non-synonymous mutations

All SNPs identified within the coding regions of the five fecundity genes (BMP15, FOXL2, GDF9, MHCDQA1 and MTNR1A) were classified according to whether they alter the encoded amino acid. The results show distinct patterns among genes (Table 2).

Table 2. Functional Classification of SNPs in the protein coding regions of some key fertility genes from three deer species reared at PTH Lenggong, Perak, Malaysia.

Gene	Total SNPs	Synonymous	Non-synonymous
BMP15	23	17 (73%)	6 (27%)
FOXL2	56	46 (83%)	10 (17%)
GDF9	14	11 (83%)	3 (17%)
MHCDQA1	20	5 (25%)	15 (75%)
MTNR1A	24	5 (19%)	19 (81%)
Total	137	88 (64%)	49 (36%)

These results demonstrate that while BMP15, FOXL2, and GDF9 are dominated by synonymous substitutions, MHCDQA1 and MTNR1A show a higher proportion of non-synonymous mutations, suggesting stronger potential functional effects in these loci.

3.3. Protein-level analysis and functional implications

To evaluate the functional consequences of non-synonymous SNPs, the nucleotide sequences were translated into amino acid sequences and compared with the wild-type proteins. This analysis allowed the identification of variant proteins and amino acid substitutions across the studied deer population.

BMP15: Three protein variants were identified with amino acid substitutions at positions 70, 71 and 125, involving residues E, R, Q, G and H.

FOXL2: Four variants were detected with substitutions affecting seven amino acid positions (25, 30, 50, 60, 100, 120 and 180) involving residues K, P, S, Q, G, A, N, Y, E and H.

GDF9: Two variants were observed with a single amino acid substitution at position 113 (Val and Ile).

MHCDQA1: Three variants were identified with amino acid substitutions involving I, D, V, T and N at positions 2, 3, 58 and 64.

MTNR1A: This gene showed the greatest degree of protein variation. Four variants were detected with thirteen amino acid substitutions across twelve positions (3, 16, 33, 128, 138, 206, 209, 210, 213, 214, 215 and 218).

Hardy–Weinberg equilibrium (HWE) tests conducted within each of the three recognised species indicated that most SNPs conformed to equilibrium expectations. However, a subset of loci exhibited significant deviations (p < 0.05), particularly within R. timorensis and R. unicolor populations (Table 3). These deviations are consistent with the presence of population substructure and recent admixture rather than genotyping artefacts, and therefore SNPs were retained for subsequent analyses.

3.4. Determination of optimal ancestral populations

Model-based admixture analysis was conducted across a range of ancestral population numbers (K = 1–5), with model fit evaluated using 10-fold cross-validation. Cross-validation error declined from K = 1 to K = 4, with a minimum value of CV error of 1.124 at K = 4 and overfitting at k = 5 (CV error = 1.303) as shown in Figure 1. These results support K = 4 as the optimal number of ancestral genetic clusters that contribute to the observed variation. By exceeding the number of sampled species, this suggests a more complex genetic structure than simple species-level differentiation.

3.5. Admixture patterns among deer breeds

Ancestry proportions inferred under the K = 4 model revealed extensive and heterogeneous admixture patterns across the sampled individuals as illustrated in Figure 2. To enhance interpretability, individuals were ordered according to their dominant ancestry component and relative contribution within each species.

3.5.1. Axis axis

On average, Ancestry 2 contributes about 49-51%, Ancestry 4 contributes about 32- 35%, Ancestry 3 contributes about 13-15% while Ancestry 1 contributes less than 1% as shown in Figure 2. This result indicates that Axis axis individuals in the population may not be genetically homogeneous and a likely historical introgression or shared ancestral polymorphism exists.

3.5.2. Rusa timorensis

Members of this group exhibited substantial heterogeneity in their ancestry composition. Some members of the group showed near-complete assignment to a single ancestral component (>99%), while others showed highly admixed profiles involving multiple components. Several members showed large contributions of about 42% to nearly 90% from Ancestry 1. Others, however, exhibited combinations of Ancestry 3 and Ancestry 4 (Figure 2). This diversity suggests the presence of multiple genetic lineages and ongoing or recent gene flow within the members of the R. timorensis group.

3.5.3. Rusa unicolor

Like the R. timorensis, members of this group also demonstrated pronounced admixture, though with slightly lower variability compared to R. timorensis. Some individuals showed near-exclusive assignment to Ancestry 4 (>99%), while others exhibited mixed ancestry profiles which involved Ancestry 1, Ancestry 2, and Ancestry 4 (Figure 2). The presence of shared ancestry components between R. unicolor and R. timorensis supports extensive historical or contemporary gene flow between the two species in captive environments.

3.6. Gene-specific admixture patterns

Separate admixture analyses conducted for each gene revealed notable locus-specific differences in genetic structure (Table 4). Optimal K values varied among genes, ranging from K = 2 for MHCDQA1 to K = 4 for FOXL2 and MTNR1A. These differences indicate that individual fertility-associated genes have experienced distinct evolutionary and demographic histories.

Table 4. Optimal K values and differentiation metrics for BMP15, FOXL2, GDF9, MTNR1A and MHCDQA1 genes from Axis axis, R. timorensis and R. unicolor deer species reared at PTH Lenggong, Perak, Malaysia.

Gene	Optimal K	FST between species	Predominant ancestry components by breed
BMP15	3	0.21	A. axis: 2, 3; R. timorensis: 1, 2, 3; R. unicolor: 1, 3
FOXL2	4	0.24	A. axis: 2, 4; R. timorensis: 1, 2, 3, 4; R. unicolor: 1, 4
GDF9	3	0.18	A. axis: 2, 3; R. timorensis: 1, 2, 3; R. unicolor: 1, 3
MTNR1A	4	0.31	A. axis: 2; R. timorensis: 1, 3; R. unicolor: 4
MHCDQA1	2	0.09	A. axis: 1, 2; R. timorensis: 1, 2; R. unicolor: 1, 2

The MTNR1A gene exhibited the highest level of interspecific differentiation (F_st = 0.31), with ancestry components closely aligned with species boundaries. In contrast, MHCDQA1 showed minimal differentiation (F_st = 0.09), with extensive sharing of ancestral components across species. BMP15, FOXL2, and GDF9 displayed intermediate differentiation levels (F_st = 0.18–0.24), which suggests partial retention of species-specific variation alongside substantial admixture.

We determined the average number of ancestral components contributing at least 10% to individual ancestry profiles in order to measure the degree of mixing in each breed. For Axis axis, R. timorensis, and R. unicolor, the mean (± SD) number of components was 2.00 (± 0.00), 1.69 (± 0.75), and 1.58 (± 0.67), respectively. These differences indicate a tendency toward more mixing in the Axis axis relative to the other species, although they were not statistically significant (one-way ANOVA, F (2,23) = 0.37, p = 0.69). As a gauge of admixture complexity, we also computed the average Shannon entropy of ancestry proportions. The mean (±SD) entropy values were 0.94 (±0.00) for Axis axis, 0.38 (±0.36) for R. timorensis, and 0.33 (±0.33) for R. unicolor. These differences approached statistical significance (one-way ANOVA, F (2,23) = 2.89, p = 0.076), suggesting greater admixture complexity in Axis axis compared to the other species.

3.8. Association between genetic admixture and reproductive performance

Correlation analyses were performed using reproductive records from 18 females with complete data (Table 5). A significant negative relationship was detected between admixture complexity (Shannon entropy) and conception rate (Pearson’s r = −0.58, p = 0.012). This indicates reduced fertility in individuals with more complex ancestry profiles (Figure 3).

Table 5. Reproductive parameters of deer hinds from Axis axis, R. timorensis and R. unicolor deer species reared at PTH Lenggong, Perak, Malaysia.

Animal ID	Age (years)	Parity (n)	Conception Rate (%)	Gestation Length (days)	Offspring Survival (%)
F1	5	3	80	240	100
F2	6	4	75	242	90
F3	4	2	90	238	100
F4	7	5	60	245	80
F5	5	3	85	241	100
F6	8	6	50	248	70
F7	4	2	95	237	100
F8	6	4	70	243	85
F9	5	3	80	240	100
F10	7	5	65	244	75
F11	4	2	90	239	100
F12	8	6	55	247	70
F13	6	4	75	242	90
F14	5	3	80	241	100
F15	7	5	60	246	80
F16	4	2	90	238	100
F17	8	6	50	249	65
F18	6	4	70	243	80

No significant associations were observed between admixture metrics and gestation length or offspring survival (p > 0.05). Analysis of individual ancestry components revealed contrasting effects: Ancestry 2 was positively correlated with conception rate (r = 0.61, p = 0.007), whereas Ancestry 1 exhibited a negative association (r = −0.47, p = 0.049).

4. Discussion

This study investigated genetic admixture within five fecundity-associated genes in a captive deer population which comprise Axis axis, R. timorensis, and R. unicolor. The aim was to understand potential genetic contributors to an observed fertility decline. The results showed extensive admixture within and between the three species, and notably, this is more pronounced among R. timorensis and R. unicolor, thereby providing some evidence that link admixture complexity to reduced conception rates. Furthermore, the findings highlight the risks associated with unmanaged interbreeding in captive deer populations and underscore the importance of integrating genetic data into reproductive management strategies. At the molecular level, protein-level analysis in silico confirms that several SNPs identified in the coding regions lead to amino acid substitutions that produce multiple protein variants across the three deer species. In particular, MTNR1A and MHCDQA1 contain the highest proportions of non-synonymous mutations, suggesting stronger selective pressure or functional diversification in these genes, while BMP15, FOXL2 and GDF9, although containing mostly synonymous substitutions, still exhibit specific amino acid changes that generate distinct protein variants. Collectively, these findings support the hypothesis that genetic admixture among the deer populations may contribute to functional diversity in key fertility-related proteins, potentially influencing reproductive physiology and adaptation [29, 30].

Building on this, the main finding of this study is the identification of four ancestral genetic clusters (K = 4) which contribute to the observed variation across the sampled individuals, in spite of the fact that the population is made up only three recognised species. This discrepancy therefore indicates that genetic structure within the population cannot be fully explained by species designations alone, and indeed, such patterns have been reported in cases involving historical divergence among lineages, undocumented founder effects, or introgression from unsampled or extinct source populations [2]. Moreover, these processes are often made worse by management practices that facilitate close contact among species with incomplete reproductive isolation, in captive systems. Consistent with this, the pronounced heterogeneity observed within R. timorensis and R. unicolor, ranging from near-pure ancestry profiles to highly admixed individuals, strongly suggests ongoing or recent gene flow between the species. This interpretation is further supported by the substantial sharing of ancestral components and aligns with previous reports that documented hybridisation among cervids in captive or semi-managed systems [5]. Similarly, the admixed genetic profile observed in Axis axis further suggests either historical introgression or the persistence of ancestral polymorphisms predating species divergence. Consequently, the presence of this complex structure makes breeding management very difficult and challenges assumptions of genetic purity based solely on phenotypic identification.

At a finer scale, at the gene-level, analyses showed marked differences in differentiation and admixture patterns across loci, and this reflects the distinct evolutionary pressures which act on genes with different biological functions. Among the genes examined here, MTNR1A exhibited the highest interspecific differentiation (Fst = 0.31) [4], with the ancestry components here closely aligned with species boundaries. This pattern is therefore consistent with strong divergent selection acting on photoperiodic regulation of reproduction, which is a trait that is particularly critical in seasonally breeding cervids [9]. As a result, disruption of such finely tuned regulatory mechanisms through hybridisation may lead to mismatches between reproductive timing and environmental cues, potentially contributing to reduced fertility [13]. In contrast, differentiation across species in MHCDQA1 was minimal (Fst = 0.09) and there was extensive sharing of ancestral components, a pattern which is characteristic of MHC genes that are often subject to balancing selection that maintains polymorphism across species boundaries [10]. While increased MHC diversity can enhance immune competence, admixture at these loci may also interfere with species-specific mate choice mechanisms or immune–reproductive interactions, with unpredictable consequences for reproductive success [10]. Meanwhile, the intermediate differentiation observed for BMP15, FOXL2, and GDF9 suggests partial retention of species-specific variation alongside substantial admixture, and given that these genes play central roles in ovarian development, folliculogenesis, and ovulation [6–8], their coordinated function depends on finely balanced regulatory and epistatic interactions. Therefore, hybridisation may disrupt these interactions by introducing alleles that have not co-evolved within the same genomic background, potentially impairing reproductive efficiency even in the absence of overt pathological mutations [14].

Importantly, one of the most significant findings of this study is the significant negative association between admixture complexity, measured using Shannon entropy, and conception rate, and notably, this relationship remained robust after controlling for age and parity, indicating that genetic admixture itself contributes independently to reduced fertility. Such a pattern is therefore consistent with outbreeding depression, a phenomenon in which hybridisation disrupts co-adapted gene complexes and favourable epistatic interactions, leading to reduced fitness [11]. Although hybridisation can sometimes produce heterosis for growth or survival traits, reproductive traits are often more sensitive to genetic incompatibilities [12]. In line with this, the contrasting effects observed for individual ancestry components further support this interpretation, as the positive association between Ancestry 2 and conception rate suggests that this genetic background may harbour allele combinations better suited to the captive environment or more compatible with the broader genomic context of the population. Conversely, the negative association between Ancestry 1 and fertility, particularly given its strong representation in Rusa hybrids, points to the accumulation of genetic incompatibilities consistent with Dobzhansky–Muller interactions [13]. Even though the present analysis focused on a limited number of genes, their direct involvement in reproductive function strengthens the inference that admixture within these loci reflects broader genomic disruptions affecting fertility, and thus, the observed patterns likely represent the cumulative effects of multiple interacting genes rather than the influence of any single locus in isolation. Nevertheless, it is important to recognise that reproductive performance in cervids is influenced by a complex interaction of genetic and environmental factors [31, 32], and accordingly, in addition to female genetic background, conception success may be affected by male fertility, nutritional status, environmental conditions, disease burden, and farm management practices [33-36]. Consequently, the association observed between admixture complexity and conception rate should be interpreted as evidence of a genetic correlation rather than direct causation [37, 38], and therefore, further studies incorporating controlled breeding data would be valuable to isolate genetic effects more precisely [39, 40].

From a practical standpoint, the findings of this study have immediate relevance for the management of captive deer populations and broader implications for conservation breeding programmes, as the clear link between admixture complexity and reduced conception rates highlights the risks associated with maintaining multiple closely related species in shared or poorly segregated breeding environments. Without such control, hybridisation may progressively erode reproductive performance, undermining both economic viability and conservation objectives. To address this, genetically informed management strategies are essential, and these may include routine genetic screening to assess individual ancestry profiles, structured mating systems designed to limit further admixture, and the strategic use of individuals with lower admixture complexity in breeding programmes. Furthermore, where the introduction of external stock is considered, careful genetic evaluation is required to avoid introducing additional incompatibilities or disrupting locally adapted gene combinations [4]. Beyond commercial farming, these results serve as a cautionary example for conservation initiatives, since captive breeding programmes intended to preserve genetic diversity may inadvertently compromise fitness if hybridisation is not carefully managed, and therefore, integrating molecular genetic monitoring into conservation protocols is critical to maintaining both genetic integrity and reproductive viability in managed wildlife populations.

Finally, several limitations should be acknowledged, as the sample size, while sufficient to detect meaningful patterns, limits the generalisability of the findings to other populations. In addition, the analysis was restricted to five candidate genes, providing a focused but incomplete view of genome-wide admixture dynamics, while the absence of confirmed purebred reference populations from outside the study farm also constrains definitive interpretation of ancestral components. In light of these limitations, future studies should prioritise genome-wide approaches using high-density SNP arrays or whole-genome sequencing to capture a more comprehensive picture of admixture and identify additional loci associated with reproductive performance. Moreover, functional studies examining gene expression or hormonal regulation in individuals with contrasting admixture profiles would further clarify the mechanistic basis of fertility decline, and longitudinal monitoring across generations would also be valuable for assessing the long-term consequences of admixture and the effectiveness of genetic management interventions.

5. Conclusion

This study demonstrates that genetic admixture within key fecundity-associated genes is significantly associated with reduced reproductive performance in a captive multi-species deer population. The identification of four ancestral genetic components across three recognised species reveals extensive historical and contemporary gene flow, particularly between R. timorensis and R. unicolor. Importantly, increased admixture complexity, quantified through Shannon entropy, was strongly and independently linked to lower conception rates, providing empirical support for outbreeding depression affecting reproductive function.

The gene-specific analyses further indicate that hybridisation does not affect all loci equally. The strong differentiation observed at MTNR1A suggests disruption of species-specific photoperiodic regulation, while the extensive sharing of MHCDQA1 alleles reflects balancing selection typical of immune-related genes. Intermediate differentiation in BMP15, FOXL2, and GDF9 genes central to ovarian physiology-suggests that admixture may interfere with coordinated reproductive pathways, thereby contributing to fertility decline.

These findings have direct implications for the genetic management of captive cervid populations. Uncontrolled interbreeding in mixed-species settings can erode reproductive efficiency and compromise long-term population viability. Integrating genetic screening, structured mating systems, and ancestry-informed breeding strategies offers a practical pathway for mitigating fertility loss while preserving genetic diversity. Beyond the study population, this work underscores the broader importance of genetic monitoring in wildlife farming and conservation breeding programs, where maintaining both genetic integrity and reproductive fitness is essential for sustainable population management.

Acknowledgments

We thank the Department of Veterinary Services (DVS) Malaysia for providing all the animals for data generation and support staff for sampling.

Authors’ contributions

MS Yahaya: Conceptualization, Investigation, Methodology, Data analysis, Writing - Original Draft. MS Salisi: Methodology, Investigation, Conceptualization. Resources, Writing - Review & Editing. NM MD Isa: Methodology, Resources, investigation. YM Goh: Methodology, Resources. Writing - Review & Editing. Abd Wahid Haron: Conceptualization, Methodology, Resources, Writing - Review & Editing.

Conflict of interest statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding

This study was supported by the Universiti Putra Malaysia Excellence grant and the Department of Veterinary Services (DVS) Malaysia with grant number 9617200.

Ethics approval

AUP-R047/2017.