Genome-wide identification and characterization of long non-coding RNAs during postnatal development of rabbit adipose tissue

The rabbit is widely used as an important experimental model for biomedical research, and shows low adipose tissue deposition during growth. Long non-coding RNAs (lncRNAs) are associated with adipose growth, but little is known about the function of lncRNAs in the rabbit adipose tissue. Deep RNA-sequencing and comprehensive bioinformatics analyses were used to characterize the lncRNAs of rabbit visceral adipose tissue (VAT) at 35, 85 and 120 days after birth. Differentially expressed (DE) lncRNAs were identified at the three growth stages by DESeq. The cis and trans prediction ways predicted the target genes of the DE lncRNAs. To explore the function of lncRNAs, Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed on the candidate genes. A total of 991,157,544 clean reads were generated after RNA-Seq of the three growth stages, of which, 30,353 and 107 differentially expressed (DE) lncRNAs were identified. Compared to the protein-coding transcripts, the rabbit lncRNAs shared some characteristics such as shorter length and fewer exons. Cis and trans target gene prediction revealed, 43 and 64 DE lncRNAs respectively, corresponding to 72 and 20 protein-coding genes. GO enrichment and KEGG pathway analyses revealed that the candidate DE lncRNA target genes were involved in oxidative phosphorylation, glyoxylate and dicarboxylate metabolism, and other adipose growth-related pathways. Six DE lncRNAs were randomly selected and validated by q-PCR. This study is the first to profile the potentially functional lncRNAs in the adipose tissue growth in rabbits, and contributes to our understanding of mammalian adipogenesis.


Background
Obesity is becoming increasingly prevalent in both developed and developing countries, and the proportion of overweight and obese individuals is expected to reach 89 and 85% in men and women respectively,by 2030 [1,2]. Since excessive weight gain is associated with diseases of the metabolic syndrome, including hyperglycemia, dyslipidemia, hypertension, atherosclerosis [3,4], diabetes and cardiovascular disease, the global rise in obesity poses a considerable threat to human health [5,6]. Obesity manifests as over-accumulation of fat in white adipose tissue (WAT), due to an increase in the size and number of adipocytes [7]. WAT is categorized into two types based on its distribution: the visceral adipose tissue (VAT) which is located within specific regions of the abdominal cavity, and subcutaneous adipose tissue (SAT), which is present below the skin [8].
The amount of fat in both VAT and SAT is closely related to the incidence of metabolic diseases [9], and excessive fat accumulation in VAT in particular is regarded as a high risk factor for metabolic disorders and cardiovascular diseases [10][11][12]. Therefore, studies on visceral adipocyte differentiation and its potential regulatory mechanisms have long been the core of obesity research. The development and function of VAT is age dependent, as one study showed that protein synthesis in bovine adipose tissue was more active in fetuses than in adults [13], as were the transcription factors PPARG, C/EBP and STAT which are known to promote the development of WAT [14].
Rabbits are an economically important livestock animal, and are raised for meat and fur production, as well as experimental models for biomedical research. Rabbit meat is protein rich, and low in cholesterol, and fat. Since the adipose tissue in rabbits has a lower deposition rate during growth, the rabbit is an ideal model to study adipose regulation [28][29][30][31][32]. Several lncRNAs involved in adipogenesis in mammals like cattle [33], pigs [34] and chicken [35], have been identified but the adipose tissue growth-related lncRNAs in rabbits have not been profiled so far.
We used the RNA-seq based approach to determine lncRNAs expression levels in rabbit VAT at 35, 85 and 120 days after birth, and identified 30,353 lncRNAs of which 107 were differentially expressed indicating their potential in rabbit adipogenesis. Our findings provide further insights into the regulatory function of lncRNAs in rabbits and for annotating the rabbit genome, as well as contribute to better understanding of adipogenesis.

Animal and sample collection
Tianfu Black rabbits (native species in Sichuan province of China) aged 35, 85 and 120 days were used in this study. Given the plasticity and maturation of rabbit VAT [36], three biological replicates of peri-renal fat were collected for 35 days (YR) and 120 days (TR), and two for 85 days (MR). The samples were snap frozen in liquid nitrogen, and stored at − 80°C until RNA extraction.

Total RNA extraction
Total RNA was isolated using Trizol Reagent (Life Technologies, Carlsbad, CA, USA). The purity and integrity of the RNA were determined using Nanodrop (Thermo Fisher Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) respectively, and RNA concentration was measured using a Qubit® RNA Assay Kit and Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). Integrity of RNA was assessed using. Only samples that had RNA Integrity Number scores > 8 were used for sequencing.

