Species composition and environmental adaptation of indigenous Chinese cattle

Indigenous Chinese cattle combine taurine and indicine origins and occupy a broad range of different environments. By 50 K SNP genotyping we found a discontinuous distribution of taurine and indicine cattle ancestries with extremes of less than 10% indicine cattle in the north and more than 90% in the far south and southwest China. Model-based clustering and f4-statistics indicate introgression of both banteng and gayal into southern Chinese cattle while the sporadic yak influence in cattle in or near Tibetan area validate earlier findings of mitochondrial DNA analysis. Geographic patterns of taurine and indicine mitochondrial and Y-chromosomal DNA diversity largely agree with the autosomal cline. The geographic distribution of the genomic admixture of different bovine species is proposed to be the combined effect of prehistoric immigrations, gene flow, major rivers acting as genetic barriers, local breeding objectives and environmental adaptation. Whole-genome scan for genetic differentiation and association analyses with both environmental and morphological covariables are remarkably consistent with previous studies and identify a number of genes implicated in adaptation, which include TNFRSF19, RFX4, SP4 and several coat color genes. We propose indigenous Chinese cattle as a unique and informative resource for gene-level studies of climate adaptation in mammals.

studies using mtDNA 7,8,16,17 and Y-linked markers 18 . A genetic diversity study using microsatellite markers clustered Chinese indigenous cattle breeds into one taurine and four indicine groups 19 .
In addition to taurine and indicine cattle, several other bovine species have been living in southern China and Southeast Asia, including banteng (Bos javanicus), gaur (Bos gaurus) or gayal (Bos frontalis), which may have been the dominant cattle species until 4500 BP 4,12 . Gayal in Yunnan province of China carried indicine or taurine mtDNA but gaur Y chromosome, indicating its hybrid origin 20 . Meanwhile, in Tibetan Autonomous Region (TAR) of China, bidirectional introgression between yak and cattle has also been reported [21][22][23][24] . Genetic admixture has been identified between zebu and Bali cattle (domestic banteng) in Indonesia 22,25 . A previous study on hair color and blood protein polymorphism provided evidence of banteng introgression into Hainan cattle in southeastern China 26 , which was confirmed by genomic SNP array data 25 .
Genomic SNP array has become a powerful tool for population genomics studies in animals. A recent genomic variation study revealed a worldwide pattern of genetic admixture in domestic cattle 25 . Other studies focused on Creole 27 , American 28 , East African zebu 29 and Korean cattle 30 . These advanced approaches also allow the genomic localization of genes involved in the adaptation to natural or artificial selective constraints 27,[31][32][33][34][35][36][37][38] . In the current study, we generated 50 K SNP genotypes to infer the fine-scale characterization of unique species composition of highly diverse Chinese cattle. In addition, we performed a whole genome-scan for adaptive differentiation and association analyses with environmental and morphological population-specific covariables to detect genes that responded to adaptive constraints.

