Integrated analysis of lncRNA and mRNA reveals novel insights into cashmere fineness in Tibetan cashmere goats

Animal ethics

All animal experiments were strictly performed according to the guidelines established by the Animal Care and Use Committee of Xinjiang Academy of Animal Science (Approval number 2019009).

Sample collection

All sample Tibetan Cashmere goats were selected from the same flock and raised under the same nutritional conditions on the Tibet Naqu Nima Cashmere Goat Breeding Farm (longitude: 91°12′–93°02′E, latitude: 30°31′–31°55′N), China. We collected cashmere samples of 50 ewes at the age of 15-month-old in anagen stage from the left mid-side regions of each animal body, these ewes were artificially inseminated with fresh sperm from a single ram. The wool samples were sent to a commercial laboratory (Quality Standards Institute of Animal Husbandry, Xinjiang Academy of Animal Husbandry, Urumqi, China) for measuring mean fiber diameter (MFD) values. The measurements were following the standardized methods set by China Fiber Inspection Bureau (CFIB). According to the measured results, we selected 4 ewes with the lowest MFD values (average 12.04 ± 0.03 µm) as the experimental fineness type (F) group, and another 4 ewes with the highest MFD values (average 14.88 ± 0.05 µm) as the experimental coarse type (C) group, respectively. Scapular skin tissues were collected from these 8 Tibetan cashmere goats in vivo. From each piece of skin tissue, approximately two cm² was immediately frozen by liquid nitrogen and stored at −80 °C until analysis.

RNA isolation, library construction and sequencing