Library construction and sequencing
The libraries construction and sequencing were performed by Mega Genomics Co.,Ltd., (Beijing, China). Briefly, 1 μg RNA was taken per sample and rRNA was removed using a NR603-VAHTS Total RNA-seq (HMR) Library Prep Kit (Vazyme Biotech Co.,Ltd., Nanjing, China). First-strand cDNA was synthesized using random hexamer primers, followed by the second-strand synthesis using DNA polymerase I and RNase H. The resulting double-stranded DNA was purified by AMPure XP beads, and a poly A tail was ligated to the sequencing joint. The correct-sized fragments were purified by AMPure XP beads. The USER enzyme was used to degrade the cDNA strands containing U instead of T, and the first-strand cDNA was sequenced, thereby preserving the direction of the RNA. Finally, PCR amplification was conducted and the products were purified (AMPure XP beads) for constructing the cDNA libraries.
The quality of the latter was assessed using Agilent BioAnalyzer 2100 system and qPCR. The libraries were sequenced on an Illumina HiSeq X Ten platform and 150-bp long paired-end reads were generated.

Transcriptome assembly and lncRNAs detection
To ensure the accuracy of information analysis, the adapter and low-quality reads were removed from the raw reads, and only the high quality clean reads were used for subsequent analysis. They were aligned to the rabbit reference genome (GCF_000003625.3_OryCun2.0_genomic.fa), along with annotated genes (GCF_000003625.3_Ory-Cun2.0_genomic.gff) using histat2 (2.0.5) software with the parameters '-dta-rf-p 1-x-1-1-S File for SAM output (default:stdout)' [37]. The StringTie program [38,39] was used to splice the Mapped Reads. LncRNA identification consisted of two steps: basic screening and potential coding ability screening. The transcripts longer than 200 bp, containing two or more exons, and fragments per kilobase of transcript per million fragments mapped (FPKM) [40] ≥ 0.1 were first selected by basic screening. Subsequently, CPC [41], CNCI [42], CPAT [43], and Pfam protein structure domain analysis [44] were used to screen for potential coding transcripts. Stringtie (1.3.3) was used to quantify transcripts and normalize the expression values (FPKM).

Screening of differentially expressed lncRNAs
Differentially expressed (DE) lncRNAs between any two libraries were identified by DESeq (1.26.0), with padjust < 0.01 and an absolute value of the |log2(fold change)| ≥ 2.0 or ≤ 1/2.0 as the threshold. All DE lncRNAs in eight libraries were clustered using the Heatmaps in R software package.

Target gene prediction and functional enrichment analysis
Most of the lncRNAs deposited in current databases have not yet been functionally annotated. Therefore, prediction of their functions is based on the functional annotations of their related cis and trans mRNAs. Coding genes located at the distance of 100-kb were considered potentially cis-regulated target genes, while LncTar [45] was used to predict the potentially trans-regulated target genes of the DE lncRNAs. GO enrichment and KEGG pathway analyses of the candidate DE lncRNA target genes were then performed using GOseq [46] and R package, respectively. Significance was calculated using the Expression Analysis method, and P value < 0.05 was considered significant.

Validation of DE lncRNAs by q-PCR
Primers for the lncRNAs and internal controls (Additional file 1) were designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). Total RNA was converted to cDNA using a Prime-Script™ RT Reagent Kit containing gDNA Eraser (TAKARA, Dalian, China), and oligo (dT) and random hexamer primers. The q-PCR was performed using SYBR Premix Ex Taq™ II (TAKARA) according to the manufacturer's instructions. The reaction mix consisted of 5 μl SYBR Premix Ex Taq™ II, 1 μl template cDNA, 0.4 μl of 10 μM forward and reverse primers, and 3.2 μl dH2O to a final volume of 10 μl. The reactions were performed on a Rotor gene 6000 PCR System (QIAGEN, Hiden, Germany) as follows: 95°C for 10s, followed by 40 cycles of 95°C for 5 s, and 20s at the Tm (Additional file 1). Melting curve analysis was performed from 65°C to 95°C with increments of 1.5°C. The expression levels of lncRNAs were normalized to HPRT1 and GAPDH. Relative gene expression levels were calculated using the 2 -ΔΔCt method [47], and data were expressed as mean ± standard error of the mean (SEM).

Statistical analysis
Statistical analysis was performed using the SPSS Statistics 20.0 (SPSS Inc., Chicago, IL, USA), and P < 0.05 was considered statistically significant.

