Whole-genome sequences of 89 Chinese sheep suggest role of RXFP2 in the development of unique horn phenotype as response to semi-feralization

Abstract

Background: Animal domestication has been extensively studied, but the process of feralization remains poorly understood.

Results: Here, we performed whole-genome sequencing of 99 sheep and identified a primary genetic divergence between 2 heterogeneous populations in the Tibetan Plateau, including 1 semi-feral lineage. Selective sweep and candidate gene analysis revealed local adaptations of these sheep associated with sensory perception, muscle strength, eating habit, mating process, and aggressive behavior. In particular, a horn-related gene, RXFP2, showed signs of rapid evolution specifically in the semi-feral breeds. A unique haplotype and repressed horn-related tissue expression of RXFP2 were correlated with higher horn length, as well as spiral and horizontally extended horn shape.

Conclusions: Semi-feralization has an extensive impact on diverse phenotypic traits of sheep. By acquiring features like those of their wild ancestors, semi-feral sheep were able to regain fitness while in frequent contact with wild surroundings and rare human interventions. This study provides a new insight into the evolution of domestic animals when human interventions are no longer dominant.

Figure 1:

Genetic relationships and population structure in Chinese sheep. a) Geographic distribution of the Chinese indigenous sheep breeds (PT, Prairie Tibetan; OL, Oula; VT, Valley Tibetan; BY, Bayinbuluke; WZ, Wuzhumuqin; T, Tan; CB, Cele Black; STH, Small-tailed Han; H, Hu) and a European-originated breed (AM, Australian Merino) sampled in the present study. The background color of the sheep pictures represents their lineages (red: TBS, Tibetan sheep; blue: MGS, Mongolian sheep; green: EUS, European sheep). b) Neighbor-joining tree of the 10 breeds based on F_ST distances. c) Principal component plot. The first (PC1) and second (PC2) principal components are shown. d) Population structure analysis of 99 sheep, where number of ancestral clusters were set from K = 2–4.

Figure 2:

Manhattan plot of genome-wide selective sweep signals (F_ST and log-scaled H_P ratio) in 4 sheep breeds. For each metric, a 30-kb sliding window with a step size of 15 kb was applied. F_ST distances were calculated between each of the 4 breeds (PT, OL, VT, or BY) vs MGS (WZ, T, STH, H, and CB). The log-scaled H_P ratio was calculated as −log₂(H_{P|PT, OL, VT or BY}/H_P|MGS), a positive value of which suggests reduction of variability in the breed.

Figure 3:

Candidate genes associated with selective sweeps in semi-feral sheep. a) A Venn plot showing numbers of overlapping candidate genes among 4 breeds (PT, OL, VT, and BY). b) Sweep signal metrics for genes selected from feralization-related categories as well as 3 genes associated with hypoxic adaptation. c) A summary of feralization-related adaptation observed in semi-feral sheep. Affected functional terms were manually summarized based on gene ontology enrichment analysis of the candidate genes, as well as literature mining. Numbers denote the count of candidate genes within each major category.

Figure 4:

Selective sweep over the horn-related gene RXFP2. a) Statistics plotted over a approximately 400-kb region surrounding RXFP2, including 1) population differentiation (F_ST) between PT, OL, VT, and BY vs MGS; 2) intrapopulation heterozygosity in PT, OL, VT, and BY, calculated as Z-transformed log₂(H_{P|PT, OL, VT or BY}/H_P|MGS); 3) haplotypic length measured by Z-transformed XP-EHH_{PT, OL, VT or BY vs. MGS}. b) Haplotypic distributions among 99 sheep of a local region of RXFP2 (chromosome 10: 29 400 000–29 550 000 bp). Biallelic SNPs are shown in blue and yellow. c) Alignment of the RXFP2 protein sequences from 9 vertebrate species. Two protein variants (RXFP2: 627 and 641) with top F_ST in PT and OL are indicated in red. For 627, PT and OL have the variant allele, whereas for 641, they have the reference allele. The dots in the alignment denote amino acids that are identical with those in PT and OL. d) Distribution of the haplotype frequency of 2 protein-altering variants (RXFP2: 627 and 641) in 1155 sheep. “Haplotype1” corresponds to V627 + E641 (OAR10_29 461 968: C + OAR10_29 462 010: C) and “Haplotype2” corresponds to M627 + K641 (OAR10_29 461 968: T + OAR10_29 462 010: T).

Figure 5:

RXFP2 haplotype is correlated with horn shape and size. a) Features of SHE-type and TCF-type horns. b) Association between 8 SNPs and horn phenotypes (size and shape) analyzed in 182 PT sheep. After testing all combinations of genetic models and confounding effects (Supplementary Fig. S12), an additive model (assume A as major allele, a as minor allele, we have code 2 for AA, 1 for Aa, and 0 for aa) was applied for horn size, and a recessive code (1 for AA, 0 for Aa and aa) was applied for horn shape; pair-wise LD between SNP pairs were plotted at the bottom, where numbers represent D΄ statistics. c) Box plot of individual horn sizes among different OAR_29 461 968 genotypes; P value was calculated by linear regression based on additive genetic model, and the fitting line is shown in red. d) Distribution of OAR10_29 461 968 genotypes among PT sheep with different horn shapes.

Figure 6:

Gene expression patterns of RXFP2. a) Expression of RXFP2 and β-actin in 13 tissue samples from PT sheep: 1, heart; 2, liver; 3, spleen; 4, lung; 5, kidney; 6, muscle; 7, brain; 8, ovary; 9, corpus uteri; 10, adipose; 11, thyroid; 12, soft horn; and 13, horn periosteum. b) Expression pattern of RXFP2 in SHE-type, TCF-type, scurred soft-horn tissues examined by RT-PCR (left) and real-time PCR (right); error bars denote SD of the mean; groups with significant differences (*, P < 0.05; **, P < 0.001) are indicated. c) Scatter plot on RXFP2 expression and horn size; the fitting line of linear regression is shown in blue. d) Western blot analysis of soft-horn tissues with different horn types, using antibodies of RXFP2 and β-actin.