Whole-genome sequencing of Zhoushan cattle and Wenling cattle populations
We collected seven individuals of Zhoushan cattle (Fig. 1a, upper panel). We also collected nine individuals of Wenling cattle (Fig. 1a, lower panel). Wenling cattle have a prominent hump on the back, dewlap, and larger ears, suggesting that its genetic background is largely B. indicus (Fig. 1a, lower panel). We performed whole-genome sequencing of these samples. To resolve their phylogenetic positions and interrelationships within domesticated cattle, we combined our data of 16 cattle individuals with publicly-available whole-genome sequencing data of five individuals from the Angus breed, a typical B. taurus in Europe, and 33 individuals from nine breeds with genetic backgrounds similar to B. indicus3, giving a total of 54 individuals (Fig. 1b, c; Table S1). We performed read trimming and aligned the trimmed reads to the UOA_Brahman_1 assembly of the cattle genome11. This assembly represents the maternal haplotype of an F1 hybrid of Brahman cattle (dam) and Angus (sire)11. After variant calling and filtering, we identified 32,970,327 single-nucleotide polymorphisms (SNPs) and 3,331,322 small indels. Based on this genomic variant information, we conducted the population genomic analyses.
Phylogenetic analysis of Zhoushan cattle and other cattle breeds. (a) Gross appearance of Zhoushan (upper panel) and Wenling cattle (lower panel). Note that Zhoushan cattle have a dark black coat colour. The arrow indicates the curving horn of Zhoushan cattle. (b) Geographic map indicating the origins of Zhoushan (green dot) and Wenling (orange dot) cattle analysed in this study. We also examined other Chinese cattle (red dots) whose genome sequencing data were available. (c) Regional map around the Zhoushan islands. Wenling, Wannan, and Guangfeng are mainland regions close to the Zhoushan islands. (d) Neighbour-joining tree of the 54 domesticated cattle. The scale bar represents pairwise distances between different individuals. The maps were constructed by R38 and R packages of maps v3.3.0 (https://cran.r-project.org/web/packages/maps) and mapdata v2.3.0 (https://cran.r-project.org/web/packages/mapdata).
Genetic relationship between Zhoushan cattle and other domesticated cattle
To reveal the phylogenetic positions and interrelationships of Zhoushan and other domesticated cattle, we performed population genomic analyses on 54 cattle individuals. First, we calculated the pairwise evolutionary distance between individuals and generated a neighbour-joining (NJ) tree to reconstruct the phylogenetic relationships between individuals of Zhoushan and other domesticated cattle (Fig. 1d). In the NJ tree, cattle clustered consistently with their geographical location (Fig. 1d). Angus individuals formed a sister group to all other individuals, including Zhoushan cattle, Wenling cattle, and other B. indicus (Fig. 1d). The individuals of Zhoushan and Wenling cattle formed monophyletic groups and were sisters to each other (Fig. 1d). The cattle in Guangfeng formed another monophyletic group and were sisters to both Zhoushan and Wenling cattle (Fig. 1d). Cattle in Wannan, Ji’an, and Leiqiong formed a single group, sister to the cattle of Zhoushan, Wenling, and Guangfeng (Fig. 1d). Zhoushan, Wenling, Guangfeng, Wannan, and Ji’an are geographically close to each other (Fig. 1b, c). The cattle of Dianzhong and Wenshan, which are in the south part of China, were distant from them (Fig. 1d). Cattle in Pakistan and India were located near the root of the phylogenetic tree (Fig. 1d). The branch lengths of Zhoushan cattle were shorter than other B. indicus cattle, suggesting the reduced genetic diversity of Zhoushan cattle (Fig. 1d).
To estimate the relatedness between Zhoushan and other domesticated cattle, we performed unsupervised clustering analysis with ADMIXTURE v1.3.0 software (https://dalexander.github.io/admixture/index.html)12. At K = 2, Angus cattle were distinct from all other cattle (Fig. 2a). At K = 3, Zhoushan and Wenling cattle were newly segregated from other cattle, suggesting that these two cattle breeds are genetically close to each other (Fig. 2a). The cattle of Guangfeng, Wannan, Ji’an, Leiqiong, and Wenshan had intermediate genetic structures between Zhoushan cattle and Dianzhong cattle (Fig. 2a). At K = 4, Zhoushan cattle and Wenling cattle were separated from each other (Fig. 2a).
Admixture and principal component analysis of Zhoushan cattle and other cattle breeds. (a) Admixture plot (K = 2, 3, 4) for the 54 cattle individuals. Each individual is shown as a vertical bar divided into K colours. (b) PCA plot showing the genetic structure of the 54 cattle individuals. The degree of explained variance is given in parentheses. Colours reflect the geographic regions of sampling in Fig. 1d. The cluster composed of cattle in Wenling, Guangfeng, Wannan, Ji’an, and Leiqiong is highlighted in the black dotted ellipse. (c) Estimate of the effective population sizes of Zhoushan (green) and Wenling (orange) cattle over the past 100 generations.
To infer the population structure of cattle individuals analysed in this study, we conducted principal component analysis (PCA). The top three principal components accounted for 21.1% of the total variance (Fig. 2b). In the first component of PCA, Angus individuals were separated from all other cattle (Fig. 2b). Additionally, cattle of Wenling, Guangfeng, Wannan, Ji’an, and Leiqiong formed a cluster (dotted ellipse in Fig. 2b). In the second component of PCA, individuals of Zhoushan cattle were separated from all other cattle (Fig. 2b). In the third principal component, Wenling cattle individuals were separated from all other cattle (Fig. 2b).
We estimated the trends of the effective population size of Zhoushan and Wenling cattle over the past 100 generations (Fig. 2c). Both populations showed decreasing trends of effective population sizes (Fig. 2c). The effective population size of Zhoushan cattle was estimated to be smaller than that of Wenling cattle, suggesting the effect of island isolation on the genetic diversity of Zhoushan cattle (Fig. 2c).
Detection of candidate genes associated with dark black coat colour of Zhoushan cattle
To identify putative genes associated with the dark black coat colour of Zhoushan cattle, we searched genomic regions where the same mutations were shared between Zhoushan cattle and Angus cattle. To achieve this, we calculated the average fixation index (Fst) values in 40 kb windows with 10 kb steps (Fig. 3a). We identified four peaks of Fst at chromosomes 2, 4, 8, and 18 (Fig. 3a). Among these peaks, the highest peak of Fst was identified in the region from 51.05 to 51.35 Mbp on chromosome 18 (Fig. 3a, b). This region contains 18 genes (Fig. 3c). We searched for genes that have mutations altering the amino acid sequence and have been reported to be involved in the regulation of coat colour. Among these 18 genes, only the gene of melanocyte-stimulating hormone receptor (MC1R) is known to involved in the regulation of coat colour13,14,15. Therefore, we regarded MC1R as a strong candidate gene associated with the dark black coat colour of Zhoushan and Angus cattle (Fig. 3c). This gene is located in the region between 51,094,227 bp and 51,095,177 bp on chromosome 18. MC1R is expressed in the skin melanocyte and plays a crucial role in regulating animal coat colour formation16. Mutations of MC1R have been reported to be associated with black coat colour in some animals, such as cattle17, sheep16, pigs18, reindeer19, and geese20. In the protein-coding region of MC1R, we identified one missense mutation (c.583T > C, p.F195L) and one synonymous mutation (c.663C > T) (Figs. 3d, 4a). The missense mutation is located in the fifth transmembrane region of MC1R (Fig. 4b). All seven Zhoushan cattle were homozygous for the missense mutation (Figs. 3d, 4a). Four of five Angus individuals were homozygous for the missense mutation, and the remaining one was heterozygous for the missense mutation (Figs. 3d, 4a). Conversely, only 19% (8/42) and 33% (14/42) of B. indicus individuals were homozygous or heterozygous, respectively, for the missense mutation (Figs. 3d, 4a). The remaining 48% (20/42) of individuals of B. indicus were homozygous for the wild-type allele (Figs. 3d, 4a). We also found that the p.F195L mutation is also present in MC1R of Black Angus (accession number: ABX83563.1) in the NCBI Protein database (Fig. S1). Furthermore, we identified 15 upstream variants and three downstream variants in the intergenic regions between neighbouring genes (Table S2).
Genomic regions associated with dark black coat colour of Zhoushan cattle. (a) Manhattan plot for average Fst values in 40 kb windows with 10 kb steps between Zhoushan cattle plus Angus and other B. indicus. A region with an average Fst of more than 0.6 is coloured in green. The arrow indicates the highest peak. The x-axis represents chromosomal positions, and the y-axis represents the average Fst values. (b) Manhattan plot on chromosome 18 for average Fst values in 40 kb windows with 10 kb steps between Zhoushan cattle, Angus, and other B. indicus. (c) Regional plot around the MC1R gene. The genotype of each individual at each variant site is shown. The genotype homozygous for the reference allele is coloured grey. Heterozygous variants are coloured blue. The homozygous genotype for alternative alleles is coloured light blue. Note that homozygous genotypes for alternative alleles are enriched in Zhoushan and Angus cattle in this region. (d) Regional plot showing the mutations around MC1R gene.
Secondary structure of MC1R and protein sequence alignment of MC1R orthologs. (a) Regional highlight of the c.583 T > C mutation of MC1R. The genomic region from 51,094,590 to 51,094,598 bp on chromosome 18 is shown. Note that MC1R is located on the reverse strand. (b) Secondary structure of MC1R. MC1R is a seven-transmembrane receptor. The p.F195L mutation is located in the 5th transmembrane region and enclosed by the red circle. This figure is generated by using the Protter server application39. (c) Multiple sequence alignment of MC1R orthologs. The black rectangle highlights the 195th phenylalanine residues. The red rectangle encloses the p.F195L mutation in Zhoushan cattle. The cladogram of the species is shown to the left of the species name. The cladogram topology is derived from a previous study40.
To characterise the missense mutation of MC1R (c.583T > C, p.F195L) found in Zhoushan and Angus cattle, we estimated the degree of evolutionary conservation of the 195th phenylalanine of MC1R. We obtained various MC1R orthologs of vertebrates from eight eutherian mammals, two marsupial mammals, four reptiles, two birds, two amphibians, one lobe-finned fish, one polypterus fish, four teleost fish, and two cartilaginous fish (Table S3). We aligned these 26 sequences with MC1R of Zhoushan cattle and B. indicus (Fig. 4c). This analysis revealed that the 195th phenylalanine of MC1R is highly conserved among vertebrates (Fig. 4c).
Furthermore, we verified whether any larger structural variants are spanning the MC1R region (chr18:51,058,185–51,148,307 bp) of Zhoushan cattle and Angus. If there are large structural variants in this region for these breeds, we should see regions where the read depth distributions are different among the groups. We assessed the integrated read depth distributions of Wenling cattle (n = 9), Zhoushan cattle (n = 7) and Angus (n = 5) (Fig. 5a). The read depth distribution was very similar among the three groups suggesting that there are not large structural variants spanning the MC1R region in these breeds (Fig. 5a). We also collected the sequence reads mapped to this region, and performed BreakDancer to detect structural variants21. However, no structural variants were detected in this region in any breeds. Moreover, we compared the reference genome sequence in MC1R region of the UOA_Brahman_1 assembly and that of the UOA_Angus_1 assembly11. The UOA_Brahman_1 assembly represents the maternal haplotype of an F1 hybrid of Brahman cattle (dam) and Angus (sire), and the UOA_Angus_1 assembly represents its paternal haplotype11. The results showed that the genome sequence in the MC1R region are highly preserved between these two assemblies (Fig. 5b).
Read depth distribution, genome alignment and admixture analysis of the MC1R region. (a) Read depth distributions in the MC1R region. The left panel shows the read depth distributions in the region from 51,058,185 to 51,148,307 bp on chromosome 18. The right panel shows the read depth distributions in the region from 51,090,618 to 51,099,796 bp on chromosome 18. For each breed, the sequencing reads were integrated. The first track represents read depth distribution in each breed, and the second track represents read alignments to the reference genome. For a given base position, if the base call in the sequencing read and the corresponding base in the reference genome are different, adenine is shown in green, thymine in red, guanine in orange, and cytosine in blue. (b) Dot plots showing the genome alignments of the MC1R regions of the UOA_Angus_1 assembly (chr18:49,477,288–49,566,766 bp) and the UOA_Brahman_1 assembly (chr18:51,058,185–51,148,307 bp). The left panel shows the genome alignment by minimap2 aligner and the right one shows the genome alignment by LASTZ aligner. The region corresponding to the MC1R gene body is highlighted in red. (c) Admixture analysis of the MC1R region. The SNPs located in the MC1R region (chr18:51,058,185–51,148,307 bp) were collected and subjected to admixture analysis. The order of the samples is the same as in Fig. 2a.
Finally, we deduced the origin of the MC1R haplotype in Zhoushan cattle. We collected the SNPs located in the MC1R region (chr18:51,058,185–51,148,307 bp) from all individuals and performed admixture analysis using these SNPs. The result showed that Zhoushan cattle and Angus shared highly similar genetic components (Fig. 5c). However, the other individuals of B. indicus showed genetic components that differed from both Zhoushan cattle and Angus (Fig. 5c). These results suggest that the MC1R haplotype in Zhoushan cattle is derived from B. taurus, even though the genome of Zhoushan cattle as a whole is that of B. indicus.
Source: Ecology - nature.com