Overview of RNA-Seq
To identify lncRNAs expressed during the growth of rabbit adipose tissue, we constructed 8 cDNA libraries (YR-1,YR-2,YR-3, MR-1,MR-2,TR-1,TR-2 and TR-3) from the peri-renal adipose tissues of 35 (YR), 85 (MR), and 120 (TR) day-old rabbits. The libraries were sequenced using the Illumina HiSeq X Ten platform, and 1,016,066,842 raw reads were generated. After filtering adaptor sequences and low-quality reads, we finally obtained 991,157,544 clean reads. The percentage of clean reads and GC content of libraries ranged from 97.29-98.12% and 49.5-51.65%. After mapping the clean reads to the rabbit reference genome, approximately 89.53% of the clean reads were selected for further analysis (Additional file 2).

Identification of DE lncRNAs
The expression levels of the lncRNAs were calculated by FPKM using DESeq. A total of 107 DE lncRNAs were identified (p < 0.05) during the growth of adipose tissue (Fig. 4), of which 42.7% were up-regulated and 57.3% were down-regulated. The DE lncRNAs with similar expression levels across the different libraries were screened and clustered using the systematic cluster analysis (Fig. 5a) (Fig. 5b).

Enrichment analysis of the target genes of DE lncRNAs
The potential cis-regulated target genes of lncRNAs were first predicted and 43 out of the 107 DE lncRNAs corresponded to 72 protein-coding genes. Gene ontology (GO) analysis [48] of the cis lncRNA targets showed significant enrichment of 192 GO terms (P < 0.05), of which U2-type spliceosomal complex (GO:0005684), nucleic acid transmembrane transporter activity (GO:0051032) and regulation of DNA ligation (GO:0051105) were the top listed GO terms involved in cellular component (CC), molecular function (MF) and biological process (BP) respectively (Fig. 6a). In addition, several terms associated with adipose accumulation like cell morphogenesis (GO:0000902), macromolecule metabolic process (GO:0043170) and lytic vacuole (GO:0000323) were also highly enriched. Furthermore, 20 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched, of which several were related to adipocyte growth such as oxidative phosphorylation (KO00190), glyoxylate and dicarboxylate metabolism (KO00630), and calcium signaling pathway (KO4323) (Fig. 6b). Interestingly, we also found target genes including GCSH, NDUFS5, HNRNPA3, GSN, CYSLTR2, and TFIP11, which were annotated with adipose growth-related GO terms and pathways.
Prediction of the trans-regulated target genes showed that, 64 DE lncRNAs corresponded to 20 protein-coding genes, with significant enrichment of 124 GO terms (P < 0.05), including carbon-oxygen lyase activity (GO:0016835), regulation of catabolic process (GO:0009894) and cell cycle process (GO:0022402) (Fig. 7a). KEGG analysis revealed 5 enriched pathways (Fig. 7b) including the Jak-STAT signaling (KO04630) and Purine metabolism pathways (KO00230). Taken together, the majority of target genes of Fig. 2 Comparison between rabbit lncRNAs and protein-coding genes. a Transcript size distribution of rabbit lncRNAs and protein-coding genes. b Exon numbers of rabbit lncRNAs and protein-coding genes

Validation of DE lncRNAs
To validate the RNA-Seq results, we randomly selected six DE lncRNAs and examined their expression patterns at the three growth stages by q-PCR. The results showed that all six lncRNAs (MSTRG.118114.1, MSTRG.93684.2, MSTRG.13808.1, MSTRG.20558.1, MSTRG159118.2, MSTRG.67280.1) were differentially expressed at different stages. In addition, the six lncRNAs exhibited a similar trend between the results of RNA-seq and q-PCR (Fig. 8). Therefore, the FPKM obtained from RNA-seq can be reliably used to determine lncRNAs expression.

