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 (Fu et al., 2020). 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 (Wang et al., 2020). 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.
Materials and Methods
The whole process of experiments was based on guidelines from the Animal Experimental Committee of Shenyang Agricultural University (Shenyang, China, 201906099).
The sample collection was a crucial step for Ribo-seq since it was the starting point of library construction. Two groups of three Liaoning Cashmere Goats (LCGs) differing in cashmere fineness were used here, including three fine LCGs (14.32, 14.69, and 14.77 μm) and three coarse LCGs (17.23, 17.63, and 17.91 μm). These 2-year-old adult female Liaoning Cashmere Goats were reared under the same (sheep house, environment, management, nutrition) conditions. Skin samples were collected from the front edge of the left scapula, with a sample size of 2 cm2. Anesthetics were used to relieve the pain of goats.
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.
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.
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 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 (Ahmed et al., 2019). Novogene used the Ribocode package to analyze the P-site using Ribo-seq data (Xiao et al., 2018).
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.
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 (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).
FIGURE 1. Statistical chart of clean reads length enrichment of Ribo-seq in LCG. The abscissa is the enrichment length and the ordinate is the percentage of enrichment frequency. (A) Sample of fine LCG. (B) Sample of coarse LCG.
FIGURE 2. Distribution of reads in different regions of the reference genome. Red indicates the coding area. Purple indicates the UTR area. Green indicates the intergenic area. Blue indicates the intron area. (A) Sample of fine LCG. (B) Sample of coarse LCG.
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).
FIGURE 3. Distribution of gene expression in fine and coarse type LCGs. The abscissa in the graph is the sample name, and the ordinate is log2 (RPKM+1). The fine LCG sample is shown in pink, and the coarse LCG sample is shown in blue. The box graph of each region has five statistics (the maximum, upper quartile, median, lower quartile, and minimum from top to bottom).
The Statistics of Differential Genes Were Carried out, and the Volcanic Map and Cluster Map Were Drawn
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).
FIGURE 4. Statistical chart of the number of different genes between fine type and coarse type LCG. (A) The abscissa of the histogram is divided into two groups: fine type and coarse type, the ordinate is the gene counts. (B) In the volcanic map, the abscissa represents the log2foldchange of gene expression in the treatment and control groups, and the ordinate represents the significant level of gene expression difference between the treatment and control groups, in which the upregulated genes are indicated by red dots and the downregulated genes are indicated by green dots. (C) Gene expression cluster map of fine and coarse type LCG. In the graph, red indicates high gene expression and blue indicates low gene expression. The abscissa is the sample name, and the ordinate is the normalized value of RPKM.
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).
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.
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).
FIGURE 6. P-site map of fine type and coarse type LCG. The abscissa is the distance from the start code or stop code, and the ordinate is the reads on the comparison. (A–C) P-site map of three fine type samples of LCG. (D–F) P-site map of three coarse samples of LCG.
ORF Analysis of Six Samples of LCGs
In the samples of fine LCG, 9925 genes were transformed into protein during the translation process, accounting for 95.5% of the total. Among them, 47 uORF (0.5%) and 9 dORF (0.1%) could be translated into protein. There were 69 overlapping genes (0.7%), 286 novel-PCGs (2.8%), and 59 novel-NonPCGs (0.6%) (Figure 7A). There were 8540 genes (96.3%) transformed into protein in coarse LCG samples during translation, among which 23 genes (0.3%) were uORF and 6 genes (0.1%) were dORF. There were 35 overlapping genes (0.4%), 229 novel-PCGs (2.6%), and 36 novel-NonPCGs (0.4%) (Figure 7B).
LCG is unique in China, which produces a large number of high-quality 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 (Zheng et al., 2020). 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 (Zheng et al., 2019). 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 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 (Wang et al., 2016). 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.
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.
The animal study was reviewed and approved by Animal Experimental Committee of Shenyang Agricultural University.
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.
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 multi-omics, 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.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
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.
Aeschimann, F., Xiong, J., Arnold, A., Dieterich, C., and Großhans, H. (2015). Transcriptome-Wide Measurement of Ribosomal Occupancy by Ribosome Profiling. Methods 85, 75–89. doi:10.1016/j.ymeth.2015.06.013
Ahmed, N., Sormanni, P., Ciryam, P., Vendruscolo, M., Dobson, C. M., and O’Brien, E. P. (2019). Identifying A- and P-Site Locations on Ribosome-Protected mRNA Fragments Using Integer Programming. Sci. Rep. 9 (1), 6256. doi:10.1038/s41598-019-42348-x
Bai, W. L., Dang, Y. L., Wang, J. J., Yin, R. H., Wang, Z. Y., Zhu, Y. B., et al. (2016). Molecular Characterization, Expression and Methylation Status Analysis of BMP4 Gene in Skin Tissue of Liaoning cashmere Goat during Hair Follicle Cycle. Genetica 144 (4), 457–467. doi:10.1007/s10709-016-9914-1
Bazzini, A. A., Lee, M. T., and Giraldez, A. J. (2012). Ribosome Profiling Shows that miR-430 Reduces Translation before Causing mRNA Decay in Zebrafish. Science 336 (6078), 233–237. doi:10.1126/science.1215704
Cai, Y., Fu, W., Cai, D., Heller, R., Zheng, Z., Wen, J., et al. (2020). Ancient Genomes Reveal the Evolutionary History and Origin of Cashmere-Producing Goats in China. Mol. Biol. Evol. 37 (7), 2099–2109. doi:10.1093/molbev/msaa103
Caputo, V., Strafella, C., Termine, A., Campione, E., Bianchi, L., Novelli, G., et al. (2020). RNAseq-Based Prioritization Revealed COL6A5, COL8A1, COL10A1 and MIR146A as Common and Differential Susceptibility Biomarkers for Psoriasis and Psoriatic Arthritis: Confirmation from Genotyping Analysis of 1417 Italian Subjects. Int. J. Mol. Sci. 21 (8), 2740. doi:10.3390/ijms21082740
Courtes, F. C., Lin, J., Lim, H. L., Ng, S. W., Wong, N. S. C., Koh, G., et al. (2013). Translatome Analysis of CHO Cells to Identify Key Growth Genes. J. Biotechnol. 167 (3), 215–224. doi:10.1016/j.jbiotec.2013.07.010
Dai, B., Zhang, M., Yuan, J.-L., Ren, L.-Q., Han, X.-Y., and Liu, D.-J. (2019). Integrative Analysis of Methylation and Transcriptional Profiles to Reveal the Genetic Stability of Cashmere Traits in the Tβ4 Overexpression of Cashmere Goats. Animals 9 (12), 1002. doi:10.3390/ani9121002
Duan, Y., Liu, G., Sun, Y., Wu, J., Xiong, Z., Jin, T., et al. (2020). Collagen Type VI α5 Gene Variations May Predict the Risk of Lung Cancer Development in Chinese Han Population. Sci. Rep. 10 (1), 5010. doi:10.1038/s41598-020-61614-x
Fu, X., Zhao, B., Tian, K., Wu, Y., Suo, L., Ba, G., et al. (2020). Integrated Analysis of lncRNA and mRNA Reveals Novel Insights into cashmere Fineness in Tibetan cashmere Goats. PeerJ 8, e10217. doi:10.7717/peerj.10217
Gerashchenko, M. V., Lobanov, A. V., and Gladyshev, V. N. (2012). Genome-Wide Ribosome Profiling Reveals Complex Translational Regulation in Response to Oxidative Stress. Proc. Natl. Acad. Sci. 109 (43), 17394–17399. doi:10.1073/pnas.1120799109
Gross, S. M., Martin, J. A., Simpson, J., Abraham-Juarez, M. J., Wang, Z., and Visel, A. (2013). De Novo transcriptome Assembly of Drought Tolerant CAM Plants, Agave Deserti and Agave Tequilana. BMC Genomics 14, 563. doi:10.1186/1471-2164-14-563
He, H., Suryawanshi, H., Morozov, P., Gay-Mimbrera, J., Del Duca, E., Kim, H. J., et al. (2020). Single-Cell Transcriptome Analysis of Human Skin Identifies Novel Fibroblast Subpopulation and Enrichment of Immune Subsets in Atopic Dermatitis. J. Allergy Clin. Immunol. 145 (6), 1615–1628. doi:10.1016/j.jaci.2020.01.042
Hou, J., Wang, X., McShane, E., Zauber, H., Sun, W., Selbach, M., et al. (2015). Extensive Allele‐Specific Translational Regulation in Hybrid Mice. Mol. Syst. Biol. 11 (8), 825. doi:10.15252/msb.156240
Hui, T. Y. (2020). Correlation Analysis of cashmere Goat Skin ceRNA Regulatory Network and KRT26 and cashmere Fineness [D]. China: Shenyang Agricultural Uinversity. doi:10.27327/d.cnki.gshnu.2020.000662
Jia, B., Xi, J. F., Zhang, S. Y., Zhao, Z. S., Zhao, R. Q., and Chen, J. (2006). The Developmental Patterns of GH-R, IGF-1 and IGF-IR Gene Expression in Sheep Skin. Yi Chuan 28 (9), 1078–1082. Chinese. doi:10.1360/yc-006-1078
Jin, M., Wang, L., Li, S., Xing, M. X., and Zhang, X. (2011). Characterization and Expression Analysis of KAP7.1, KAP8.2 Gene in Liaoning New-Breeding cashmere Goat Hair Follicle. Mol. Biol. Rep. 38 (5), 3023–3028. doi:10.1007/s11033-010-9968-6
Jin, M., Cao, Q., Wang, R., Piao, J., Zhao, F., and Piao, J. (2017). Molecular Characterization and Expression Pattern of a Novel Keratin-Associated Protein 11.1 Gene in the Liaoning cashmere Goat (Capra Hircus). Asian-australas J. Anim. Sci. 30 (3), 328–337. doi:10.5713/ajas.16.0078
Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for Linking Genomes to Life and the Environment. Nucleic Acids Res. 36 (Database issue), D480–D484. doi:10.1093/nar/gkm882
Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: Accurate Alignment of Transcriptomes in the Presence of Insertions, Deletions and Gene Fusions. Genome Biol. 14 (4), R36. doi:10.1186/gb-2013-14-4-r36
Li, J., Wang, X., Zheng, K., Liu, Y., Li, J., Wang, S., et al. (2019). The Clinical Significance of Collagen Family Gene Expression in Esophageal Squamous Cell Carcinoma. PeerJ 7, e7705. doi:10.7717/peerj.7705
Liu, G., Liu, R., Tang, X., Cao, J., Zhao, S., and Yu, M. (2015). Expression Profiling Reveals Genes Involved in the Regulation of Wool Follicle Bulb Regression and Regeneration in Sheep. Int. J. Mol. Sci. 16 (5), 9152–9166. doi:10.3390/ijms16059152
Mahata, B., Sundqvist, A., and Xirodimas, D. P. (2012). Recruitment of RPL11 at Promoter Sites of P53-Regulated Genes upon Nucleolar Stress through NEDD8 and in an Mdm2-Dependent Manner. Oncogene 31 (25), 3060–3071. doi:10.1038/onc.2011.482
Martinelli-Boneschi, F., Colombi, M., Castori, M., Devigili, G., Eleopra, R., Malik, R. A., et al. (2017). COL6A5 Variants in Familial Neuropathic Chronic Itch. Brain 140 (3), 555–567. doi:10.1093/brain/aww343
McManus, C. J., May, G. E., Spealman, P., and Shteyman, A. (2014). Ribosome Profiling Reveals Post-Transcriptional Buffering of Divergent Gene Expression in Yeast. Genome Res. 24 (3), 422–430. doi:10.1101/gr.164996.113
Mendoza, M. N., Raudsepp, T., Alshanbari, F., Gutiérrez, G., and Ponce de León, F. A. (2019). Chromosomal Localization of Candidate Genes for Fiber Growth and Color in Alpaca (Vicugna Pacos). Front. Genet. 10, 583. doi:10.3389/fgene.2019.00583
Rufaut, N. W., Pearson, A. J., Nixon, A. J., Wheeler, T. T., and Wilkins, R. J. (1999). Identification of Differentially Expressed Genes during a Wool Follicle Growth Cycle Induced by Prolactin. J. Invest. Dermatol. 113 (6), 865–872. doi:10.1046/j.1523-1747.1999.00775.x
Sabatelli, P., Gara, S. K., Grumati, P., Urciuolo, A., Gualandi, F., Curci, R., et al. (2011). Expression of the Collagen VI α5 and α6 Chains in Normal Human Skin and in Skin of Patients with Collagen VI-Related Myopathies. J. Invest. Dermatol. 131 (1), 99–107. doi:10.1038/jid.2010.284
Sendoel, A., Dunn, J. G., Rodriguez, E. H., Naik, S., Gomez, N. C., Hurwitz, B., et al. (2017). Translation from Unconventional 5′ Start Sites Drives Tumour Initiation. Nature 541 (7638), 494–499. doi:10.1038/nature21036
Ständer, S., Schäfer, I., Phan, N. Q., Blome, C., Herberger, K., Heigel, H., et al. (2010). Prevalence of Chronic Pruritus in Germany: Results of a Cross-Sectional Study in a Sample Working Population of 11,730. Dermatology 221 (3), 229–235. doi:10.1159/000319862
Strafella, C., Caputo, V., Minozzi, G., Milano, F., Arcangeli, M., Sobhy, N., et al. (2019). Atopic Eczema: Genetic Analysis of COL6A5, COL8A1, and COL10A1 in Mediterranean Populations. Biomed. Res. Int. 2019, 3457898. doi:10.1155/2019/3457898
Tian, Q., Stepaniants, S. B., Mao, M., Weng, L., Feetham, M. C., Doyle, M. J., et al. (2004). Integrated Genomic and Proteomic Analyses of Gene Expression in Mammalian Cells. Mol. Cell Proteomics 3 (10), 960–969. doi:10.1074/mcp.m400055-mcp200
Wang, L., Feng, Z., Wang, X., Wang, X., and Zhang, X. (2010). DEGseq: An R Package for Identifying Differentially Expressed Genes from RNA-Seq Data. Bioinformatics 26 (1), 136–138. doi:10.1093/bioinformatics/btp612
Wang, X., Pandey, A. K., Mulligan, M. K., Williams, E. G., Mozhui, K., Li, Z., et al. (2016). Joint Mouse-Human Phenome-Wide Association to Test Gene Function and Disease Risk. Nat. Commun. 7, 10464. doi:10.1038/ncomms10464
Wang, S. H., Hsiao, C. J., Khan, Z., and Pritchard, J. K. (2018). Post-Translational Buffering Leads to Convergent Protein Expression Levels between Primates. Genome Biol. 19 (1), 83. doi:10.1186/s13059-018-1451-z
Wang, Z.-Y., Leushkin, E., Liechti, A., Ovchinnikova, S., Mößinger, K., Brüning, T., et al. (2020). Transcriptome and Translatome Co-Evolution in Mammals. Nature 588 (7839), 642–647. doi:10.1038/s41586-020-2899-z
Wang, Y., Li, G., Zhang, X., Zheng, Y., Guo, S., Guo, D., et al. (2021). Analysis of m6A Methylation in Skin Tissues of Different Sex Liaoning Cashmere Goats. Anim. Biotechnol. 25, 1–11. doi:10.1080/10495398.2021.1962897
Xiao, Z., Huang, R., Xing, X., Chen, Y., Deng, H., and Yang, X. (2018). De Novo Annotation and Characterization of the Translatome with Ribosome Profiling Data. Nucleic Acids Res. 46 (10), e61. doi:10.1093/nar/gky179
Yang, J.-B., Gan, S.-Q., Yang, Y.-L., Zhang, H.-L., Song, T.-Z., Feng, J., et al. (2012). Cloning and Expression in Follicle Anagen of ILK Gene in Sheep. Hereditas (Beijing) 34 (6), 719–726. Chinese. doi:10.3724/sp.j.1005.2012.00719
Yang, F., Liu, Z., Zhao, M., Mu, Q., Che, T., Xie, Y., et al. (2020). Skin Transcriptome Reveals the Periodic Changes in Genes Underlying cashmere (Ground Hair) Follicle Transition in cashmere Goats. BMC Genomics 21 (1), 392. doi:10.1186/s12864-020-06779-5
Yu, Z., Wildermoth, J. E., Wallace, O. A. M., Gordon, S. W., Maqbool, N. J., Maclean, P. H., et al. (2011). Annotation of Sheep Keratin Intermediate Filament Genes and Their Patterns of Expression. Exp. Dermatol. 20 (7), 582–588. doi:10.1111/j.1600-0625.2011.01274.x
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). ClusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. (clusterfile). doi:10.1089/omi.2011.0118
Zeng, X.-C., Chen, H.-Y., Jia, B., Zhao, Z.-S., Hui, W.-Q., Wang, Z.-B., et al. (2011). Identification of SNPs within the Sheep PROP1 Gene and Their Effects on Wool Traits. Mol. Biol. Rep. 38 (4), 2723–2728. doi:10.1007/s11033-010-0416-4
Zhang, B., Chang, L., Lan, X., Asif, N., Guan, F., Fu, D., et al. (2018). Genome-Wide Definition of Selective Sweeps Reveals Molecular Evidence of Trait-Driven Domestication Among Elite Goat (Capra Species) Breeds for the Production of Dairy, cashmere, and Meat. Gigascience 7 (12), giy105. doi:10.1093/gigascience/giy105
Zheng, Y. Y., Sheng, S. D., Hui, T. Y., Yue, C., Sun, J. M., Guo, D., et al. (2019). An Integrated Analysis of Cashmere Fineness lncRNAs in Cashmere Goats. Genes 10 (4), 266. doi:10.3390/genes10040266