De Novo Reconstruction of Transcriptome Identified Long Non-Coding RNA Regulator of Aging-Related Brown Adipose Tissue Whitening in Rabbits

Simple Summary Brown adipose tissues (BATs) undergo the conversion to white adipose tissues (WATs) with age. Long non-coding RNAs (lncRNAs) were widely involved in adipose biology. Rabbit is an ideal model for studying the dynamics of the transformation from BATs to WATs. However, our knowledge of lncRNAs that mediate the transformation remains unknown in rabbits. By histological analysis and sequencing, we found rabbit interscapular adipose tissues (iATs) from BATs to WATs within two years and identified a total of 631 differentially expressed lncRNAs (DELs) during the transformation process. Several signal pathways were involved in the transformation from BAT to WAT. A novel lncRNA that was highly expressed in iATs of aged rabbits was validated to impair brown adipocyte differentiation in vitro. Our study provided a comprehensive catalog of lncRNAs involved in the transformation from BATs to WATs in rabbits, facilitating a better understanding of adipose biology. Abstract Brown adipose tissues (BATs) convert to a “white-like” phenotype with age, which is also known as “aging-related BAT whitening (ARBW)”. Emerging evidence suggested that long non-coding RNAs (lncRNAs) were widely involved in adipose biology. Rabbit is an ideal model for studying the dynamics of ARBW. In this study, we performed histological analysis and strand-specific RNA-sequencing (ssRNA-seq) of rabbit interscapular adipose tissues (iATs). Our data indicated that the rabbit iATs underwent the ARBW from 0 days to 2 years and a total of 2281 novel lncRNAs were identified in the iATs. The classical rabbit BATs showed low lncRNA transcriptional complexity compared to white adipose tissues (WATs). A total of 631 differentially expressed lncRNAs (DELs) were identified in four stages. The signal pathways of purine metabolism, Wnt signaling pathway, peroxisome proliferator-activated receptor (PPAR) signaling pathway, cyclic guanosine monophosphate (cGMP)/cGMP-dependent protein kinase (cGMP-PKG) signaling pathway and lipid and atherosclerosis were significantly enriched by the DELs with unique expression patterns. A novel lncRNA that was highly expressed in the iATs of aged rabbits was validated to impair brown adipocyte differentiation in vitro. Our study provided a comprehensive catalog of lncRNAs involved in ARBW in rabbits, which facilitates a better understanding of adipose biology.


Introduction
The mammals contain white adipose tissues (WATs) that function as energy depositories and brown adipose tissues (BATs) that function in energy expenditure [1]. BAT is densely packed with mitochondria that express high levels of uncoupling protein 1 (UCP1), which allows proton leakage to uncouple respiration from ATP synthesis, playing an important role in thermogenesis for temperature homeostasis [1,2]. The BAT content is negatively correlated with the total fat deposition of the body [3,4]. In humans, BAT activity has beneficial metabolic effects on obesity, insulin resistance and atherosclerosis [5]. Although BAT can be found in most new-born animals, there is a difference in BAT development among different species [6][7][8]. The hibernating animal and some rodents persist in their BATs into adult life [9,10]. In most mammals such as ruminants, rabbits and humans, BATs undergo a poorly elucidated conversion to a "white-like" phenotype with age, which is known as "aging-related BAT whitening (ARBW)" [11]. However, the regulatory mechanisms underlying ARBW remain unknown.
Rabbit (Oryctolagus cuniculus) is an economically important domestic animal due to its high-quality meat, fur and hair [22,23]. Rabbits have low-fat deposition than other mammals, such as swine, cattle and sheep [24]. The BAT content and activity may account for the natural low-fat deposition of rabbits. Thus, rabbit, as an animal model, can be used in studying the dynamics of ARBW. The BAT development is important for the temperature homeostasis of newborn rabbits, especially for those were born in the cold ambient temperature [10]. Investigation of the molecular mechanisms underlying the process of ARBW could contribute to improving the welfare and survival rate of newborn rabbits. However, our knowledge of lncRNAs that mediate ARBW remains largely unknown in rabbits.
Interscapular adipose tissue (iAT) is the major BAT depot of rabbits [10]. In this study, we carried out a histological analysis to determine the process of ARBW of iATs in rabbits. We performed strand-specific RNA sequencing (ssRNA-seq) to identify lncRNAs involved in the process. Our data revealed the development of ARBW in rabbits and a total of 2281 novel lncRNAs were identified in our samples. Furthermore, classical rabbit BAT was a low lncRNA transcriptional complexity tissue. We identified lncRNAs that were significantly correlated with their flanking or reference genes. Clustering analysis of differentially expressed lncRNAs (DELs) revealed that lncRNAs with different expression patterns in stages played different roles in a cis-regulating way during ARBW. One DEL, MSTRG.2316.1, was validated to impair brown adipocyte differentiation in vitro. Our work is the first report of the dynamics of lncRNA regulatory mechanisms underlying the ARBW in rabbits and facilitates a better understanding of adipose biology.