Discussion
LncRNAs play important roles in genomic imprinting, cellular trafficking and organization, as well as regulation of gene expression and dosage compensation [49]. However, in contrast to mice [50,51] and humans [52], limited information is available regarding rabbit lncRNAs. We identified 30,353 lncRNAs in the rabbit peri-renal adipose tissues, which is more than the number reported for chicken [53,54] and pigs [34]. This may be the result of the differences between the species. To the best of our knowledge, this is the first report to systematically identify lncRNAs at the three growth stages of rabbit VAT: 35, 85, and 120 days after birth using RNA-Seq. Six of these DE lncRNAs were validated using q-PCR, and showed similar trends as the RNA-Seq results, indicating that our approach was reliable and accurately reflected the differential lncRNA expression during adipose growth in rabbit.
Compared to previous studies on other species, the rabbit lncRNAs were shorter, and had fewer exons, shorter ORFs and lower expression levels compared to the protein-coding genes [55][56][57]. These characteristics of lncRNAs are common across mammals, and are likely significant to their regulatory function, in addition to acting as a template for identifying the lncRNAs of different mammalian species.
In this study, we detected a total of 107 DE lncRNAs in the pairwise comparisons of different growth stages of rabbit VAT. The numbers of down-regulated lncRNAs Fig. 3 The distribution of open reading frame(ORF) size and expression level of the lncRNAs and protein-coding genes. a Distribution of ORF sizes in the lncRNAs and protein-coding genes. b Expression level analysis in the lncRNAs and protein-coding genes Fig. 4 Number of up-regulated and down-regulated lncRNAs in the rabbit peri-renal adipose at three growth periods was higher than that of the up-regulated lncRNAs, indicating that most lncRNAs were involved in the early adipose growth. In addition, 54 lncRNAs were differentially expressed in both the YR vs MR and TR vs MR comparisons, although their expression levels were different. This indicates that the DE lncRNAs likely play different role in the different stages of VAT growth. Taken together, the DE lncRNAs identified in our study are important regulators of rabbit adipose growth.
Our and others' studies were not able to infer the function of lncRNAs from their sequences and structures [58,59]. LncRNAs regulate the expression of nearby or distal genes, respectively known as cis or trans regulation [60]. We predicted the potential cis and trans regulated target genes of the 107 DE lncRNAs, and determined their potential biological relevance using GO and KEGG analysis.
The cis-regulated target genes showed significant enrichment in GO terms associated with lipid metabolism and functions of adipose tissue, such as mRNA metabolic process, cellular protein modification process, cell part morphogenesis, mast cell cytokine production, histone H2B ubiquitination, response to corticotropin-releasing hormone, mast cell activation involved in immune response, macromolecule metabolic process, growth and lytic vacuole organization. Pathway analysis showed that the cis-target genes of the DE lncRNAs were mainly involved in oxidative phosphorylation, glyoxylate and dicarboxylate metabolism, glycine, serine and threonine metabolism, and calcium signaling pathway. Calcium is a key intracellular signaling mediator regulating numerous cellular processes [61]. Activation of CaSR, an extracellular Ca 2+ sensor, in the visceral white adipose progenitor cells increases proliferation and adipocyte differentiation [62]. Some of the cis target protein-coding genes are  and NEAT1 directly or indirectly regulate expression of PPARγ [63][64][65][66]. The PU.1 antisense lncRNA was discovered within the mouse PU.1 locus [67], and promoted adipogenesis by transcriptionally repressing the adjacent PU.1 mRNA [63]. Some lncRNAs can also regulate their target genes in the trans mode [68], and the trans-regulated target genes predicted in this study were enriched in MF terms like immune system, signal transduction, and nucleotide metabolism and BP terms including cell cycle process, coenzyme metabolic process, regulation of cell differentiation and regulation of multicellular organismal development. Some of the trans target protein-coding genes were involved in VAT growth such as IL10, PDE6D, MAFB, and PCBD2. Large scale genomic studies in recent years have furthered our understanding of lncRNAs as functional adipogenic regulators [69]. IL-10, an anti-inflammatory cytokine is known to contribute to childhood obesity [70,71], and Liu et al found that IL-10 and the downstream JAK-STAT pathway were down -regulated in obese children with hypertriglyceridemia and in HFD obese rats [72]. This, indicated a protective effect of IL-10 on lipid metabolic disorders, especially hypertriglyceridemia. Therefore, the DE lncRNAs might trans-regulate adipogenic differentiation and lipid metabolism by targeting genes involved in the relevant signaling pathways. In conclusion, we identified lncRNAs that regulate adipose growth and apoptosis through cis or trans-acting mechanisms. In future studies, we will investigate the functions of some of these DE lncRNAs in order to elucidate the regulatory mechanisms and potential therapeutic targets for obesity.

Conclusions
This is the first report of the lncRNA profile of rabbit VAT at different stages of growth, which identified 107 DE lncRNAs associated with adipogenetic pathways like oxidative phosphorylation, cell part morphogenesis and glyoxylate and dicarboxylate metabolism. These DE lncRNAs are thus involved in the development and metabolism of adipose tissue in rabbits.