Total RNA was isolated using the RNeasy mini kit (Qiagen, Germany), after which the RNA concentration and quality were determined by the Qubit^®2.0 Fluorometer (Life Technologies, USA) and the Nanodrop One spectrophotometer (Thermo Fisher Scientific Inc., USA). The integrity of total RNA was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies Inc., USA), and samples with RNA integrity number (RIN) values above 7.0 were used for sequencing. One microgram of RNA was used as input material for the RNA sample preparations. Subsequently, 8 cDNA libraries for two types of skin samples were constructed using rRNA-depleted RNA with the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, USA) following the manufacturer’s recommendations, each sample type contained four biological replicates. The libraries were sequenced on an Illumina NovaSeq 6000 platform (Illumina, USA), and 150 bp paired-end reads were generated. The raw data were assessed for quality using FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Clean reads were obtained by removing adapters and poly-N>10%, and low-quality reads in the raw data using the Fastx_toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html). Then, clean reads were mapped to goat reference sequences (Capra hircus ARS1.93) using Hisat2 (Kim, Langmead & Salzberg, 2015). Reference genome and annotation files were downloaded from the Ensembl browser (https://www.ncbi.nlm.nih.gov/genome/?term=capra+hircus+ars1). According to the mapping results, Cufflinks (Trapnell et al., 2012) programs were used to assemble transcripts. The raw data has been deposited in the Sequence Read Archive (SRA) database (BioProject ID: PRJNA643003).

Identification of lncRNAs

Based on the splicing results of Cufflinks, lncRNA was selected as the final candidate lncRNA for subsequent analysis through a workflow: (1) The class_codes annotated by Cuffcompare (http://cole-trapnell-lab.github.io/cufflinks/cuffcompare/) (Trapnell et al., 2012) was used to screen the candidate transcripts. Only those annotated as “i”, “u” and “x” were retained, which represent potential novel intergenic, intronic and cis-antisense transcripts, respectively. Transcripts with class code “=” annotated by Cuffcompare were considered as known genes. (2) Transcripts with a transcript length > 200 bp and exon number > 2 were selected. (3) The read coverage of each transcript was calculated by Cufflinks, and the transcript with coverage > 3 in at least one sample was selected. (4) The class-code information originating from Cuffcompare was used to screen the candidate transcripts. Finally, the candidate lncRNA coding potential was predicted by the Coding-Non-Coding Index (CNCI) (score < 0) (Sun et al., 2013), Coding Potential Calculator (CPC) (0.9-r²) (score < 0) (Lei et al., 2007), Pfam-scan (E-value < 0.001) (Robert et al., 2016).

Analysis of differentially expressed (DE) lncRNAs and mRNAs

Gene abundance was expressed as fragments per kilobase of exon per million reads mapped (FPKM). Stringtie software (Kovaka et al., 2019) was used to count the fragment within each gene, and TMM algorithm was used for normalization. Differential expression analysis for mRNA and lncRNA were performed using R package edgeR (Robinson, McCarthy & Smyth, 2010). DE RNAs with p-values < 0.05 and —log₂FoldChange— ≥ 1, considered as significantly modulated between the two groups, were retained for further analysis.

Enrichment analysis and protein-protein interaction (PPI)

All DE lncRNA target genes and DE mRNAs were analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) pathway enrichment analyses using the GOseq R package (Young et al., 2010) (v1.18.0). The Benjamini–Hochberg (BH) method was used to adjust significant p-values. GO terms and pathways with p-values < 0.05 were considered to be significantly enriched. We investigated candidate gene-encoded functional PPI network by means of the STRING Genomics 11.0 database (Szklarczyk et al., 2015), and these interactions predicted active interaction sources from text mining, experiments, databases, co-expression, neighborhood, gene fusion, or cooccurrence. A global PPI network was built that is limited to interactions exhibiting high confidence, with scores over 0.7 and no more than 10 interactors. We further showed the network by Cytoscape 3.6.1 software (Shannon et al., 2003).

Interaction and co-expression analysis of lncRNAs and mRNAs

We constructed a lncRNA–target gene network combining two instances of DE lncRNAs with their corresponding target genes, and the target genes of DE lncRNAs were predicted in cis. We further investigated the relationship between lncRNA and mRNA expression levels based on pearson correlations. If the expression levels of two transcripts exceeded a preselected threshold and were similar in the Pearson analysis, they would be considered to have a co-expression relationship and connected linearly, indicating a tight correlation. lncRNA or mRNA importance in the network depends on the degree of correlation. The absolute value of Pearson’s correlation coefficient (r_P) ≥ 0.90 and p-value < 0.05 were reserved for further construction of the lncRNA and mRNA network. Two networks were created by Cytoscape 3.6.1 software (Shannon et al., 2003).

Quantitative real-time PCR validation

To identify the DE lncRNAs and DE mRNAs between fineness type and coarse type, we randomly selected six DE lncRNAs and six DE mRNAs to verify the expression levels by quantitative real-time PCR (qRT-PCR). The GAPGH gene was used as an internal control. Primers for real-time PCR were designed with Primer 5.0 (Table S1). qRT-PCR was carried out with a Bio-Rad CFX96 Real-Time PCR system (Bio-Rad, CA, USA) using TB Green^® Premix Ex Taq™ II (TAKARA, Dalian, China). The 2^−ΔΔCt method (Livak & Schmittgen, 2001) was used to analyze the relative expression levels.

Overview of high-throughput sequencing

We comprehensively profile the lncRNAs and mRNAs expressed in the skin of cashmere goats (Table S2). Ultimately, a total of 2,059 lncRNAs were identified, and we further analyzed their genomic context based on their orientation to local protein-coding genes and classified them into 66 exonic-sense lncRNAs, 50 intronic sense lncRNAs, 195 exonic-antisense lncRNAs, 55 intronic antisense lncRNAs, 92 bidirectional lncRNAs, and 1,601 long intergenic RNAs (lincRNAs). After excluding candidate lncRNAs with coding potential using the software CNCI, CPC and Pfam-scan protein domain analysis, 470 lncRNAs were identified (Fig. 1A), including 333 lincRNAs, 40 intronic lncRNAs and 97 anti-sense lncRNAs (Fig. 1B). In addition, 29,119 mRNAs were identified.

Comparison of genomic features between lncRNAs and mRNAs

We analyzed the characteristics of expression level, length, and exon number between lncRNAs and protein-coding genes to help us understand the structural features of these genes in cashmere goats (Fig. 2). The results showed that the expression levels of most detected lncRNAs were lower than those of protein-coding genes (Fig. 2A). LncRNAs tended to contain fewer exons than protein-coding transcripts, most of them were composed of 2 or 3 exons (Fig. 2B), and they were shorter than mRNAs due to having fewer exons (Fig. 2C).

Differential expression analysis

We identified 80 DE lncRNAs (Fig. 3A) and 384 DE mRNAs (Fig. 3B) between the fine-type cashmere and coarse-type cashmere skin tissues in Tibetan Cashmere goats. Of these, there were 36 upregulated and 44 downregulated lncRNAs (Fig. 3A), together with 180 upregulated and 204 downregulated mRNAs (Fig. 3B) shown in the Volcano plots. Cluster analysis of DE lncRNAs and DE mRNAs was depicted in a heatmap, respectively (Figs. 3C, 3D). The results showed similar expression patterns and relationships in each group.

The top 10 DE lncRNAs and their potential target genes in fine-type cashmere goats and coarse-type goats are shown in Table 1. A total of 285 targets of 80 DE lncRNAs were predicted by in cis and trans principles, including KRT26, KRT28, KRT39, IFT88 and COL4A3BP, which are related to the hair follicle fiber growth. We also found that some DE genes might be involved in skin and hair follicle differentiation and development, including KRT85, KRT33A, NOTCH2, NOTCH3, FOS, FOSB, CXCL9, CXCL911 and COL4A3BP, and the details are shown in Table 2. Of all the target genes and DE genes, a total of 5 common genes were obtained, including ENSCHIT00000009853 targets were PCCB, DLG1, SELE, COL4A3BP and MSTRG.11347.5 target was DCAF11.

GO and KEGG functional enrichment analysis of DE lncRNAs and DE mRNAs

The top10 significant GO terms (p-value < 0.01) of DE lncRNA target genes and DE mRNAs are shown in Figs. 4A, 4B, respectively. GO terms of these lncRNA potential targets were highly enriched for several developmental processes, such as positive regulation of GTPase activity, generation of precursor metabolites and energy, and microtubule cytoskeleton organization involved in mitosis (Fig. 4A). The DE mRNAs were highly enriched for regulation of protein processing, positive regulation of the canonical Wnt signaling pathway, positive regulation of actin filament polymerization, maintenance of localization in cells and smooth muscle contraction (Fig. 4B).

The significant KEGG pathways (p-value < 0.05) of DE lncRNA target genes and DE mRNAs are shown in Figs. 4C, 4D, respectively. For example, the propanoate metabolism, valine, leucine and isoleucine degradation, DNA replication and carbon metabolism were enriched in DE lncRNA target genes and might regulate hair follicle and skin development (Fig. 4C). KEGG analyses of DE mRNAs were also enriched in two of the four pathways mentioned above, including the propanoate metabolism pathway and valine, leucine and isoleucine degradation (Fig. 4D). PPI analyses were performed for the 265 most plausible candidate genes. The interaction network of proteins encoded by these genes was more extensive and significant than expected (72 edges identified; PPI enrichment p-value = 1.92 e−08; Fig. 5). The network revealed that PMSD3, PMSD11, PMSD14, and SRC were the core nodes. Several potential key DE mRNAs were enriched in the TNF signaling pathway (FOS, CXCL10), hippo signaling pathway (DLG1, DLG2, CDH1), adherens junction (SRC, CDH1, WASL), and focal adhesion (SRC, TLN1).

Network of lncRNAs and mRNAs

LncRNAs can regulate mRNAs by targeting their functions; therefore, we constructed an interaction network based on two interested DE lncRNAs and their target genes (Fig. 6A). ENSCHIT00000009853 targets were predicted to interact with keratins (KRT26, KRT28, KRT39), TNF signals (SELE, CASP7), notch signals (PSEN2), PI3K-Akt signals (PSEN2), and focal adhesion (VAV1), and MSTRG.16794.17 was predicted to interact with TNF signals (CASP7).

Meanwhile, we selected 65 significantly DE lncRNAs and 64 DE mRNAs to construct a co-expression network (Fig. 6B). This result demonstrates these DE lncRNAs were similar to the expression patterns of adjacent genes, constituted a complex regulatory network to regulate the cashmere fineness. For instance, MSTRG.17532.2 with high correlation with Notch2 and Nothch3.

Validation of lncRNAs and mRNAs by qRT-PCR

To validate the reliability of the RNA-seq data, six DE mRNAs, including COL4A3BP, SELE, KRT85, KRT33A, FOSB, and FOS, and six DE lncRNAs, including ENSCHIT00000009853, MSTRG.6712.3, MSTRG.5970.13, MSTRG.1727.25, MSTRG.17103.10 and MSTRG.9227.3, were randomly selected for qRT-PCR analysis. Figure 7 presents the expression levels of DE lncRNAs and DE mRNAs by qRT-PCR and were found to be consistent with the RNA-seq data. This indicated that our transcript identification and abundance estimation were highly reliable.