Selection of Cashmere Fineness Functional Genes by Translatomics

Cashmere fineness is an important index to evaluate cashmere quality. Liaoning Cashmere Goat (LCG) has a large cashmere production and long cashmere fiber, but its fineness is not ideal. Therefore, it is important to find genes involved in cashmere fineness that can be used in future endeavors aiming to improve this phenotype. With the continuous advancement of research, the regulation of cashmere fineness has made new developments through high-throughput sequencing and genome-wide association analysis. It has been found that translatomics can identify genes associated with phenotypic traits. Through translatomic analysis, the skin tissue of LCG sample groups differing in cashmere fineness was sequenced by Ribo-seq. With these data, we identified 529 differentially expressed genes between the sample groups among the 27197 expressed genes. From these, 343 genes were upregulated in the fine LCG group in relation to the coarse LCG group, and 186 were downregulated in the same relationship. Through GO enrichment analysis and KEGG enrichment analysis of differential genes, the biological functions and pathways of differential genes can be found. In the GO enrichment analysis, 491 genes were significantly enriched, and the functional region was mainly in the extracellular region. In the KEGG enrichment analysis, the enrichment of the human papillomavirus infection pathway was seen the most. We found that the COL6A5 gene may affect cashmere fineness.


INTRODUCTION
Cashmere goat is a unique animal husbandry resource in China (Cai et al., 2020;Jin et al., 2020), which has made outstanding contributions to the needs of the livestock textile industry and human society (Yang et al., 2020). The characteristics of cashmere show significant differences due to different goat species and regions (Jin et al., 2017). Cashmere goat has a double coat, with the primary hair follicles producing coarse hair and the secondary hair follicles producing fine hair (Dai et al., 2019). The cashmere thickness produced by the secondary hair follicles varies among individuals. Cashmere fiber is famous for being slender and soft. To improve the economic value of cashmere, the fineness of its fibers needs to be reduced, making cultivating high-yield and high-quality cashmere goat varieties the core element to improve the economic value of cashmere fiber. Cashmere fineness is a quantitative characteristic, which is determined by micro effective genes. Molecular breeding research in domestic and foreign studies have proved that the single nucleotide polymorphisms (SNPs) of KAP (Jin et al., 2011), RPL (Mahata et al., 2012), FGF (Liu et al., 2009), KRT (Hui, 2020), PROP (Zeng et al., 2011), MAF70 (Arranz et al., 2001), and KIFI (Li, 2009) genes are correlated with cashmere fineness. Many studies have identified some genes related to goat wool fiber growth characteristics, such as DSG1, IGF-IR, KRTAPs, ILK, and KRTAP genes (Rufaut et al., 1999;Jia et al., 2006;Yu et al., 2011;Yang et al., 2012;Liu et al., 2015). Full transcriptomics also proved that NFKBIA, TCHH, COL1A1, CXCL8, and LTBP2 genes are closely related to cashmere fineness, and researchers have been looking for ways to improve the quantity and quality of cashmere (Bai et al., 2016;Zhang et al., 2018), as well as analyzing the key genes, signal pathways, and expression regulation level under different cashmere state diameters . Despite the efforts, all the genes regulating cashmere fineness from the translation level remain unknown. Due to the gap between the transcriptome and proteome in omics, only 20-40% of proteins in mammalian cells are determined at the transcription level (Tian et al., 2004;Cox et al., 2005). Current translatomics can make up the gap between the two. In the study of cashmere fineness, translatomics can be used to identify genes that are not translated into proteins at the transcriptional level. Translatomics provides new insights into gene expression (Courtes et al., 2013).
Translatomics is the sequencing and analysis of translated RNA molecules, which can accurately quantify the genes being translated, and compare the translation amount of different samples of genes under different physiological and pathological conditions or different treatments. At the same time, translation control is the key determinant of protein abundance, which in turn determines cell state (Sendoel et al., 2017). The gene expression process reveals that there are specific phenotypic characteristics among species, but their evolutionary process is uncertain outside the transcriptome. There are studies on coevolution at the level of mammalian transcriptome and translatome. Ribosome analysis and RNA sequencing are used to analyze the three organs of five mammals (human, macaque, mouse, opossum, and platypus) and birds (chicken), coevolutionary analysis (brain, liver, and testis) shows that translation regulation widely exists in different organs, especially in spermatogenic cell types of testis, and some genes evolve faster at the translatome level . Ribosomal analysis can also evaluate translation efficiency on a genome-wide scale, which has been previously proved in yeast (McManus et al., 2014), nematodes (Stadler and Fire, 2013), primate (Wang et al., 2018) cells, and hybrid mouse cells, it was also found that translation efficiency was a momentous predictor of protein level in mouse fibroblasts (Hou et al., 2015). These studies provide preliminary insights into the coevolution model of the transcriptome and translatome. Therefore, protein abundance seems to be mainly regulated by ribosomes, highlighting the importance of translation control (Gebauer and Hentze, 2004;Sonenberg and Hinnebusch, 2009). Using microRNA (miR-430) in zebrafish to investigate its translational repression and mRNA decay, we found that translation repression occurs before mRNA decay, which is induced by reducing the translation initiation rate, and that mRNA decay is induced by deadenylation. Besides, microRNA has been proposed to affect protein translation by reducing the rate of translation initiation (Bazzini et al., 2012). On the other hand, because the translatome studies RNA molecules being translated, which includes the RNA molecules that are traditionally considered noncoding, it can provide direct translation evidence for the study of these new translatable molecules (Gerashchenko et al., 2012). In addition, the lack of a strict direct correlation between gene and protein levels limits translation studies by combining the transcriptome and proteome. Considering the high cost associated with protein synthesis, the dominant role of translation regulation is meaningful. Therefore, it promotes the progress of translatomics technology (King et al., 1996;Kirkpatrick et al., 2005).
According to Tian Wenliang, the standard of cashmere fineness in China is: the diameter of coarse hairs: 16.0-18.5 μm and the diameter of fine hairs: 15.5-16.0 μm. In 2010, the cashmere of LCG was the thickest among all cashmere goats, and its diameter was 20.32 μm (Tian, 2015). It can be found that the main disadvantage of cashmere in LCG is that the cashmere is thicker. In this study, we found the key genes regulating cashmere fineness by using translatomics, and understood the regulatory relationship of related genes. We used Ribo-seq to test the fine skin samples and coarse skin samples of LCG, and found the key genes, differential genes, and co-expression genes related to cashmere fineness through GO function enrichment analysis and KEGG pathway enrichment analysis. These findings pave the way for the study of the regulation mechanism of cashmere fineness and the protection and cultivation of cashmere varieties.