Results
Genomic variation. Observed heterozygosity ( Table 1) ranged from 0.145 to 0.327 in Chinese cattle populations ( Fig. 1). Mongolian cattle (MG, NM) and Kazakh cattle had the highest values but southern and southwestern Chinese cattle populations were the lowest. This is most likely explained by the ascertainment bias, by which the heterozygosity of indicine cattle is underestimated. Indeed, the observed heterozygosity correlates negatively (r 2 = 0.96) with the zebu ancestry. A similar trend has previously been observed in West-African cattle 39 . Bali cattle (0.026), gayal (0.059) and yak (0.029) also have relatively low levels of heterozygosity as normally observed with SNP panels designed for a different species 40 . Population structure. Five methods were implemented in this part to explore the population structure.
PCA. Figure 2a shows a scatter plot of the first two principal components (PCs), allowing to assess the structuring of genetic diversity across all the 45 sampled populations, including outgroup species (Bali cattle, yak, gaur and gayal). The first PC accounting for 17.1% of total variation separates taurine and indicine cattle as well as other bovine species. The predominant taurine breeds from northern China and Tibet Autonomous Region of China (TAR) cluster with European breeds while populations from southeastern and far-southwestern China are closer to Indian zebu. The second component accounting for 3.1% of all variation displays the contrast of Chinese cattle to other bovine species (Bali, gayal and yak) and also differentiates European cattle breeds. In an analysis of only Asian cattle (Fig. 2b), the first PC that explains 14.5% of the genetic variation again corresponds to the taurine-indicine separation, but the second PC that explains 1.45% of the genetic variation represents a gradient from India via Southeast Asia to Southeastern China.
Model-based clustering. Figure 3a shows an unsupervised hierarchical clustering performed on the whole data set (i.e., including outgroup species) with different values of K, the number of clusters. K = 2 reproduces the first PCA coordinate with taurine and indicine clusters and a taurine-indicine gradient from north to south. Higher values of K differentiate European from Asian taurine breeds (K = 3), zebu from other bovine species (K = 4), and southern Chinese from Asian indicine cattle (K = 5). Increasing K further generates separate clusters for European cattle breeds whereas Chinese cattle populations tend to show admixed ancestries. The cross-validation error gave the lowest value at K = 17 ( Supplementary Fig. S1), which differentiates Bali cattle, gayal, yak, the European taurine breeds and indicine cattle from India/Pakistan, Southeast Asia and southeastern China, respectively. Figure 3b shows a supervised clustering with prior population information for taurine breeds (SH, HOL, BSW, SIM), indicine breeds (GIR, SAHW) and other bovine species. In this analysis, Bali cattle represents the banteng and the domestic gayal replaces the wild gaur since the latter was too inbred for this analysis. This analysis reproduces the zebu-banteng hybridization in Indonesia 41 and yak introgression into Tibetan cattle, but also suggests minor banteng and gayal into southeastern China ( Supplementary Fig. S2). Supplementary Figure S3 shows the geographic distribution of the inferred species components.
NeighborNet network. In a NeighborNet graph constructed from the matrix of Reynolds' distances between populations using Splitstree (Fig. 4), European and Indian cattle are at the extreme ends of the network, which is entirely in agreement with both the first PCA coordinate and the K = 2 clustering. Figure 4 also reproduces the intermediate positions of the predominantly taurine-indicine breeds from TAR or central China and also of the predominantly indicine breeds from central or southwestern China. In addition, the network confirms the affinity of Indonesian and southeastern Asian continental breeds with Bali cattle and gayal.
Population mixture. The taurine-indicine mixed composition of Chinese cattle was confirmed by the sensitive f4 ratio test (Supplementary Table S1). The estimated taurine ancestry ranges from close to 0.0 in southeastern China to close to 1 25,41 . Although Bali cattle are relatively closely related to gaur, replacing Bali cattle by gaur as source of admixture generates only a moderately negative value for the Indonesian breeds (Fig. 5b). This is even observed for Bali cattle as a test breed and may reflect the inbreeding of the gaur samples. However, the same plot shows relatively low (i.e., negative, indicating gene flow) values for the southern Chinese breeds. In combination with the supervised model-based clustering (Fig. 3b), this may suggest that these breeds have been introgressed by gaur and/or gayal in addition to banteng, the wild ancestor of Bali cattle. A f4(SH,X;GAY,YAK) plot ( Supplementary Fig. S4A) generated negative values for the indicine GIR, SAHW and BD as test breeds. Since this is not observed with the wild gaur instead of the domestic gayal ( Supplementary  Fig. S4B), this indicates indicine introgression into the domestic gayal population. This is confirmed by two other observations: • The statistic f4(GAU,GAY;X,YAK) is negative for all test breeds, but clearly more negative for indicine than for taurine breeds ( Supplementary Fig. S4C). The same patterns are also observed if GIR is replaced by SH ( Supplementary Fig. S4A,B), which is phylogenetically close to the indicine GIR. Apparently the allele sharing of GAY with GIR or SH outweighs the allele sharing expected because of the phylogenetic relationship of gayal and Bali cattle, which most likely is an effect of the ascertainment bias of the SNP panel towards taurine breeds 42 .
Uniparental markers. Supplementary Tables S2 and S5 show mtDNA and Y-chromosomal haplotype distributions of Chinese cattle breeds, based on our data supplemented by literature data (Supplementary Table S3). These uniparental markers show a north-to-south taurine-indicine gradient that resembles closely the autosomal cline ( Fig. 6). Haplotype diversity (Supplementary Table S2 and Fig. S5C) shows that the diversity of taurine mtDNA hardly decreased from north to south, while indicine mtDNA clearly decreases from southeastern and southwestern China to northern China.
Detection of genomic regions subjected to adaptive constraints. For a genome-wide scan for adaptive differentiation, the XtX differentiation statistics was estimated for each SNP under the so-called core model and visualized in a Manhattan plot (Fig. 7a). Among 84 significant SNPs, SNP Hapmap28985-BTA-73836 on BTA5 was the most significant. The three synthetic environmental covariables are associated with 185 significant SNPs (Supplementary  Table S4), from which SNP Hapmap44345-BTA-119580 on BTA11 was the most significant (Fig. 7b). From 30  SNPs associated with morphological covariables (Supplementary Table S4), SNP ARS-BFGL-NGS-67505 on BTA25 was most significant (Fig. 7c).
By applying the sliding window approach, 27 significant genomic regions were obtained (in Table 2), of which eight showed significantly differentiated SNPs; 22, 0 and 4 displayed SNPs associated with the first, the second and the third environmental co-variable, respectively; and 4 displayed SNPs associated with the morphological covariable. Note that only 3 regions (out of the 8) containing significantly differentiated SNPs did not contain any SNPs associated with the population-specific covariables studied. Finally, a total of 28 candidate genes were annotated in the significant regions using UCSC (https://genome.ucsc.edu/) ( Table 2).

Discussion
We have investigated the species composition and genomic, mitochondrial as well as Y-chromosomal variation in Chinese cattle populations. A clear north-south gradient of taurine and indicine cattle ancestries combined with banteng, gayal and yak introgressions into southern and southwestern Chinese cattle populations defines the pattern of admixture among Chinese indigenous cattle, which are supposed to underlie local breeding objectives and   adaptation to different agro-ecological environments. Genome scan for adaptive differentiation and association with population-specific covariable identify regions and candidate genes relevant for environmental adaption and morphological differentiation in Chinese cattle.
Previous genetic diversity studies using mtDNA and Y-linked markers have characterized a north-south gradient of taurine and indicine admixture in Chinese cattle, which is consistent with the transition from humpless  to humped morphology 7,8,[16][17][18] . A microsatellite study differentiated five groups of Chinese indigenous cattle breeds 19 .
Our PCA (Fig. 2) and model-based clustering (Fig. 3, Supplementary Fig. S3) patterns as well as the NeighborNet graph (Fig. 4) reveal a clear transition from taurine cattle in the north to zebu in the south with consistent admixture levels within the breeds and a clear demarcation: • Cattle from Manchuria, Inner Mongolia and northwestern China with 4 to 8% indicine admixture. This group corresponds to the taurine type 1 19 (this type also comprises Tibetan cattle) and to Group 9B in the Felius classification 43 . • Taurindicine cattle in TAR and northern China above the Yellow or Wei River with 28 to 35% indicine admixture, denoted as zebu type 2 19 and belonging to the Huanghuai group of central Chinese cattle (Group 10A) 43 . • Taurindicine cattle below the Yellow River with 61 to 73% indicine genome, corresponding to indicine types 1 and 4 19 . The northernmost LX, JX and NY breeds belong to the Huanghai group (Group 10A) 40 , but the other breeds to the Changzu group (Group 10B) 43 . • Predominant indicine cattle in southern China with an indicine genomic component of 79 to 87% (zebu type 3 19 , Group 10B 43 ). • Indicine cattle with a >90% indicine genome (also indicine type 3 19 , Group 10B 43 ).
Cattle from Northeast Asia and northern China have typical taurine morphological features. Its genetic distinctiveness and old origin evolved from the unique mtDNA haplogroup T4 found in modern breeds from Japan, Korea, Mongolia 6,44 , Siberia (Yakut cattle) 9 , northern China 7,8 ( Supplementary Fig. S5A) and in ancient cattle in northern China from 2300-4500 BP 5 . In addition, a separate position of Japanese and Korean cattle was found in previous study 25 . In this study, we found that northern Chinese, South Korean and Japanese cattle share genetic ancestry with the Siberian Yakut and the indigenous northeastern China close to Korea. These breeds represent the eastern range of Turano-Mongolian cattle, which have retained their original dark-brown coat color pattern in Mongolian and Korean cattle.
The admixture pattern, however, suggests European cattle influence to northern Chinese cattle diversity. Possibly migrations of nomads in the steppes of Central Asia and Mongolia, which in the Middle Ages led to the establishment of the Mongolian empire, facilitated eastward as well as westward gene flow across Eurasia. Additionally, in the past few decades programs have been implemented in China to upgrade productivity by crossing local breeds with European breeds 3,12 . Influence of European cattle was captured by model-based clustering analysis, e.g. Brown Swiss and Simmental in Kazakh, Holstein in LS, and Limousine in LX, JN and WL (Fig. 3). From north to south, levels of taurine autosomal, mitochondrial and Y-chromosomal DNA decrease, but are still appreciable in southwestern China ( Supplementary Figs S3 and S5). The occurrence of indicine mtDNA in mixed taurine-indicine cattle is a unique feature of Chinese breeds 11 . The high taurine mtDNA diversity in southern China (Supplementary Fig. S5C) indicates an absence of a major founder effect. This is more compatible with immigration before the arrival of indicine cattle around 3000 BP than with a later introgression into an existing indicine population. The immigration of zebus from the present Myanmar and Indochina to the north with a plausible contribution of an eastward gene flow in western China 15 resulted in significant indicine components in northern Chinese and Mongolian cattle ( Supplementary Figs S3 and S5). The pattern of indicine mtDNA diversity ( Supplementary Fig. S5C) suggests a population bottleneck when they crossed the Pearl River in southeastern China.
Y-chromosomal and autosomal indicine components correlate well (Fig. 6). Exceptions are the Guanlin (GL) and Nanyang (NY) breeds with relatively low indicine Y-chromosomal frequencies. Remarkably, the indicine autosomal component has a discontinuous distribution (Figs 2-4; Supplementary Fig. S3) with the largest gap across the Yellow River separating JN (29% zebu) from ES (61%) and NY (64%). This river might be a physical barrier of gene flow, but this cannot explain the absence in our panel of cattle with an indicine component between 35% and 61%. It might be hypothesized that in cattle with equal taurine and indicine components, outbreeding depression outweighs heterosis, for instance by incompatibility of genes from different origins conferring fitness 45 . It is again remarkable that cattle at both sides of the river are categorized as belonging to the Huanghuai group and resemble western Asian and African cattle because of their similar cervico-thoracic hump 2 .
The NeighborNet shows that indicine cattle from far southwestern China (DH, BN), which are found in the west of the Mekong River, are more closely related to Thai cattle than to southeastern Chinese cattle. This is also supported by the PCA and admixture patterns (Figs 2 and 3). The separate position of Thai cattle confirms the results of previous study 14 . The Mekong River acted also as a genetic barrier for swamp buffalo 46 .  Another potential source of diversity in southern Chinese cattle are the introgressions from other bovine species living in China and Southeast Asia, including yak, gayal and its wild ancestor gaur 43 , and banteng, represented by its domestic relative Bali cattle. (Fig. 3; Supplementary Fig. S2). There are several examples of hybridization of different bovine species. Yak mtDNA has been detected in indigenous cattle distributed on the Qinghai-Tibetan plateau and in Diqing cattle (DQ) of Yunnan province 21,23 . We did not detect yak mtDNA in our cattle panel, but we detected influence of yak in the Tibetan LZ population based on the model-based clustering analysis.
A study of blood protein polymorphism 26 suggested banteng ancestry in Hainan cattle. Using mtDNA and microsatellite genotyping, Mohamad et al. 47 characterized banteng admixture in Indonesian zebu breeds. This was confirmed by 50 K SNP analysis 25 , which also detected a low level of banteng introgression in southern Chinese breeds. However, analyzing Chinese cattle together with both banteng (represented by its domestic derivative Bali cattle) and gaur (the wild ancestor of gayal) with model based clustering (Fig. 3) and f4-statistics (Fig. 5) provided consistent evidence of both gayal and banteng introgression into WL, WN and HN breeds from southeastern China and also into Thai zebu at a relatively low proportion.
Conversely, similar f4-statistics ( Supplementary Fig. S4) suggested introgression of zebu into gayal, which has been confirmed in the gayal population from Yunnan carrying both indicine and taurine mitochondrial genomes 20 . Similarly, cattle introgression has been detected in yak populations 24 .
The genetic variation described above reflects the combined effect of prehistoric immigrations of taurine and indicine cattle, subsequent gene flow between populations, local selection objectives and environmental adaptation. Indigenous Chinese cattle with indicine-taurine ratios varying between zero and one and subject to a broad range of climates is a valuable resource to identify potential genomic regions and functional genes underlying the environmental adaption. By combining signals of population differentiation (XtX) and association with three synthetic environmental covariables and one synthetic morphological covariable ( Supplementary Fig. S6), we identified 27 genomic regions and 28 candidate genes targeted by natural or artificial selection (Table 2).
Interestingly, 12 out of the 27 regions overlap with core selective sweep (CSS) regions 48 , while 20 and 23 regions overlap with breed-wise and breed group-wise hotspots of selective sweeps, respectively 49 (Supplementary  Table S5). However, we report for the first time the region BTA12: 34.0-35.5 Mb, which harbors TNFRSF19. This member of the TNF-receptor superfamily 50 is highly expressed during embryonic development 51 .
In other studies, strong signals of selection in tropical cattle have been detected on BTA5 32,39,52,53 . Notably, Porto-Neto et al. 32 identified a 20 Mb region on BTA5 with effects on parasite resistance, yearling weight, body condition score, coat color and penile sheath score. We found a significant signature of selection for XtX and the environmental PC1 and PC3 as well as the morphological PC1 on BTA5 region 69.5-71.0 (Table 2), which contains the candidate gene RFX4. This gene is a member of Regulatory Factor X (RFX) family of transcriptional regulators that influence MHC class II expression 54 and play a critical role in brain development 19,55 . It was also found to affect heifer fertility in tropical composition breed Brangus 56 .
Coat color is an important target of selection in many domestic animals. The common denotation of yellow cattle for the indigenous Chinese cattle refers to its predominant light to dark brown color. In current study, selection signatures were identified near several known color genes, including KITLG (near SNP BTA-74300-no-rs on BTA5) 57,58 , and LEF1 32,59 (here indicated by a peak in the XtX GWAS on BTA6). These genes and another candidate gene MCM6 (near ARS-BFGL-NGS-92772 on BTA2, also identified by Hudson et al. 53 ) overlap with pigmentation QTL regions underlying UV-protection 60 . The environmental PC1 signal near IGFBP7 and the combined XtX-morphological signal near ADRA1D (Table 2) are close to the coat color genes KIT 60,61 and ATRN 60,62 , respectively.
We further detected an environmental PC1 association signal near SP4 (Sp4 transcription factor) as novel candidate gene on BTA4. Interestingly, this signal is overlapped with a selection signature region in African cattle 34 . SP4 is a member of the Sp1-family of zinc finger transcription factors and is required for normal murine growth, viability, and male fertility 63 . In cattle, SP4 was suggested to have effect on body size and testicular growth from birth to yearling age 64 . In human, diseases associated with SP4 include bipolar disorder 65 and schizophrenia 66 .
It is interesting to note that Chinese and African cattle have developed independently a variable taurine-indicine ancestry following a gradient from tropical to temperate climates. An attractive opportunity is a detailed comparison of gene variants involved in climate adaptation by using whole genome sequence data 34 . It may be anticipated that in both regions adaptation to agro-ecological constraints is mediated by recruiting and combining gene variants from taurine and indicine origins with possible original contributions in Chinese indigenous cattle from the indicine mtDNA, and the minor gayal, and banteng and yak genomic ancestry.

Methods
Ethics statement. The protocols for collection of the blood and hair samples of experimental individuals were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) at China Agricultural University. All experiments were performed in accordance with approved relevant guidelines and regulations.
Samples collection and genotyping. We collected samples of 437 animals from 24 populations ( Table 1), twenty of which are indigenous cattle populations from northeastern China, central China, southeastern China, southwestern China, far southwestern China and or TAR (Fig. 1). We also examined Bangladeshi cattle and German Simmental, and two related bovine species, the gayal and yak. Samples were genotyped with Illumina BovineSNP50 BeadChip using standard procedures 67 . Genotypes are accessible via the WIDDE repository (http://widde.toulouse.inra.fr/widde/).
We compared these newly generated data with published genotypes of some European and Asian cattle breeds, Bali cattle and gaur ( cattle, zebu, and three species related to cattle (Bos javanicus -banteng, Bos gruniens -yak and, Bos frontalis -gayal). Gaur (Bos gaurus) 69 was used instead of gayal in f4 analysis (see below) of the species composition of Chinese cattle because of the indicine zebu introgression into gayal.
Population genetic analysis. We used PLINK 68 to calculate the observed homozygosity for each population. Three complementary methods were used to analyze the genetic diversity among populations. First, a Principal Component Analysis (PCA) was carried out to investigate the pattern of genetic differentiation among populations and individuals using the R package SNPRelate 70 , which performs eigen-decomposition of the genetic covariance matrix to compute the eigenvalues and eigenvectors. Second, population structure was evaluated by unsupervised and supervised model-based hierarchical clustering implemented in the Admixture software 71 . The results were visualized using the program Distruct 72 . Third, a NeighbourNet network was constructed using Reynold's distances between populations using Splitstree 4.13 73 .
To investigate species composition of Chinese cattle, we used the four-population test (f4-statistics) implemented in ADMIXTOOLS 74 . Additionally, taurine and indicine ancestries in Chinese cattle populations were quantified via the f4 ratio estimation in ADMIXTOOLS, which allows inference of the admixture proportions without access to accurate surrogates for the ancestral populations 74  Genome-scan for adaptive differentiation and association with environmental and morphological covariables. Whole genome-scans for adaptive differentiation and association with population-specific co-variables were performed with BayPass 2.1 57 . The underlying models explicitly account for the covariance structure among the population allele frequencies, which make the approach particularly robust to complex demographic histories 57 . Identification of overly differentiated SNPs was based on the XtX statistics 57,78 estimated under the core model of BayPass. To calibrate the XtX's, a pseudo-observed data set (POD) containing 250,000 SNPs simulated under the inference model with hyperparameters equal to those estimated on the real data set was generated and further analyzed under the same conditions following the procedure described in Gautier 57 . In particular, we ensured that the posterior estimate of the scaled covariance matrix of population allele frequencies (Omega) obtained with the POD was similar to that obtained on the real data since the FMD distance between the two matrices was found equal to 0.28 57 . Similarly, the posterior means of the two hyperparameters a and b for the Beta distribution of across population allele frequencies obtained on the POD (a = b = 1.02) were almost equal to the ones obtained in the original data (a = b = 1.00). Taken together, these sanity checks indicated that the POD faithfully mimics the real data set, allowing us to define a 0.1% significance threshold on the XtX statistics (XtX = 32.3) to identify genomic regions harboring footprints of selection.
We collected values for six environmental covariables, i.e. average temperature, average relative humidity, sunshine, average air pressure, wind speed, and precipitation from China Meteorological Administration (http:// data.cma.cn/). Values for 10 morphological covariables, i.e. male body weight, female body weight, male height, female height, male body length, female body length, male heart girth, female heart girth, male fore-shank circumference, and female fore-shank circumference were provided for 17 out of the 20 studied indigenous breeds in National Commission of Animal Genetics Resources 3 , but these data were not available for the populations LZ, BN and HH. We carried out a PCA on scaled variables for environmental and morphological co-variables separately. The first three environmental PCs and the first morphological PC were retained as uncorrelated co-variables for association studies (Supplementary Tables S6 and S7).
Genome-wide analysis of association with population specific co-variables was carried out using the default options of the AUX model parameterized with the scaled covariance matrix (Omega) obtained on the real data set as described above. This model allows to account explicitly for multiple testing issue by integrating over (and estimating) the unknown proportion of SNPs actually associated with a given covariable. The support for association of each SNP with each co-variable was evaluated by computing Bayes Factor (BF) and a BF > 15 is considered as decisive evidence for association 57 .
As a matter of expedience, we applied a sliding window approach to identify the main genomic regions of interest as described in previous study 57 . Briefly, the UMD3.1 bovine genome assembly 79 on which all the SNPs were mapped was first split into 4718 consecutive 1-Mb windows (with a 500-kb overlap). For each window, we counted the number n s of SNPs that were either significantly differentiated at the 0.1% threshold (i.e., with XtX > 32.3) or associated (BF > 15) with at least one of the four population-specific covariables. Windows with n s > = 2 were deemed significant and overlapping windows were further merged.