Ethics Approval
All surgical procedures involving rabbits were performed according to the approved protocols of the Biological Studies Animal Care and Use Committee, Sichuan Province, China. Rabbits had free access to food and water under normal conditions and were humanely sacrificed as necessary to ameliorate suffering.

Tissue Sample Preparation, Histological Analysis and Immunohistochemistry (IHC)
In this study, the Tianfu Black rabbits (native species in Sichuan province of China) were raised at the breeding center of Sichuan Agricultural University, Ya'an, China. These rabbits were given ad libitum access to a standard diet and water as described previously [20]. To determine the lncRNA dynamics during BAT development in rabbits, a total of 12 samples were collected from iATs at four growth stages of 0 days (D0, infant stage), 15 days (D15, early whitening stage), 85 days (D85, puberty stage) and 2 years (Y2, aged stage) under sterile condition and 3 individuals were set at each stage. The samples for ssRNA-seq were snap-frozen in liquid nitrogen and stored at −80 • C until RNA extraction. The samples for histological assay were fixed using 4% paraformaldehyde at 4 • C overnight, embedded in paraffin and sliced. Haematoxylin and eosin (H&E) staining was carried out after deparaffination according to standard protocols as described in our previous study [25]. For IHC, the slices were incubated in the primary UCP1 antibody (1:100, Sangon Biotech, Shanghai, China) at 4 • C overnight and washed three times using phosphate buffer saline (PBS). Then, the slices were incubated in the secondary antibody (1:500, rabbit anti-mouse IgG, Sangon Biotech, Shanghai, China) for 45 min. All slices including in the assays of H&E staining and the IHC, were observed using an Olympus BX-50F light microscope (Olympus Optical, Tokyo, Japan).