Ethics Statement
The whole process of experiments was based on guidelines from the Animal Experimental Committee of Shenyang Agricultural University (Shenyang, China, 201906099).

Sample Preparation
The sample collection was a crucial step for Ribo-seq since it was the starting point of library construction.

Library Construction for Ribo-Seq
To digest RNA other than RPFs, cell or tissue lysate was treated with unspecific endoribonuclease RNase |. Isolation of monosomes was performed by size-exclusion chromatography with MicroSpin S-400HR columns. The RNA samples were then treated with an rRNA depletion kit to deplete the samples of as much rRNA contamination as possible before PAGE purification of the relatively short (20-38 nt) RPFs. Following PAGE purification, both ends of the RPF were phosphorylated and ligated with 5′ and 3′ adapters, respectively. Then the fragments were reversely transcribed to the cDNAs and amplified by PCR (Aeschimann et al., 2015). After library construction, the concentration of the library was measured by The Qubit ® 2.0 Fluorometer and adjusted to 1 ng/ uL. An Agilent 2100 Bioanalyzer was deployed to examine the insert size of the acquired library. At last, the accurate concentration of the cDNA library was again examined using qPCR. Once the insert size and concentration of the library became identical, the samples could then be subjected for sequencing.

Sequencing
After library preparation and pooling of different samples, the samples were subjected to Illumina sequencing. Commonly, the Ribo-seq uses PE150 (paired-end 150 nt) sequencing for 15 G raw data.

Quality Control for Raw Data
Firstly, the initial data (in the format of FASTQ) and the adapter were processed to delete the 3′ ends sequence and obtain the clean data of Q20. The following analysis was based on the clean data.

