Animals and sample collection
Nine Leizhou black goats were divided equally into three groups according to the age of 0, 3 and 6 months old. All of them were raised under the same conditions with free access to food and water in natural lighting. All animals were slaughtered in accordance with animal welfare procedures, and after slaughter, we collected longissimus dorsi (LD) muscle samples of three growth stage (M0, M3 and M6). The tissue samples were immediately frozen in liquid nitrogen and stored at -80℃ before the analysis.
RNA isolation, library preparation and sequencing
The total RNA was isolated with the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and treated with DNase I (Qiagen, Beijing, China). Then 1% agarose gel electrophoresis was used to assess the degradation and contamination of the RNA and RNA Nano6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) were used to check the integrity of RNA.
The total RNA (3 μg) was used as input material for each sample preparation. First, Epicentre Ribo-zeroTMrRNA Removal Kit (Epicentre, Madison, WI, USA) was applied to remove the ribosomal RNA (rRNA) and ethanol precipitation was used to clean up the rRNA-free residue. Subsequently, sequencing libraries were generated using the rRNA-depleted RNA with the NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina® (NEB, Beverly, MA, USA), according to the manufacturer’s instructions. Briefly, fragmentation was carried out in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNaseH). The second strand cDNA synthesis was then performed using DNA polymerase I and RNaseH. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After the adenylation of the 30-ends of the DNA fragments, NEBNext Adaptors containing a hairpin loop structure were ligated to prepare for the hybridization. To preferentially select cDNA fragments of 150–200 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Miami, FL, USA). Subsequently, 3 μL USER enzyme (NEB, Beverly, MA, USA) was incubated with size-selected, adaptor-ligated cDNA at 37℃ for 15 min, followed by 5 min at 95℃. A PCR amplification was then performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, the products were purified using the AMPure XP system (Beckman Coulter, Miami, FL, USA), and the library quality was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).
Data analysis
The constructed libraries were sequenced on an Illumina HiSeq 4000 platform, and 150-bp paired-end reads were generated. After removing the sequence containing adapter and the reads containing ploy-N and low-quality reads through in-house Perl scripts, clean data were obtained. All the downstream analyses were based on the clean data with high quality. To obtain the high-quality reads, we performed the following filtering process: we removed the reads containing more than 10% unknown nucleotides and the reads containing more than 50% low quality nucleotides with Phred with a quality under 20. Mapping to the reference genome was the next step. Reads that passed the quality control were then mapped to the Ovis aries reference genome (Oar_v3.1). The index of the reference genome was built using bowtie2 v2.2.8, and paired-end clean reads were aligned to the reference genome using HISAT2 (v2.0.4). HISAT2 was run with “–rna-strandness RF”, and other parameters were set as default. Next was the transcriptome assembly. The mapped reads of each sample were assembled by StringTie (v1.3.1).
Before the screening, Cu merge was used to create the set of transcripts. Then, the lncRNA screening was carried out through the following steps: Step1: select the number of transcripts with 2 exons; Step2: out of the results from step1, select the transcripts with a length >200 bp; Step3: annotate the above transcripts using the Cu compare software; Step4: calculate the expression level of each transcript by Cu quant, and select the transcripts with FPKM 0.1; Step5: coding the potential screening: the coding potential of the transcript was predicted by three softwares: CNCI (Coding-Non-Coding-Index) (v2), CPC (Coding Potential Calculator) (0.9-r2), and PFAM (Pfam Scan) (v1.3); the intersection of the transcripts without a coding potential screened through the above three softwares with the default parameter was predicted as the lncRNA dataset.
The Ballgown was utilized to perform the straightforward linear-model-based differential expression analyses within a default statistical modeling framework. The transcripts with p-adjust < 0.05 were assigned as being differentially expressed.For each lncRNA locus, the 100-kb upstream and downstream regions were chosen to screen the co-located genes through the UCSC Genome Browser.
GO and KEGG enrichment analysis
In this study, the mRNAs within a 100-kb window upstream or downstream of DE-s were served as a cis-target mRNAs dataset of DE-lncRNAs.
Gene Ontology (GO) enrichment analysis was used on the target gene candidates of differentially expressed mRNAs and lncRNAs. GO-seq based Wallenius non-central hyper-geometric distribution[33], which could adjust for gene length bias, was implemented for GO enrichment analysis.
KEGG[34] is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS[35] software to test the statistical enrichment of the target gene candidates in KEGG pathways.
Construction of DEmRNAs-DElncRNAs Interaction Network
The combinations with pearson correlation lower than 0.60 and negative correlation were excluded. The relationship between DE-lncRNAs and DE-mRNAs was visualized using Cytoscape (V3.5.1).
qRT-PCR Verification
The cDNA for qRT-PCR was synthesized using PrimeScript RT Reagent Kit With gDNA Eraser (TaKaRa, Dalian, China) and qRT-PCR was performed using 2×Ultra SYBR Green qPCR Mix (Cistro, Shanghai, China). Capra hircus β-actin served as the endogenous control for mRNA and lncRNA expression analyses.