SsRNA-seq and lncRNA Identification
The ssRNA-seq was done following our previous study [20]. Briefly, the total RNA of these samples was extracted using TRIzol reagent (Invitrogen, Hong Kong, China) and 1 µg RNA was used to construct a strand-specific library using the deoxyuridine triphosphate (dUTP) method. All purified libraries were sequenced on an Illumina NovaSeq 6000 platform. Finally, 150 bp oriented paired-end reads were generated. The quality of ssRNA-seq reads was checked using the Fastqc program (v0.11.8) [26]. Sequencing adapters and low-quality reads were removed using Cutadapt (v3.2) [27] software with parameters of '-a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCG-GAAGAGCGTCGTGTAGGGAAAGAGTGT -e 0.1 -m 100 -cut 0 -O 13'. Clean reads were aligned to the rabbit reference genome (OryCun2.0, Ensembl release 101) using HISAT2 (v2.1.0) [28] with the strand-specific parameter of '-rf' and the parameter of '-dta' for the downstream transcriptome reconstruction. Reconstruction of the transcriptome was conducted using Stringtie (v2.1.1) [29]. All transcripts generated from all samples were merged using Stringtie with a parameter of 'merge' to generate a consensus transcriptome. The transcripts with the length size of <200 bp were discarded and only transcripts containing multiple exons and having a value of transcripts per kilobase of exon model per million mapped reads >0.1 (TPM > 0.1) at least in one sample were collected to downstream analyses. To identify novel transcripts, the consensus transcriptome was compared to the rabbit reference transcriptome (OryCun2.0, Ensembl release 101) using Gffcompare (v0.11.2) [30]. The transcripts were located in the intergenic regions or a single intron, overlapping on an opposite strand of annotated transcripts and inadequately overlapping with an annotated exon were retained. For the retained transcripts, CPC2 (v0.1) [31], CPAT (v2.0.0) [32] and CNCI (2.0) [33] were used to check their protein-coding capacity and PfamScan [34] was used to align them to the Pfam database to remove potential protein-coding transcripts (http://pfam.xfam.org/, accessed on 10 November 2021). Only the non-coding transcripts identified by CPC2, CPAT, CNCI and PfamScan were considered credible lncRNAs. Most of the rabbit lncRNAs deposited in current databases had not yet been functionally annotated and lncRNAs shared poor conservation across species, which raised the challenges in inferring the functions of lncRNAs [35]. Recent studies demonstrated that lncRNAs could regulate the flanking genes to play their roles by cisregulation [36]. We performed lncRNA category-based cis-regulation prediction. For AS lncRNA, SO lncRNA and intronic lncRNA, the protein-coding genes located at the loci of the lncRNA were considered potential cis-regulated target genes. For lncRNA, the protein- coding genes located within 100 kb of lincRNA were considered potential cis-regulated target genes. The potential cis-regulated targets were used to predict lncRNA functions.

Transcriptomic Quantification and Differential Expression Analysis
To identify DELs during ARBW of the rabbit iAT, we quantified gene expression levels of lncRNAs. The stringtie '-eB' and '-A' were used to estimate raw read counts and TPM. The raw read counts were used as inputs to identify DELs in pairwise comparisons over the time courses of D0, D15, D85 and Y2 using DESeq2 [37]. The P values of hypothetical tests were adjusted using the method of false discovery rate (FDR). The lncRNAs with the thresholds of |log2(fold-change)| > 1.5 and FDR < 0.01 were considered DELs. Based on the TPM value, all DELs were clustered using R software (v4.4.1) and the clustering result was visualized using R package complexHeatmap [38]. Kyoto Encyclopedia of Genes and Genomes enrichment (KEGG) pathway analysis was conducted using R package clusterProfiler [39] and the enriched pathways that had a p value < 0.05 were considered significant.

Cell Culture and Plasmid-Mediated Overexpression
The brown preadipocytes (BPAs) were isolated from iATs of new-born Tianfu black rabbits using the collagenase I (Gibco, Carlsbad, CA, USA) method as described in our previous study [21]. The BPAs were placed into 12-well plates at a density of 3 × 10 5 cells per plate in the complete medium (Dulbecco's Modified Eagle Medium (DMEM) with high glucose, supplemented with 10% fetal bovine serum (FBS)) (Gibco, Carlsbad, CA, USA) and BPAs were placed in a humidified incubator at 37 • C and 5% CO 2 . Upon reaching approximately 80% confluence, the cells were induced using differentiation medium I (DMEM with high glucose, supplemented with 10% FBS, 500 µM 1-methy1-3-iosbutylxanthine [IBMX], 10 µg/mL insulin, 50 nM T3 and 5 µM dexamethasone) for 2 days. The cells were cultured in differentiation medium II (DMEM with high glucose, supplemented with 10% FBS, 500 µM IBMX, 2.5 µg/mL insulin, 50 nM T3 and 1 µM rosiglitazone) for 2 days. Finally, cells were cultured in differentiation medium III (DMEM with high glucose, supplemented with 10% FBS, 500 µM IBMX, 50 nM T3 and 1 µM rosiglitazone) for an additional 1 day. The cells were treated with 10 µM isoproterenol for 4 h before harvesting. The IBMX, insulin, T3, dexamethasone and rosiglitazone were purchased from Sigma-Aldrich (Shanghai, China). The accumulated lipid droplets were measured by Oil red O staining. The lipid-combined Oil red O dye was extracted by 2.5 mL isopropanol. The optical density (OD-510 nm) of the Oil red O elution was used to quantify the degree of lipid accumulation.
For lncRNA overexpression (OE), the vector of pcDNA3.1(+) loaded MSTRG.2316.1 and empty vector were purchased from Sangon Biotech, Shanghai, China. BPAs growing in complete medium at approximately 80% confluence were transfected with 1.6 µg/mL concentration of MSTRG.2316.1 vector or empty vector using Lipofectamine 3000 reagent (ThermoFisher, Carlsbad, CA, USA) according to the manufacturer's instructions. After eight hours, the cells were induced to differentiate (set day 0). The cells were harvested to detect OE efficiency at 48 h.

Quantitative Real-Time PCR (qRT-PCR)
Total RNA was extracted using Trizol reagent according to the manufacturer's instructions. The qRT-PCR primers used were designed using online Primer3 software (Table S1). Total RNA was reversed transcribed to complementary DNA (cDNA) using PrimeScripts RT Reagent Kit containing gDNA Eraser (TAKARA, Dalian, China). Then, the cDNA template used for qPCR was determined using the SYBR II master mix kit (TAKARA, Dalian, China). The qPCR was performed on a Bio-Rad CFX manager according to the manufacturer's instructions. The amplification reaction was conducted under the following program: pre-denaturation at 95 • C for 10 s, followed by 40 cycles of denaturation at 95 • C for 5 s and annealing/extension at 59 • C for 20 s. The melting curve analysis was performed Biology 2021, 10, 1176 5 of 15 from 65 to 95 • C with an increment of 0.5 • C. All the qRT-PCR Ct-values were normalized to the Ct-value of the ACTB gene and group D0 using the 2 −∆∆Ct method.

Statistical Analysis
Statistical analyses, including t-test and one-way analysis of variance (ANOVA), were conducted on R software. The p value < 0.05 was considered significant.

Histological Dynamics of ARBW in Rabbits
According to the histological analysis, the brownish color of iATs gradually faded during the development of rabbit iAT ( Figure 1A upper panel). There was an increase in cell diameter from D0 to Y2 and obvious heterogeneity of adipose tissue was found at D15 and D85. The cells in D0 were composed of multiple small triglyceride droplets (multilocular adipocytes), while cells in Y2 were composed of a single large lipid droplet (unilocular adipocytes) ( Figure 1A middle panel). The ratio of multilocular adipocytes to unilocular adipocytes gradually decreased from D0 to Y2 and the iATs in D85 contained a very low proportion of multilocular adipocytes ( Figure 1A middle panel). The IHC assay showed that the expression levels of UCP1 protein gradually declined during the development of iATs ( Figure 1A bottom panel). Furthermore, the results of the qRT-PCR assay showed that transcriptional copy numbers of the mitochondrial genes, including CYTB, COX2 and ND1, were dramatically decreased from D0 to D15 and then gradually decreased from D15 to Y2 in the tissue level ( Figure 1B). the cDNA template used for qPCR was determined using the SYBR II master mix kit (TAKARA, Dalian, China). The qPCR was performed on a Bio-Rad CFX manager according to the manufacturer's instructions. The amplification reaction was conducted under the following program: pre-denaturation at 95 °C for 10 s, followed by 40 cycles of denaturation at 95 °C for 5 s and annealing/extension at 59 °C for 20 s. The melting curve analysis was performed from 65 to 95 °C with an increment of 0.5 °C. All the qRT-PCR Ctvalues were normalized to the Ct-value of the ACTB gene and group D0 using the 2 −ΔΔCt method.

Statistical Analysis
Statistical analyses, including t-test and one-way analysis of variance (ANOVA), were conducted on R software. The p value < 0.05 was considered significant.

Histological Dynamics of ARBW in Rabbits
According to the histological analysis, the brownish color of iATs gradually faded during the development of rabbit iAT ( Figure 1A upper panel). There was an increase in cell diameter from D0 to Y2 and obvious heterogeneity of adipose tissue was found at D15 and D85. The cells in D0 were composed of multiple small triglyceride droplets (multilocular adipocytes), while cells in Y2 were composed of a single large lipid droplet (unilocular adipocytes) ( Figure 1A middle panel). The ratio of multilocular adipocytes to unilocular adipocytes gradually decreased from D0 to Y2 and the iATs in D85 contained a very low proportion of multilocular adipocytes ( Figure 1A middle panel). The IHC assay showed that the expression levels of UCP1 protein gradually declined during the development of iATs ( Figure 1A bottom panel). Furthermore, the results of the qRT-PCR assay showed that transcriptional copy numbers of the mitochondrial genes, including CYTB, COX2 and ND1, were dramatically decreased from D0 to D15 and then gradually decreased from D15 to Y2 in the tissue level ( Figure 1B). The expression levels of mitochondrial genes were detected at four growth stages during rabbit iAT whitening using qRT-PCR. The expression was normalized to the ACTB gene and D0. The data shows the means of three independent experiments. Two technical replicates were set for one individual experimental replicate. The "Rel." represents "Relative".

Identification and Characterization of lncRNAs in Rabbits iATs
To better define active lncRNAs during the ARBW of iATs, we reconstructed credible ted the transcriptome of iAT at the growth stages of D0, D15, D85 and Y2 in rabbits (n = 3 The expression levels of mitochondrial genes were detected at four growth stages during rabbit iAT whitening using qRT-PCR. The expression was normalized to the ACTB gene and D0. The data shows the means of three independent experiments. Two technical replicates were set for one individual experimental replicate. The "Rel." represents "Relative".

Identification and Characterization of lncRNAs in Rabbits iATs
To better define active lncRNAs during the ARBW of iATs, we reconstructed credible ted the transcriptome of iAT at the growth stages of D0, D15, D85 and Y2 in rabbits (n = 3 per stage). We performed paired-end ssRNA-seq for each tissue and obtained an average of 122.33 million clean reads from each sample, of which an average of 94.01% was properly mapped to the rabbit genome (Table S2). A total of 2281 novel lncRNAs were identified using different machine-learning models (Figure 2A). The novel lncR-NAs were then classified into 4 categories, including lincRNAs (1058), SO lncRNA (291), intronic lncRNAs (548) and AS lncRNAs (384) ( Figure 2B upper panel). By integrating the 1640 Esembl annotated lncRNAs (all were classified lincRNAs category), which were expressed in our samples, we obtained a total of 3921 lncRNAs expressed in the rabbit iATs ( Figure 2B bottom panel). Analysis of rabbit iAT lncRNA structure found that most lincRNAs, intronic lncRNAs and AS lncRNAs contained one or two exons, especially for intronic lncRNAs, while more SO lncRNAs contained multiple exons ( Figure 2C). When compared to protein-coding genes, lncRNAs were less expressed ( Figure 2D), which were in line with our previous lncRNA study in rabbits [20,21] and other animals [2]. Analysis of lncRNA transcriptional complexity of the iATs, visceral white adipose tissues (vWATs) and skeletal muscle tissues revealed that the complexities of the of iATs and skeletal muscle tissues were lower than those of vWATs. The iATs in D0, which represented the classical BAT, were the least complex tissue (top 10 expressed lncRNAs, accounting for approximately 80% of total lncRNA expression). Our data showed that the top 10 expressed lncRNAs in BATs were novel lncRNAs, such as MSTRG.11968.1 (an intronic lncRNA at ZNF777), MSTRG.17113.1 (an AS lncRNA of MAPK8) and MSTRG.14188.65 (a lincRNA located at an unplaced genome scaffold) ( Figure 2E). We characterized a catalog of credible novel lncRNAs expressed in the rabbit iATs and indicated that classical BATs were low lncRNA transcriptional complexity tissues.
We explored the expression correlation between lncRNAs and the protein-coding genes located at the corresponding lncRNA loci. For the AS lncRNAs, the expression patterns of 58 ones were significantly positively correlated with their anti-sense proteincoding genes across all samples (Pearson's correlation coefficient, p < 0.05). For instance, the AS lncRNAs located at the loci of FTX3, ENSOCUT00000058972, ID2, KCNA7, ENSO-CUT00000034065 were the top five lncRNA that had the highest correlation coefficients with their corresponding protein-coding genes ( Figure 2F). On the other hand, 10 AS lncRNAs were significantly negatively correlated with their anti-sense protein-coding genes, such as the AS lncRNAs located at the loci of Ubiquitin regulatory X (UBX) domain protein 2B (UBXN2B), sex determining region (SRY) transcription factor 5 (SOX5) and zinc finger and BTB domain containing 21 (ZBTB21) ( Figure 2G). For the intronic lncRNAs, the expression patterns of 91 and only 2 were significantly positively and negatively correlated with their corresponding protein-coding genes, respectively ( Figure S1A,B). For the SO lncRNAs, the expression patterns of 77 and only 2 were significantly positively and negatively correlated with their corresponding protein-coding genes, respectively ( Figure S1C,D). Alluvial diagram analysis of our AS lncRNAs, intronic lncRNAs and SO lncRNAs showed that approximately 25% of these lncRNAs were significantly correlated (p < 0.05) with their corresponding protein-coding genes. AS lncRNA that were negatively correlated with their corresponding protein-coding genes were prone to enrich in the negative strand of the genome ( Figure 2H). We allocated the lincRNAs to the flanking regions of annotated protein-coding genes and found that 2388 lincRNAs resided in the flanking regions of 7733 annotated protein-coding genes. The expression patterns of 731 lincRNA were significantly positively correlated with their flanking protein-coding genes and those of 92 lincRNA were significantly negatively correlated with their flanking protein-coding genes (Table S3). We identified that approximately 27% of total expressed lncRNAs were significantly correlated with their flanking (for lincRNAs) or reference (for intronic lncRNAs, SO lncRNAs and AS lncRNAs) genes. The Pearson's correlation coefficients and P values of significance tests were marked following protein-coding gene symbols. The "*" represents p < 0.05 and "**" represents p < 0.01. (H) AS lncRNA that were negatively correlated with their corresponding protein-coding genes were prone to enrich in the negative strand of genome. The "pos" represents "a positive correlation with protein-coding genes" and "neg" represents "expression of lncRNAs negatively correlated with protein-coding genes" and n.s. represents "a negative correlation with protein-coding genes". Most of the lncRNAs in "neg" came from AS lncRNAs that transcribed from negative strand (red stream).  The left heatmap depicts the lncRNA expression pattern and the right heatmap depicts the anti-sense protein-coding gene expression patterns across samples. The Pearson's correlation coefficients and p values of significance tests were marked following protein-coding gene symbols. The "*" represents p < 0.05 and "**" represents p < 0.01. (H) AS lncRNA that were negatively correlated with their corresponding proteincoding genes were prone to enrich in the negative strand of genome. The "pos" represents "a positive correlation with protein-coding genes" and "neg" represents "expression of lncRNAs negatively correlated with protein-coding genes" and n.s. represents "a negative correlation with protein-coding genes". Most of the lncRNAs in "neg" came from AS lncRNAs that transcribed from negative strand (red stream).  Figure 3C and Table S4). Therefore, the results of reconstruction of transcriptome and quantification of lncRNA expression were reliable.

Dynamics of lncRNA Expression during ARBW of Rabbit iATs
significantly enriched pathways of DELs in BATR3 and BATR4, respectively. The DELs in BATR1, BATR2 and BATR8 were upregulated from D0 to D15 and downregulated in D85. KEGG pathway enrichment showed that DELs in BATR1, BATR2 and BATR8 were significantly enriched in the white adipose development-related pathways, such as insulin resistance, glucagon signaling pathway and peroxisome proliferator-activated receptor (PPAR) signaling pathway. The DELs in BATR5 and BATR6 were upregulated from D0 to D85. KEGG pathway enrichment showed that the cyclic guanosine monophosphate (cGMP)/cGMP-dependent protein kinase (cGMP-PKG) signaling pathway and sphingolipid metabolism were the most significantly enriched pathways of DELs in BATR5 and BATR6, respectively. cGMP-PKG was significantly enriched by the DELs in BATR6. The DELs in BATR7 were constantly expressed from D0 to D85 but dramatically upregulated in Y2. The KEGG pathway enrichment showed that the top3 significantly enriched pathways of DELs in BATR7 were adrenergic signaling in cardiomyocytes, Wnt signaling pathway and lipid and atherosclerosis ( Figure 3D). The fatty acid metabolism and PPAR signaling pathway were also enriched by lncRNAs in BATR7. K-means clustering approach based on the TPM was used to sort all DELs, resulting in eight clusters, namely, BATR1 to BATR8, respectively ( Figure 3D). The DELs in BATR3 and BATR4 were expressed in D0 and downregulated in D15. The KEGG pathway enrichment showed that purine metabolism and Wnt signaling pathway were the most significantly enriched pathways of DELs in BATR3 and BATR4, respectively. The DELs in BATR1, BATR2 and BATR8 were upregulated from D0 to D15 and downregulated in D85. KEGG pathway enrichment showed that DELs in BATR1, BATR2 and BATR8 were significantly enriched in the white adipose development-related pathways, such as insulin resistance, glucagon signaling pathway and peroxisome proliferator-activated receptor (PPAR) signal-Biology 2021, 10, 1176 9 of 15 ing pathway. The DELs in BATR5 and BATR6 were upregulated from D0 to D85. KEGG pathway enrichment showed that the cyclic guanosine monophosphate (cGMP)/cGMPdependent protein kinase (cGMP-PKG) signaling pathway and sphingolipid metabolism were the most significantly enriched pathways of DELs in BATR5 and BATR6, respectively. cGMP-PKG was significantly enriched by the DELs in BATR6. The DELs in BATR7 were constantly expressed from D0 to D85 but dramatically upregulated in Y2. The KEGG pathway enrichment showed that the top3 significantly enriched pathways of DELs in BATR7 were adrenergic signaling in cardiomyocytes, Wnt signaling pathway and lipid and atherosclerosis ( Figure 3D). The fatty acid metabolism and PPAR signaling pathway were also enriched by lncRNAs in BATR7.

Selection of lncRNA Candidates and Functional Validation of lncRNA MSTRG.2316.1
To focus on the efforts of lncRNA function validation, we ranked candidate lncRNAs by their abundance, differential expression during ARBW and significant correlation with their flanking or reference genes. A total of 10 lncRNA candidates were selected (Table S5). When searching the 10 lncRNAs to NONCODE database using BLAST, we found that 7 out of the 10 lncRNAs had significant sequence hits (Evalue < 1 × 10 −6 ) in human and mouse lncRNAs, which might indicate the sequence conservation of lncRNAs between rabbits, mice and humans (Table S6). The histological analysis in this study ( Figure 1A) and the BAT master marker UCP1 read coverage ( Figure 4A) during the ARBW indicated that the BATs existed from D0 to D85. Only the samples in Y2 completed the whitening process, which spurred our interest in the cluster BATR7 and contained DELs expressed from D0 to D85, but dramatically upregulated in Y2 ( Figure 3D). One lncRNA candidate, MSTRG.2316.1, in BATR7 with the highest TPM value in our lncRNA candidates was validated and its expression changed little from D0 to D85 and markedly upregulated 35 folds in Y2 in the qRT-PCR validation ( Figure 4B). The 1449 bp length lincRNA MSTRG.2316.1 was located in chromosome 12, contained 2 exons on the positive strand of DNA and its read coverage dramatically increased in Y2 ( Figure 4C). MSTRG.2316.1 was upregulated from D85 to Y2 during the ARBW of iATs and visceral white adipocyte differentiation, suggesting a potential negative regulator of brown adipocyte development ( Figure 4D).
To validate the MSTRG.2316.1 functions, we first established the BPA differentiation model. Our results showed that BPAs isolated from iAT of rabbits in D0 were fusiform or triangular ( Figure 4E upper panel). The Oil red O staining of the induced mature BATs showed that many lipid droplets had accumulated after differentiation for 5 days ( Figure 4E bottom panel). The expression levels of common adipose markers, such as peroxisome proliferator activated receptor gamma (PPARG), CCAAT enhancer binding protein alpha (CEBPA), adiponectin C1Q and collagen domain containing (ADIPOQ) and fatty acid binding protein 4 (FABP4) and BAT-specific markers of cell death inducing DFFA like effector a (CIDEA), ELOVL fatty acid elongase 6 (ELOVL6), PPARG coactivator 1 alpha (PGC1A), peroxisome proliferator activated receptor alpha (PPARA) and UCP1 were significantly upregulated during differentiation of rabbit BPAs ( Figure 4F). The expression levels of MSTRG.2316.1 were significantly downregulated after the induced differentiation of rabbit BPAs ( Figure 4F).
To validate the function of MSTRG.2316.1, we performed gain-of-function analysis using plasmid vector-mediated overexpression (OE) of MSTRG.2316.1 in cultured BPAs. OE led to a significant 6 times increase of MSTRG.2316.1 level, compared to the cells treated with empty vector ( Figure 4G). Oil red O results showed that OE of MSTRG.2316.1 impaired brown adipose differentiation, compared to transfection of empty plasmid vector ( Figure 4H). Quantification of Oil red O of cells indicated that OE of MSTRG.2316.1 significantly decreased the lipid accumulation (p < 0.01 and Figure 4I). Additionally, the qRT-PCR analysis showed that OE of MSTRG.2316.1 significantly decreased FABP4, UCP1, CIDEA, ELOVL6, PCG1A, AIPOQ and PPARA expression (p < 0.05, Figure 4J). Our data showed that MSTRG.2316.1 was a negative regulator of BPA differentiation in vitro.  Figure 4B,F,G,J show the means of three independent experiments. Two technical replicates were set for one individual experimental replicate. The "*" represents p < 0.05 and "**" represents p < 0.01.

Discussion
The BAT histological characteristics had been depicted in some species, such as humans, mice and goats [6,11,40]. Our histological assays of rabbits suggested that the characteristics of rabbit BATs were in line with those of species showing the characteristics, The qRT-PCR data in (B,F,G,J) show the means of three independent experiments. Two technical replicates were set for one individual experimental replicate. The "*" represents p < 0.05 and "**" represents p < 0.01.

Discussion
The BAT histological characteristics had been depicted in some species, such as humans, mice and goats [6,11,40]. Our histological assays of rabbits suggested that the characteristics of rabbit BATs were in line with those of species showing the characteristics, including containing multiple small triglyceride droplets, the smaller size of cells and the highly expressed UCP1 and mitochondrial genes. BAT of most mammals undergoes ARBW. Previous studies have revealed the difference in the speed of ARBW among different species, such as hibernating animals persisting BAT into an adult, humans persisting race amounts of BAT into an adult and ruminants persisting BAT in 30 days after birth [6,10,40,41]. Our data revealed that the puberty rabbit contained a low proportion of multilocular adipocytes, suggesting that the ARBW of rabbits was slower than that of hibernating animals and ruminants and similar to that of humans. The long-time persisting BATs may explain the less fat-deposition in rabbits.
We identified 2281 credible novel lncRNAs by reconstructing the transcriptome, indicating the importance of lncRNAs in regulating the ARBW of rabbits. Intronic lncRNAs were transcribed from a single intron of protein-coding genes. Our data showed that intronic lncRNAs contained a few exons, which suggested that the exon number of lncRNAs might be affected by genomic elements of annotated protein-coding genes. Transcriptional complexity, expressed as the fraction of total RNAs, accounted for 10 or 100 most frequently expressed genes [42]. Our analysis revealed that classical BAT (iAT in D0) was a low lncRNA transcriptional complex tissue. On the other hand, the lncRNA transcriptional complexity between iATs and skeletal muscle tissues was similar and vWATs had a higher lncRNA transcriptional complexity than iATs and skeletal muscle tissues, which might suggest the transcriptomic difference between energy expenditure-and storage-related tissue.
Except for lincRNA, the AS lncRNA, intronic lncRNA and SO lncRNA were all transcribed from the loci of annotated protein-coding genes. Previous studies indicated that lncRNA transcribed from the loci of protein-coding genes might regulate the corresponding protein-coding genes [43,44]. For lincRNAs, we predicted their flanking protein-coding genes. By analyzing the expressional relationship between these three types of lncRNA and their corresponding protein-coding genes, we provided a valuable resource of the lncRNAs that were significantly correlated with their corresponding protein-coding genes, which can be used to identify functional lncRNAs efficiently.
In this study, analyses of lncRNA dynamics during rabbit iAT whitening were conducted. By comparing different stages during rabbit iAT development, we detected 631 DELs in the paired comparisons, which demonstrated that lncRNAs were widely involved in the BAT whitening in rabbits. We classified DELs into different clusters using the k-means method. The BAT is specialized for energy expenditure and heat generation, which depends on the function of mitochondria [45]. Purine nucleotides proved a constitutive inhibitor of UCP1 mediated proton conductance of the mitochondrial inner membrane and constituted the default shut-off mechanism in the absence of thermogenic demand [46]. Wnt signaling inhibits brown adipogenesis [47]. Our KEGG enrichment showed that the most significant pathway enriched by lncRNAs in BATR3 and BATR4 was purine metabolism and Wnt signaling pathway, which indicated that rapidly downregulated lncRNAs from D0 to D15 might mediate the ARBW through regulating their cis-regulated genes that were related to purine metabolism. A previous study indicated that a high-fat diet-induced BAT whitening and insulin resistance in mice [48]. In our study, the most significant KEGG pathway enriched by lncRNAs in BATR1, BATR2 and BATR8 were white adipose development-related pathways, such as insulin resistance and PPAR signaling pathway, suggesting the similarity pathway between obesity-related BAT whitening and ARBW and the lncRNA involvement. The cGMP-PKG pathway can promote the browning process of WATs. Our data indicated that the constantly upregulated lncRNAs from D0 to D85 were involved in this pathway, suggesting that these upregulated lncRNAs might play inhibitory roles during rabbit ARBW [49]. The KEGG enrichment showed that lncRNAs in BATR7 were widely involved in lipid metabolism pathways, such as lipid and atherosclerosis, fatty acid metabolism and PPAR signaling pathway, which indicated that this lncRNA might play crucial roles in white adipocyte development in the iATs of aged rabbits.
We characterized lncRNA MSTRG.2316.1 as a potential inhibitor that was upregulated in the iATs of aged rabbits and during differentiation of WATs. During the rabbit BPA differentiation, the common adipose markers and BAT-selective makers were all upregulated [50,51], which was similar to the previous mice BPA differentiation models. We established a robust rabbit BPA differentiation model. OE of MSTRG.2316.1 dramatically impaired the lipid accumulation in mature brown adipocytes and decreased the common adipose marker of FABP4 and the thermogenic marker genes (UCP1, CIDEA, ELOVL6, PGC1A and PPARA), indicating that MSTRG.2316.1 acted as an inhibitor during BAT development and might mediate the ARBW in aged rabbits. The molecular mechanisms underlying MSTRG.2316.1 regulating ARBW still need further investigation.

Conclusions
In summary, we provided comprehensive histological dynamics and detected a total of 2281 novel lncRNAs during ARBW of rabbit iATs. We revealed that classical rabbit BATs were low lncRNA transcriptional complex tissues. Dynamic analyses for lncRNAs identified 631 DELs during the ARBW. The purine metabolism pathways, Wnt signaling pathway, PPAR signaling pathway, cGMP-PKG signaling pathway and lipid and atherosclerosis were significantly enriched by the potential cis-regulated targets of the DELs with unique expression patterns during the ARBW. The novel identified lncRNA MSTRG.2316.1 can play an inhibitory role during brown adipocyte differentiation. Our work provides evidence that lncRNAs were widely involved in ARBW in rabbits, facilitating a better understanding of adipose biology.