Mapping
Ribo-seq used TopHat2 for genome mapping. TopHat2 is an enhanced version of TopHat, using short read aligner Bowtie to align the RNA-seq reads to mammalian-sized genomes and analyzing the mapping result to identify splice junctions. TopHat2 allows variable-length indels with respect to the reference genome, which give it the ability to accurately align the transcriptomes in the presence of insertions, deletions, and gene fusions (Kim et al., 2013).

Quantification of Gene Expression Level
Quantification of mapped results to gene level was carried out using HTSeq. HTSeq is a Python package that calculates the number of mapped reads to each gene (Anders et al., 2015). RPKM values were generated to represent the gene expression level of each specific gene. RPKM is the abbreviation of "Reads Per kilobase of transcript, per Million mapped reads," which normalizes both sequencing depth and gene length (Gross et al., 2013).

Differential Expression Analysis
For samples with biological replicates, the DESeq2 R package (1.14.1) was used for differential expression analysis. DESeq2 provides statistical routines for determining differential expression in digital gene expression data using a model based on negative binomial distribution (Wang et al., 2010). The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with p < 0.05 found by DESeq2 were assigned as differentally expressed (Tang et al., 2007).

GO and KEGG Enrichment Analysis
Gene Ontology (GO) is a major bioinformatics initiative to unify the representation of gene and gene product attributes across all species. GO covers three domains: cellular component, molecular function, and biological process. KEGG (Kyoto Encyclopedia of Genes and Genomes) is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances (Kanehisa et al., 2008). In the KEGG pathway database, the wiring diagram database, is the core of the KEGG resource. It is a collection of pathway maps integrating many entities including genes, proteins, RNAs, chemical Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 12 | Article 775499 compounds, glycans, and chemical reactions, as well as disease genes and drug targets, which are stored as individual entries in the other databases of KEGG.

P-Site Analysis
Identifying the A-and P-site locations on ribosome-protected mRNA fragments from Ribo-Seq experiments was a fundamental step in the quantitative analysis of transcriptome-wide translation properties at the codon level. The P-site (for peptidyl) is the second binding site for tRNA in the ribosome. During protein translation, the P-site holds the tRNA which is linked to the growing polypeptide chain. When a stop codon is reached, the peptidyl-tRNA bond of the tRNA located in the P-site is cleaved releasing the newly synthesized protein. Since translation occurs at the A-and P-sites, the identification of these sites was critical to address translation-related questions . Novogene used the Ribocode package to analyze the P-site using Ribo-seq data (Xiao et al., 2018).

uORF Analysis
An upstream Open Reading Frame (uORF) is an Open Reading Frame (ORF) within the 5′ untranslated region (5′UTR) of an mRNA. uORFs can regulate eukaryotic gene expression. Associated with mRNA-seq, all identified ORFs by Ribo-seq were classified. Ribotape was then used to analyze the motif of translated/untranslated uORFs, which can be used to study the base composition bias of uORF sequences.

RESULT Quality Control of Sequencing Data of Six LCGs
Through high-throughput sequencing, the raw reads sequences of fine LCG samples accounted to 50666392, 50286756, and 52749309. Low quality data accounted for about 0.21%. The rest were clean reads. The percentage of sequencing sequences that could be located on the genome was about 88.86%. The average percentage of sequencing sequences with unique alignment positions on the reference sequence was about 46.28%. About 99.06% of the base group had a mass value greater than 20. The average proportion of filtered rRNA to total clean reads was about 76.51%. The average proportion of filtered tRNA in the total number of clean reads was about 4.33%.
The raw reads sequences of coarse LCG samples accounted to 49894812, 51599006, and 59901380. Low quality data accounted for about 0.21%. The rest were clean reads. The percentage of sequencing sequences that could be located on the genome was about 81.64%. The average percentage of sequencing sequences with unique alignment positions on the reference sequence was about 39.10%. About 99.03% of the base group had a mass value greater than 20. The average proportion of filtered rRNA to total clean reads was about 76.66%. The average proportion of filtered tRNA in the total number of clean reads was about 3.77% ( Table 1).
The effect of experimental enrichment can be evaluated by counting the length of ribosome protected RNA fragments  Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 12 | Article 775499 6 (RPFs). The length statistics of Ribo-seq clean reads of fine and coarse samples of LCG showed that when the enrichment length was 22 nt, the enrichment frequency was the highest, 21.66 and 14.38%, respectively ( Figures 1A,B). The genomic region sequencing distribution is shown in the figure: the coding region of fine Liaoning cashmere goat accounted for 21.40%, the UTR region accounted for 59.00%, the intron region accounted for 2.80%, and the intergenic region accounted for 16.79% (Figure 2A). The coding region of coarse Liaoning cashmere goat accounted for 17.66%, the UTR region accounted for 55.69%, the intron region accounted for 2.43%, and the intergenic region accounted for 24.23% ( Figure 2B).

