Totally nine Maiwa yaks from Yak Technology Park (Hongyuan County of Sichuan Province in China) at diverse age conditions (1 day, 15 months and 5 years old) were used in this study and each developmental age stage had three yaks. The experimental animals were healthy and under the same management. Yaks were weighing 10.55 ~ 14.24 kg, 96.38 ~ 101.37 kg and 240.73 ~ 296.36 kg. All yaks were stunned with a captive bolt pistol (Cash 8000 Model Stunner, 0.22 calibre, 4.5 grain cartridge) to ameliorate the suffering of the animals prior to their humane killing, following which exsanguination via a transverse incision of the neck was carried out in the slaughterhouse. Then, the liver tissues were excised immediately and rapidly stored in liquid nitrogen until RNA isolation. Yaks with different age were established for profiling the transcriptomes in liver tissue.
RNA quantification and qualification
RNAiso Plus (TaKaRa, Japan) was used to extract total liver tissue RNA. The extent of RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). If the OD260/280 value of total liver tissue RNA is between 1.8 ~ 2.2, the quality inspection is qualified (Supplementary Table 1) and the next analysis can be carried out.
Library preparation for high-throughput sequencing
A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. A chain-specific library was constructed by removing rRNA, that was, sequencing libraries were generated using NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) . Briefly, mRNAs were purified from total RNA using poly-T oligo-attached magnetic beads. First strand cDNA was synthesized using random hexamer primer and RNase H. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. The purified double-stranded cDNA was subjected to end repair, was added poly-A tail, and was connected to a sequencing adapter. In order to select cDNA fragments of preferentially 150~200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3μL USER enzymes (NEB, USA) were used to degrade the second strand of U-containing cDNA, and perform PCR amplification to obtain the library. After the library was constructed, Qubit2.0 Fluorometer (Invitrogen, USA) was used for preliminary quantification. PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. qRT-PCR (Bio-Rad, USA) was used to accurately quantify the effective concentration of the library to ensure library quality. The library preparations were sequenced on an Illumina Hiseq platform and 150 bp paired-end reads were generated.
Transcriptome sequencing data analysis
The quality-controlled Clean Reads was compared with the yak reference genome (BosGru_v2.0: ftp://ftp.ensembl.org/pub/release-97/fasta/bos_mutus/) quickly and accurately by the HISAT2 v2.0.5 software (http://ccb.jhu.edu/software/hisat2) to obtain Mapping Reads for subsequent analysis . At the same time, the comparison results were evaluated for quality. By analyzing the different regions and chromosome distributions of Mapping Reads in the reference genome, to obtain the comparison efficiency and mapping information about Mapping Reads for each sample . The new transcripts for Mapping Reads were assembled and quantified by StringTie (1.3.3b) software . Sequencing depth and gene length were corrected using FPKM . Differential expression analysis was performed using the DESeq2 R package (1.10.1). DESeq2 provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. Genes with an adjusted padj<0.05 found by DESeq2 were assigned as differentially expressed .
ClusterProfiler (3.4.4) software was used to perform GO function enrichment analysis and KEGG pathway enrichment analysis on differential gene sets, to predict the biological processes and functions that they may participate in, and to make corresponding classification and statistics . All lncRNAs were used for target gene prediction, that was, the target genes of lncRNAs were predicted by the positional relationship (co-location) and expression correlation (co-expression) of lncRNAs and protein-coding genes . Then, functional enrichment analysis (GO / KEGG) was performed on the target genes of differential lncRNAs to predict the main function of lncRNAs.
Real-time fluorescence quantitative PCR
Measurement of gene expression with qRT-PCR has been applied in our studies. Briefly, total RNA was extracted from liver tissue of each group using a RNAiso Plus (TaKaRa, Japan). RNA was complete without degradation, good quality and high purity, which meet the requirements of subsequent experiments. qRT-PCR was performed in a total volume of 10 µL containing 5.2 µL of TB GreenTM Premix Ex TaqTM Ⅱ(TaKaRa, Japan), 1 µL of cDNA, 0.8 μL of each primer (10 mM) and 2.2 µL of double-distilled water. The reaction condition used were as follows: 95 °C for 3 min, followed by 39 cycles of 95 °C for 10 s, 53.4/52 °C (β-actin/verification genes) for 20 s and 72 °C for 20 s, with the dissolution curve increasing from 0.5 °C to 95 °C every 5 s. The assays were performed on a real-time fluorescence quantitative PCR System (Bio-Rad, USA). For each sample, the cycle threshold (CT) values were obtained from three replicates. β-actin mRNA was employed as an internal reference. The primers used for amplification of target and internal reference genes were presented in Supplementary Table 2. The relative expression levels of target genes were analyzed using the 2-ΔΔCT method.