Quantitative Analysis and Distribution of Gene Expression in Six LCGs
The translation level of a gene protein is directly reflected by the abundance of ribosome binding on its corresponding transcript. The higher the binding abundance is, the higher the level of gene translation is. We quantitatively analyzed the gene expression level of the top 30 genes in each sample, as shown in Table 2.
In addition to the true translation level, the reads count was positively correlated with the sequencing depth. Generally, the gene expression value is not expressed by reading count, but by RPKM, which corrects the sequencing depth and gene length successively. After calculating all gene expression values (RPKM) of each sample, we showed the distribution of gene expression levels of different samples by box graph. From the box graph, we can see that there were differences in the expression levels of all the genes detected among the six samples, and the box graph of sample 3 was smaller than that of other samples, so we can see that the differences between the genes detected in sample 3 and other samples were more obvious. At the same time, obvious different genes were found between the two groups. The coincidence of the two groups was the co-expressed genes, and the noncoincidence was the differentially expressed genes ( Figures 3A-C).  Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 12 | Article 775499 8 FIGURE 5 | Enrichment analysis in LCG. (A) GO enrichment analysis histogram, the abscissa is the GO term, and the ordinate is the significance level of GO term enrichment. Height was positively correlated with the significance. Different colors represent BP, CC, and MF. (B) Scatter plot of enrichment analysis, the horizontal coordinate in the Panel is the ratio of the number of different genes annotated to the GO term and the total number of different genes. The ordinate is the GO term, the size of points represents the number of genes annotated to the GO term, and the level of enrichment varies from purple to red. (C) KEGG enrichment analysis histogram, the horizontal coordinate in the panel is the KEGG pathway, and the longitudinal coordinate is the significance level of channel enrichment, and the height of the histogram indicates the degree of enrichment. (D) Scatter plot of KEGG enrichment analysis, the abscissa is the ratio of the number of differential genes annotated to the KEGG pathway to the total number of differential genes. The ordinate is the KEGG pathway. The size of the dot represents the number of genes annotated to the KEGG pathway. And the level of enrichment varies from purple to red. After the quantitative analysis of gene expression, it is necessary to conduct statistical analysis on the expression data, and screen the genes with significantly different expression levels in different states. Some differentially expressed genes are shown in Table 3.
The fine type and coarse type LCGs were compared by histogram and volcanic map, we can see that there were 27197 coexpressed genes and 529 differential genes, of which 343 were upregulated and 186 were downregulated. The clustering graph compared the gene expression of three samples of coarse LCG and three samples of fine LCG. Therefore, it can be seen that the overall gene expression trend of the three samples of coarse LCG was significantly different from that of the three samples of fine LCG. The genes with high expression in fine LCG had lower expression in coarse LCG, and the genes with low expression in fine LCG had higher expression in coarse LCG ( Figures 4A-C).

GO Enrichment Analysis and KEGG Enrichment Analysis of Differential Genes
We Used ClusterProfiler Software (Yu et al., 2012) for analysis, through GO enrichment analysis, we selected the first 30 for analysis, and found 491 significantly enriched genes, including biological process (BP, 257), cell composition (CC, 64), and molecular function (MF, 170). The function was mainly in the extracellular region ( Figures 5A,B). From the results of KEGG enrichment, the most significant 20 KEGG pathways were selected, and the most enriched pathway was the human papillomavirus infection pathway ( Figures 5C,D).

P-Site Analysis of Six LCGs
During translation, ribosomes move in the unit of codon length (3 nt) relative to RNA. Therefore, based on P-site, RPF fragments from normal translation should be periodically distributed on RNA. The translation starts at 12 nt upstream of the initiation codon, and the distance from the termination codon 15 nt disappears gradually. This is direct evidence of whether an RNA is translated. The length of RNA fragments protected by the ribosome from the initial codon to the termination codon of three samples of LCG was 26, 27, and 28 nt, respectively ( Figure 6).

DISCUSSION
LCG is unique in China, which produces a large number of highquality cashmere fibers. Moreover, China is one of the largest cashmere producing countries in the world, which has made great contributions to the fiber industry and textile industry, and plays an indispensable role in the global cashmere industry . But at present, the pursuit of cashmere fineness is increasing, and the cashmere fineness of LCG is still showing a relatively coarse trend, and quality cashmere products are still insufficient. Reducing cashmere fineness is an important issue . Ribosome profiling is a mature method to identify translation regions by high-throughput sequencing, which fills the gap between RNA sequencing and proteomics, and has become a powerful tool for gene expression (Calviello and Ohler, 2017). Ribo-seq can not only measure the translation efficiency according to the relative abundance of ribosomes on transcripts, but also reveal the dynamic and local regulation of different translation stages according to the location information of footprints on transcripts (Li et al., 2020).
However, little is known about the issue of cashmere fineness of LCG by Ribo-seq sequencing in translatomics. In this study, we selected six adult female Liaoning Cashmere Goats (LCGs) with different cashmere fineness (divided into two groups). The coarse Liaoning Cashmere Goat sample group was the reference group, and the fine Liaoning Cashmere Goat sample group was the experimental group. The classification of groups was based on the phenotypic determination of cashmere fineness (cashmere fineness analyzer, sirolan technology, Australia). A total of 529 differentially expressed genes were identified by Ribo-seq sequencing, of which 343 were upregulated and 186 were downregulated. And the enrichment length of mRNA fragments was 22 nt.
COL6A5 (formerly known as COL29A1) is a member of the collagen superfamily. The gene is located on the long arm of chromosome 3 (Strafella et al., 2019), with a length of 8742. It is a protein-coding gene and is considered to be an extracellular matrix protein, accounting for 30% of the total mammalian protein (Haq et al., 2019). COL6A5 was found in the epithelial tissues of lung and gastrointestinal tract (Fitzgerald et al., 2008), but the highest expression was found below the dermal epidermal junction and around the reticular dermal vessels (Sabatelli et al., 2011). It was found that COL6A5 fibroblasts existed in atopic dermatitis skin, but not in healthy tissues (He et al., 2020). The COL6A5 gene is not only associated with skin inflammation, but also with cancer. It has been confirmed that the COL6A5 gene is significantly associated with the overall survival rate of patients with esophageal squamous cell carcinoma (ESCC). The overall survival rate of ESCC patients with low expression of the COL6A5 gene is poor, which may become a diagnostic marker of ESCC as a collagen gene (Li et al., 2019). Abundant type ⅵcollagen in lung tissue α5 (COL6A5), rs13062453, rs1497305, and rs77123808 of COL6A5 polymorphism are associated with lung cancer risk in Chinese Han population, and the overall survival rate of patients with low expression of the COL6A5 gene is poor (Duan et al., 2020). These studies can infer that the COL6A5 gene may play a role in changing the hair follicle and cashmere fineness of LCG. Because of that and because another gene from the same family, COL1A1, is known to have an impact on changing cashmere fineness (Wang et al., 2021), we hypothesize that COL6A5 is a candidate gene for future studies regarding cashmere fineness. Studies have shown that COL6A5 is associated with familial chronic neurotrophic itching (Martinelli-Boneschi et al., 2017). It Translatomics Ribo-seq Sequencing has been found that COL6A5 expression in the papillary dermis and the surrounding dermis of the patients is reduced (Weisshaar and Dalgard, 2009;Ständer et al., 2010), and the incidence rate increases with age. This is the first time that the link between the COL6A5 gene and chronic pruritus has been revealed. Some studies have found a link between COL6A5 variants, reduced bone mass, dyspnea, and giant cell arteritis . In psoriasis (Ps) and psoriatic arthritis (PsA), bioinformatics analysis revealed that COL6A5 and COL8A1 participate in the altered proliferation and angiogenesis pathways in Ps/PsA, participate in inflammatory response together with miR-146a, and participate in the common and different biological pathways of Ps and PsA (Caputo et al., 2020). The collagen gene may be closely related to the PI3K/Akt/ mTOR pathway, p53 pathway, apoptosis, and cell cycle. COL1A1 and FGF10 genes are also enriched in the PI3K/Akt/mTOR pathway. COL1A1 can regulate the growth of alpaca wool fiber, and FGF10 can prolong the growth period of mouse hair follicles and promote hair growth (Maiese, 2015;Huang et al., 2017;Mendoza et al., 2019). Therefore, we take COL6A5 gene as a candidate gene for cashmere fineness research in the future.
In living organisms, ribosomes synthesize proteins in the process of translation, and translation regulation itself goes beyond the three processes of transcription, mRNA degradation, and protein degradation. Like other omics, translatomics analyzes all components in the translation process, and also includes the study of mRNAs, ribosomes, tRNAs, regulatory RNAs, and newborn peptide chains (Zhao et al., 2019). Meanwhile, the correlation between the transcriptome and proteome is usually very poor in omics data, because the phenomenon of post transcriptional regulation, translation control, and other factors such as frameset translation is common. Translational sequencing can accurately quantify the genes being translated, and indirectly detect the protein expression from the genomic level, indicating the real situation of gene transcription and expression in biological samples. In addition, by comparing the gene translation differences between different samples, we can reveal the molecular mechanism of related physiological and pathological differences. Translatomics is a bridge between transcriptomics and proteomics. The multi omics analysis of translatomics can better study the mechanism of translation regulation. By analyzing the correlation between translatomics and transcriptomics, we can study the change of translation rate within genes, compare the relationship between gene transcription and translation, and study the difference of gene translation efficiency under different physiological and pathological conditions by calculating gene translation efficiency, to explain its molecular mechanism. The association analysis of translatomics and proteomics can study the relationship between transient translation and protein accumulation, assist proteome identification, provide evidence of gene translation, and indirectly determine the proteins expressed in biological samples.
It is critical to accurately investigate the underlying mechanisms of miRNA translation inhibition, and analyze the effect of post transcriptional regulation and RNA modification on gene translation. It is generally believed that the gene that can encode a protein with a length of more than 100 amino acids is a protein coding gene, while other gene sequences are noncoding sequences. However, in recent years, more studies have shown that some RNA regions (including lncRNA, 5′UTR, 3′UTR, circRNA, etc.) that are traditionally considered not to encode proteins can translate some peptides with a length of less than 100 amino acids. These peptides, of less than 100 amino acids in length, also play a variety of important roles in organisms, including ontogeny, muscle contraction, and DNA repair. Because translatomics is used to sequence and quantify the RNA molecules being translated, it can provide direct translation evidence for these noncoding sequences, and help to find and identify new unknown peptides.

CONCLUSION
In conclusion, this study analyzed the process of screening cashmere fineness functional genes by translatomics through Ribo-seq sequencing, found that the COL6A5 gene may play an important role in cashmere fineness regulation, and provided some theoretical basis for future research on this gene in the field of cashmere fineness regulation.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GEO Database, GSE186959.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Experimental Committee of Shenyang Agricultural University.

AUTHOR CONTRIBUTIONS
YZ and DZ carried out data analysis and writing of the original draft. ZB, WC, and XZ were in charge of experiments. YX, MG, YQ, RC, YS, and YW were responsible for the collection of samples. ZW reviewed and edited the manuscript.

FUNDING
Our scientific research was financially aided by three projects:1. Postdoctoral Science Foundation of China: Genetic trajectory and expression localization of key genes of cashmere fineness by multiomics, Project No. 2021M693859. 2. 2021 Liaoning Province "the open competition mechanism to select the best candidates" Science and Technology Research Project: Selection and breeding of special advantageous livestock and poultry breeds in Liaoning and key technology of whole industry chain production, Project No. 2021JH1/10400033. 3. National modern agricultural industrial technology system, project No: cars-39-27.

ACKNOWLEDGMENTS
We are sincerely grateful to the coworkers who work inside our laboratory at Shenyang Agricultural University and Liaoning Cashmere Breeding base for their help in the collection of samples from the feedlot of Liaoning cashmere goats and in the massive analysis of the samples. We would also like to express our gratitude to Novogene Biotechnology Co., Ltd for providing valuable help in sequencing and bioinformatics analysis.