LncRNA and mRNA expression profiles in brown adipose tissue of obesity-prone and obesity-resistant mice.

Obesity-prone or obesity-resistant phenotypes can exist in individuals who consume the same diet type. Brown adipose tissue functions to dissipate energy in response to cold exposure or overfeeding. Long noncoding RNAs play important roles in a wide range of biological processes. However, systematic examination of lncRNAs in phenotypically divergent mice has not yet been reported. Here, the lncRNA expression profiles in BAT of HFD-induced C57BL/6J mice were investigated by high-throughput RNA sequencing. Genes that play roles in thermogenesis and related pathways were identified. We found lncRNA (Gm44502) may play a thermogenic role in obesity resistance by interacting with six mRNAs. Our results also indicated that seven differentially expressed lncRNAs (4930528G23Rik, Gm39490, Gm5627, Gm15551, Gm16083, Gm36860, Gm42002) may play roles in reducing heat production in obesity susceptibility by interacting with seven differentially expressed mRNAs. The screened lncRNAs may participate in the pathogenesis of weight regulation and provide insight into obesity therapy.

Introduction

Obesity is a worldwide public health problem, and the factors that cause it are a source of public concern. A high-fat diet (HFD) is associated with obesity, but individuals have different propensities to accumulate excess weight or to stay lean (Giles et al., 2016). Variability in obesity susceptibility is associated with many genes and epigenetic events that favor obesity or leanness (Giles et al., 2016), and therefore individuals that consume the same diet type, including HFD, may gain different amounts of weight. Moreover, rodents that had the same genetic background and were fed the same diet type gained different amounts of weight. Among them, obesity-prone (OP) or obesity-resistant (OR) phenotypes were identified (Choi et al., 2016). Transcriptomic profiles of white adipose and liver tissue in OP and OR mice and rats had common and divergent patterns of transcriptional changes (Choi et al., 2016; Li et al., 2008). However, there is currently no corresponding expression profile in brown adipose tissue in OP and OR mice.

Activation of brown adipose tissue (BAT) may reduce obesity and related diseases because active BAT has been shown to dissipate energy and contribute 2.5%-5% of the body’s resting metabolic rate, thereby promoting negative energy balance (Bostrom et al., 2012; Chiang et al., 2009; Seale et al., 2011). BAT is present in adults (Cypess et al., 2009; Van Marken Lichtenbelt et al., 2009; Virtanen et al., 2009), and high amounts of functional BAT have been associated with reduced obesity in humans (Saito et al., 2009). BAT generates energy in the form of heat in response to cold exposure or overfeeding (Rothwell and Stock, 1997; Cannon and Nedergaard, 2004). In rodents, overfeeding was found to activate BAT as a way to limit weight gain (Rothwell and Stock, 1979). Increased BAT activity or increased numbers of thermogenic adipocytes allowed mice to remain lean and healthy (Bostrom et al., 2012; Chiang et al., 2009; Seale et al., 2011). BAT is characterized by the unique expression of uncoupling protein-1 (UCP-1) in mitochondria. Upon activation, UCP-1 dissipates the electrochemical gradient for ATP synthesis, thereby generating heat (Klingenberg, 1999; Kozak and Harper, 2000). The proliferation and activity of thermogenic BAT may have potential as an obesity therapy (Bai et al., 2017). Some transcription factors have been detected in brown adipocytes, and finding other factors and pathways associated with brown adipocyte thermogenic function is of great interest.

Long non-coding RNAs (lncRNAs) are RNA transcripts that are >200-nt long. They are transcribed by RNA polymerase II and undergo splicing and polyadenylation similarly to mRNAs (Nagano and Fraser, 2011). However, lncRNAs lack functional open reading frames, and therefore their protein-coding potential is limited. Compared with mRNAs, lncRNAs tend to be expressed at lower levels (Bono et al., 2003; Ramskold et al., 2009; Babak et al., 2005), and most lncRNAs are preferentially expressed in specific tissues or cell lines (Babak et al., 2005; Mercer et al., 2008). The expression of adipose lncRNAs is low and depot-specific (Zhao et al., 2014; Divoux et al., 2014). Moreover, lncRNAs are located within the nuclear or cytosolic fractions and possess heterogeneous structures, which determine the diversity and complexity of lncRNA functions (Wei et al., 2016). LncRNAs can interact with DNA, RNA, or proteins and regulate the epigenetic state of chromatin, gene expression, tissue development, and tumorigenesis (Chen and Carmichael, 2010; Schonrock et al., 2012). To date, a few studies of lncRNAs in BAT have been reported. A previous study found that a lncRNA, brown fat lncRNA 1 (Blnc1), formed a ribonucleoprotein complex with the transcription factor early B cell factor 2 (EBF2), which regulated classic brown fat by regulating the thermogenic adipocyte differentiation (Zhao et al., 2014). Another lncRNA, LINC00473, was strongly induced by the activation of thermogenic adipocytes, acting not only through nuclear transcriptional regulation but also by shuttling into the cytoplasm to affect thermogenic adipocyte physiology (Tran et al., 2020). However, the repertoire and potentially functional characteristics of lncRNAs in BAT of individuals with OP and OR phenotypes remain poorly understood.

In this study, we systemically investigated and profiled lncRNA and mRNA expression levels in BAT of HFD-induced OP and OR mice by high-throughput sequencing. Bioinformatics analysis was used to find lncRNAs with potential thermogenic functions, and in vivo and in vitro experiments were performed to verify the importance of the candidate genes in BAT. These lncRNAs will provide preliminary insights into the pathogenesis of weight regulation and obesity.

Results

High-fat diet-induced obesity-prone and obesity-resistant mouse models

To explore the transcriptomic profile of the differences in weight gain after consuming an HFD, we performed an RNA sequencing (RNA-seq) analysis to identify lncRNAs involved in reducing the extent of weight gain in BAT. We established a mouse obesity model by feeding mice an HFD for 22 weeks beginning at 4 weeks of age and then divided the mice into two groups (OP and OR) according to their weight gain; a third group was fed a regular diet as a control (normal diet group, ND) (Figure 1A). After 22 weeks, we isolated the BAT, extracted and purified the total RNA, and performed RNA-seq to detect lncRNAs and mRNAs associated with weight and metabolism regulation. The expression level of the uncoupling protein 1 gene (Ucp1), which is involved in thermogenesis regulation and energy metabolism in BAT, was measured by real-time polymerase chain reaction (qPCR). The results showed that Ucp1 expression was significantly higher in the BAT from OR and ND mice than it was in the BAT from OP mice (Figure 1B). A significant increase in lipid accumulation was detected in the liver of OP mice compared with that in the liver of OR and ND mice by hematoxylin and eosin staining (Figure 1C).

Figure 1

Characteristics of high-fat diet (HFD)-induced obesity-prone (OP) and obesity-resistant (OR) mice

(A) Weight changes in mice fed a normal diet (ND) or an HFD. OP, obesity-prone mice; OR, obesity-resistant mice; ND, control. Values represent the mean ± SD. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001 (unpaired, two-tailed t test, n ≥ 3).

(B) Ucp1 expression levels in brown adipose tissue of mice fed an ND or an HFD. Error bars represent the mean ± SD. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001 (unpaired, two-tailed t test, n ≥ 3).

(C) Representative images of hematoxylin and eosin-stained liver of the mice fed an ND or an HFD, 200× magnification level, scale bar = 100μm.

Expression profiles and gene set enrichment analysis of the transcriptomes from the obesity-prone and obesity-resistant groups

RNA-seq analysis was performed to obtain the expression profiles of lncRNAs and mRNAs associated with weight regulation in the BAT from HFD-fed OP and OR mice. We obtained a total of 31,466 genes; 18,360 and 12,767 were annotated as mRNAs and lncRNAs, respectively (Figures 2A-2C).

Figure 2

Expression profiles and gene set enrichment analysis of the transcriptomes of mice in the obesity-prone (OP) and obesity-resistant (OR) groups

(A) Venn diagram of the numbers of genes in the OP and OR groups. The two groups have 27,710 shared genes, with 2110 and 1646 unique genes in the OP and OR groups, respectively.

(B and C) Venn diagrams of the numbers of mRNAs (B) and lncRNAs (C) in the OP and OR groups. The two groups have 17,259 shared mRNAs, with 626 and 475 unique mRNAs in the OP and OR groups, respectively, and 10,158 shared lncRNAs, with 1457 and 1152 unique lncRNAs in the OP and OR groups, respectively.

(D) Top pathway terms ranked 3: thermogenesis. A total of 243 genes were enriched in thermogenesis, and 118 of them were in the leading-edge subsets.

(E) Expression of lncRNA (Gm44502) in thermogenesis by RNA-seq. Gm44502 expression was higher in OR mice than it was in OP mice. Values represent the mean ± SD. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001 (unpaired, two-tailed t test, n = 3).

(F–H) Thermogenesis-related pathway terms: Oxidative phosphorylation (F), MAPK signaling pathway (G), and cAMP signaling pathway (H). NES, normalized enrichment score; FDR, false discovery rate.

To interpret the RNA-seq data, we performed a Gene Set Enrichment Analysis (GSEA) using Kyoto Encyclopedia of Genes and Genomes (KEGG) gene lists and pathways to determine the potential functions and interactions of all the genes in the BAT. We used stringent detection criteria (i.e., normalized enrichment score |NES| >1.0, nominal p-value <5%, false detection rate (FDR) < 25%) and found that the OR group was enriched for 26 of the 74 detected gene sets, compared to the OP group, including the Oxidative phosphorylation, Thermogenesis, Citrate cycle (TCA cycle), Cardiac muscle contraction, Retrograde endocannabinoid signaling, and Glucagon signaling pathways. Detailed information for the 26 significant enriched pathway terms is listed in Table S1. Thermogenic pathways play important roles in BAT and therefore it was not unexpected to find that thermogenesis was one of the top three significant pathways (NES = 2.20, p = 0, FDR = 0.0002) enriched using the strict criteria (Figure 2D). 243 genes were enriched in thermogenesis as candidate genes of interest, and 118 of them were in the leading-edge subsets, including a lncRNA (Gm44502) (Table S2). Genes in the leading-edge subsets are meant to play a dominant role in this pathway. The expression level of the lncRNA (Gm44502) was higher in the OR group than it was in the OP group (Figure 2E), and its characteristics were detailed in Table S3, Figures S1, andS2. The current version of lncLocator predicted Gm44502 was localized to the cytosol (85%, Figure 6E and Table S15). Besides, the oxidative phosphorylation pathway, one of the twenty-six gene sets, is associated with thermogenesis (Figures 2F and S3). The OP group was enriched for 83 of 173 detected gene sets, compared to the OR group, including osteoclast differentiation, cytokine-cytokine receptor interaction, hematopoietic cell lineage, natural killer cell mediated cytotoxicity, toll-like receptor signaling, C-type lectin receptor signaling, B cell receptor signaling, and phagosome pathways (Table S4). Both the MAPK signaling pathway and the cAMP signaling pathway are associated with thermogenesis (Figures 2G and 2H, and S3).

Figure 6

LncRNA-mRNA co-expression and functional networks

(A) Module-trait relationships. Each row represents a module, the columns are the traits (obesity-prone or obesity-resistant). Each cell contains the corresponding correlation coefficient with the p-value. The numbers of mRNAs and lncRNAs in the turquoise and yellow modules are also shown.

(B and C) Bubble patterns of KEGG pathway enrichment analysis represent 10 and 13 significant KEGG pathway terms in the turquoise (B) and yellow (C) modules, respectively. Pathway terms in red are related to thermogenesis.

(D) Graphic views of lncRNA-mRNAs co-expression network of the turquoise module.

(E) Graphic views of cytoplasmic-nuclear localization prediction of 8 screened lncRNAs.

(F) Co-expression network of hub lncRNAs and mRNAs (E). Node size indicates corresponding intramodular connectivity degrees; diamonds and circles represent lncRNAs and mRNAs, respectively. Connecting lines between nodes represent direct correlations between genes.

Protein-protein interaction network analysis and functional cluster analysis

Protein-protein interaction (PPI) network analysis was performed to find key genes enriched in thermogenesis, investigate related pathways in heat production, and obtain the potential regulatory pathways involving the lncRNA (Gm44502). Genes with a combined score >950 for interactions (scores ranged from 0 to 1000) were considered significant and possible core genes. The PPI network analysis identified 141 mRNAs from the 242 mRNAs enriched in thermogenesis (Figure 3A and Table S5). This network consists of known interactions from curated databases and interactions that were experimentally determined. The direct or indirect functional protein interactions were visualized. Four functional clusters were identified by clustering analysis of the whole network (Figure 3A). Pathway enrichment analysis of the genes in the four clusters showed that the clusters were involved not only mainly in thermogenesis but also in other thermogenesis-related pathways, namely the Oxidative phosphorylation, AMPK signaling, cAMP signaling, cGMP-PKG signaling, MAPK signaling, and mTOR signaling pathways (Figure 3B-3E and Table S6).

Figure 3

Protein-protein interaction (PPI) network construction and functional cluster analysis

(A) The whole PPI network with four clusters. The color and size of the circles indicate the degree of difference in mRNA expression between OP and OR mice.

(B-E) Bar graphs of the KEGG pathway enrichment analysis of cluster 1 (B), cluster 2 (C), cluster 3 (D), and cluster 4 (E). Pathway terms in red are related to heat production.

Co-expression analysis of screened mRNAs and lncRNAs and functional prediction

To further explore the interaction of lncRNAs in thermogenesis, we determined the correlation of the enriched lncRNA (Gm44502) and 141 candidate mRNAs from the four clusters by calculating the Pearson correlation coefficient (r) (Table S7). A strong correlation (|r| >0.9) was found between lncRNA (Gm44502) and six of the mRNAs (Adcy6, Adcy8, Ppargc1a, Sos1, Atp5g1, Smarca2) (Figure 4A and Table S8). The expression levels of four of these mRNAs (Adcy8, Ppargc1a, Adcy6, Sos1) were upregulated, and two (Smarca2, Atp5g1) were downregulated in the OR group compared with their expression levels in the OP group (Figure 4B). Moreover, in thermogenesis, Adcy6, Adcy8, and Ppargc1a were included in the leading-edge subsets by the GSEA above (Table S2).

Figure 4

Co-expression and expression of the identified mRNAs

(A) Pearson correlation between Gm44502 and six mRNAs. LncRNA (Gm44502) was strongly correlated with the six mRNAs (|r| >0.9). Red and blue indicate positive and negative correlations with Gm44502, respectively.

(B) Expression levels of six mRNAs in thermogenesis by RNA-seq. The expression of four mRNAs was increased and that of two mRNAs was decreased in the OR group compared with their expression in the OP group. Values represent the mean ± SD. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001 (unpaired, two-tailed t test, n = 3).

Expression profiles of differentially expressed lncRNAs and mRNAs and Kyoto Encyclopedia of Genes and Genomes analysis

To further characterize the differentially expressed genes (DEGs) from our RNA-seq data, we calculated the significance of 31,466 identified genes between the OP and OR groups using the DESeq2 R package (|log2fold change| ≥1, FDR ≤0.001 as the threshold) (Love et al., 2014). We have identified 1387 genes with significantly changed expression levels between the OR and OP groups (Table S9). The scatterplot in Figure 5A shows the genes that were significantly differentially expressed between the OR and OP groups. In the OP group, 123 DElncRNAs and 911 DEmRNAs were upregulated, and 105 DElncRNAs and 248 DEmRNAs were downregulated compared with their expression in the OR group (Figure S4). A530050N04Rik was the DElncRNAs with the highest expression and it was upregulated in the OR group compared with the OP group (Figure 5B). We clustered the top 80 differentially expressed lncRNAs (DElncRNAs) and visualized them as expression heatmaps by normalized raw z-scores (Figure 5C). We also clustered the top 100 differentially expressed mRNAs (DEmRNAs) (Figure 5D).

Figure 5

Expression profiles of differentially expressed genes (DEGs) and their functional analysis

Scatter plot of the DEGs between the obesity-prone (OP) and obesity-resistant (OR) groups. Red and green indicate upregulated and downregulated genes, respectively.

(A and B) Expression of DElncRNA (A530050N04Rik) between the OP and OR groups by RNA-seq. A530050N04Rik expression was higher in OR mice than it was in OP mice. Values represent the mean ± SD. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001 (unpaired, two-tailed t test, n = 3).

(C and D) Hierarchical clustering shows the top 80 distinguishable upregulated and downregulated DElncRNAs (C) and the top 100 upregulated and downregulated DEmRNAs (D) by expression profiling among samples and normalized raw z-scores. The squares on the left indicate gene-enriched level 2 KEGG pathway terms. Red indicates high relative expression; green indicates low relative expression.

(E and F) KEGG pathway enrichment network analysis shows 24 and 4 significant KEGG pathways that were enriched in the upregulated (E) and downregulated (F) genes. Red and blue indicate upregulated and downregulated KEGG pathways, respectively. The line chart shows the candidate gene number for the KEGG pathway terms. Pathway terms in red are related to thermogenesis.

To obtain the potential functions of DEGs in BAT, we performed a KEGG pathway enrichment analysis on 1159 DEmRNAs. We found that the 911 upregulated DEmRNAs in the OP group were significantly enriched in 24 out of 210 pathways (Q-value ≤0.05). The 24 KEGG pathways were divided into three level 1 KEGG categories, namely Organismal Systems, Cellular Processes, and Environmental Information Processing. The corresponding level 2 KEGG level 2 terms were shown in Figure 5E, including Hematopoietic cell lineage, Phagosome, Osteoclast differentiation, Cytokine-cytokine receptor interaction, B cell receptor signaling pathway, Chemokine signaling pathway, and ECM-receptor interaction. Among the 24 pathways, the cAMP signaling pathway is associated with thermogenesis (Figures 5E and S3). Detailed information for the 24 significant enriched pathways was provided in Table S10. Our above GSEA analysis showed these pathways were dysregulated in the obesity-susceptible phenotype. Four out of the 137 KEGG pathways were significantly downregulated in the OP group (Q-value ≤0.05) (Figure 5F and Table S11). These four KEGG pathways were divided into two level 1 KEGG categories: Organismal Systems and Environmental Information Processing. The corresponding level 2 KEGG terms (Figure 5F) included the PPAR signaling pathway, which is associated with thermogenesis.

LncRNA-mRNAs co-expression network

To further explore genes involved in important interactions associated with the OP and OR phenotypes, we built co-expression networks for the 1387 DEGs (228 DElncRNAs, and 1159 DEmRNAs, Table S9) by weighted gene co-expression network analysis (WGCNA), followed by network modules testing. We found four co-expressed modules (Figures 6A and Table S12), two of which were highly correlated with the OP and OR phenotypes, namely the turquoise and yellow modules (correlation p value <0.05, Pearson r ≥ 0.95). To determine the potential functions of the enriched genes in each module, we performed enrichment analysis for the DEmRNAs co-expressed with the DElncRNAs. The turquoise module including 45 DElncRNAs and 370 DEmRNAs was enriched with ECM-receptor interaction, Protein digestion and absorption, cAMP signaling pathway, Phagosome, Platelet activation, Hematopoietic cell lineage, Toll-like receptor signaling pathway, PI3K-Akt signaling pathway, and Focal adhesion (Figure 6B and Table S13). It was consistent with identifying dysregulated pathways in the obesity-prone phenotype by GSEA. In the above KEGG analysis, it was also found that the heat-related cAMP signaling pathway was significantly up-regulated in OP mice. Genes in the yellow module associated with the obesity-resistant phenotype were significantly enriched (Q-value ≤0.05) in 13 pathways (Figure 6C and Table S14). Among the 13 pathways, the PPAR signaling, cGMP-PKG signaling, and MAPK signaling pathways are important pathways associated with heat production. In the above KEGG analysis, the PPAR signaling pathway was found to be down-regulated in OP mice.

Next, to screen DEGs that were closely related in the turquoise, the threshold of weighted edges was set as 0.7 in the co-expression network (Figure 6D). Seven hub DElncRNAs out of the 28 DElncRNAs were Gm15551, Gm16083, Gm36860, Gm39490, Gm42002, Gm5627, and 4930528G23Rik. The characteristics of the lncRNAs were detailed in Table S3. Six of these DElncRNAs were predicted to be assigned to the cytoplasm, while the lncRNA (Gm39490) was predicted to reside in the nucleus (58%) (Figure 6E and Table S15). Seven DEmRNAs (C1qb, Lep, Aoc3, Fcer1g, Cmklr1, Ccl5, Dse) were closely associated with these seven DElncRNAs in the cluster co-expression network (Figure 6F).

Validation of important candidate genes by real-time polymerase chain reaction

To validate the candidate genes screened by our RNA-seq analysis, we measured their expression in brown adipocytes activated with Forskolin or NE in vitro, or in the BAT from cold-induced thermogenic mice in vivo.

In the BAT from cold-induced thermogenic mice, the expression levels of the lncRNA (A530050N04Rik) and three mRNAs (Adcy6, Ppargc1a, Sos1) were considerably increased, while that of the mRNA (Smarca2) was significantly decreased (Figure 7A). Besides, the expression levels of lncRNA (Gm44502) and mRNA (Adcy8) in the BAT from cold-induced thermogenic mice were up, and that of mRNA (Atp5g1) was down (Figure 7A). Furthermore, A530050N04Rik and Ppargc1a were significantly up-regulated, while Smarca2, Atp5g1, Adcy6, and Gm44502 were down-regulated in brown adipocytes activated with Forskolin or NE in vitro (Figures 7B and S5). Likewise, Sos1 and Adcy8 did not change noticeably (Figure S5).

Figure 7

Validation of selected genes by real-time PCR (qPCR)

(A and B) Validation of the expression of screened lncRNAs and mRNAs in brown adipose tissue from mice subjected to cold-induced thermogenesis in vivo by qPCR. The expression of two genes was decreased and six were increased in cold-stimulated mice. Values represent the mean ± SD. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. ∗∗∗∗p < 0.0001 (unpaired, two-tailed t test, n = 5) Further validation of the screened genes in activated brown adipocytes. Values represent the mean ± SD ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. ∗∗∗∗p < 0.0001 (one-way ANOVA, n = 3).

(C and D) Validation of the expression of screened lncRNAs (C) and mRNAs (D) in activated brown adipocytes by qPCR. Values represent the mean ± SD. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. ∗∗∗∗p < 0.0001 (one-way ANOVA, n = 3).

Next, the expression levels of seven selected DElncRNAs (4930528G23Rik, Gm39490, Gm5627, Gm15551, Gm16083, Gm36860, Gm42002) (Figure 7C) and seven selected DEmRNAs (C1qb, Lep, Aoc3, Fcer1g, Cmklr1, Ccl5, Dse), screened from thermogenesis-related pathways, were significantly down-regulated (Figure 7D) in the activated brown adipocytes.

Discussion

We conducted the first study of lncRNA profiles in BAT of obesity-prone (OP) and obesity-resistant (OR) mice via high-throughput RNA sequencing and functional predictive analysis. Previous studies of RNA expression in mice with high-fat diet-induced epigenetic differences were mainly in white adipose and liver tissue (Choi et al., 2016; Li et al., 2008); however, the repertoire and underlying functional characteristics of lncRNAs in BAT remain poorly understood. We paid attention to the uniqueness of BAT in thermogenesis and focused on the comprehensive expression, interaction network, and functional analysis of lncRNA and mRNA in thermogenesis and related pathways in BAT of mice with phenotypic differences. Moreover, we validated the expression changes in candidate genes in thermogenic function by in vivo or in vitro experiments.

First, we conducted GSEA and focused on the thermogenesis pathway, which was ranked the third most significantly enriched pathway based on the NES. Genes enriched in this gene set were closely associated with the thermogenic function in BAT. Among them, the expression level of lncRNA (Gm44502) (https://www.ncbi.nlm.nih.gov/nuccore/NR_004843.2), predicted to be in the cytosol, was higher in the OR group than it was in the OP group. A previous study has found that Gm44502 (Chkb-cpt1b) is significantly elevated in BAT compared to other adipose tissues (Rosell et al., 2014). In prediabetic and diabetic models, lncRNA (Gm44502) associated with Cpt1b was positively correlated with the lipolysis pathway (Kazeminasab et al., 2021). It suggests Gm44502 may have a potential thermogenesis regulatory role in BAT. Moreover, six mRNAs (Adcy6, Adcy8, Atp5g1, Ppargc1a, Sos1, Smarca2) were strongly correlated with the lncRNA (Gm44502) were key mRNAs in the PPI network analysis of thermogenesis. These six mRNAs were distributed in four clusters in the network and were also significantly enriched in heat-related pathways. We also found that lncRNA (Gm44502) and mRNAs (Ppargc1a, Adcy6, Adcy8) were from the leading-edge subsets in thermogenesis by GSEA, implying that they played a dominant role.

We then validated the expression changes in the lncRNA (Gm44502) and the 6 mRNAs (Adcy6, Adcy8, Ppargc1a, Sos1, Smarca2, Atp5g1) in the thermogenic function of BAT. In vivo, we found that the expression trends of lncRNA (Gm44502) and the six mRNAs in BAT of cold-induced thermogenic mice were consistent with those in the OR group by RNA-seq analysis. Ppargc1α (PGC-1α), Sos1, Adcy6, and Adcy8 were up-regulated in cold-induced mice by qPCR analysis. Previous studies suggested a dominant role for PGC-1α in its thermogenic function, brown fat development, and differentiation (Kajimura et al., 2010; Leone et al., 2005; Uldry et al., 2006). Grb2/Sos1 complexes were found to be involved in promoting adipogenic differentiation and inhibiting adipogenesis (Tao et al., 2016). And previous studies showed that Adcy6 (Salas-Perez et al., 2019) and Adcy8 (Sung et al., 2016) were strongly linked to obesity, and Adcy8 was required for hypothalamic adaptation to an HFD (Raoux et al., 2015). Meanwhile, we also found that Smarca2 and Atp5g1 were down-regulated in cold-induced mice by qPCR. Previous studies indicated that Atp5g1 affected obesity and metabolic diseases through the epigenetic regulation of gene expression (Johnson et al., 2014), and Smarca2 was involved in promoting fat development by contributing to the activation of adipogenic factors, mainly (Nguyen et al., 2015). Furthermore, in activated brown adipocytes, we found that the expression of Ppargc1a, Smarca2 was also increased, while that of Atp5g1 was decreased. It indicates that differential expression of the three genes in thermogenesis may be related to adrenergic. Thus, we speculate that Gm44502 may play a thermogenic role in obesity resistance by interacting with Adcy6, Adcy8, Ppargc1a, Sos1, Atp5g1, and Smarca2.

Second, to further characterize the DEGs from our RNA-seq data, we screened 228 DElncRNAs and 1159 DEmRNAs in BAT of HFD-induced mice. The most highly expressed DElncRNA (A530050N04Rik, also known as lncBATE10 (Bai et al., 2017)) was found in the OR group, up-regulated in the activated brown adipocyte, and induced by acute cold exposure. Previous studies have shown that A530050N04Rik is required for brown adipose differentiation and BAT thermogenesis (Bai et al., 2017; Salviano-Silva et al., 2018). The KEGG enrichment analysis of the DEGs showed the upregulated cAMP signaling and downregulated PPAR signaling pathways, which are related to heat production (Guo et al., 2019), in the OP group compared with the OR group. Besides, to further screen DElncRNAs and DEmRNAs that play key roles in phenotypic differences in mice, we performed WGCNA (correlation p-value <0.05, Pearson r ≥0.95). The yellow module, which was significantly associated with the OR phenotype, was significantly enriched for heat-related pathways such as the PPAR signaling, and MAPK signaling pathways. The 415 genes, including 45 DElncRNAs and 370 DEmRNAs, were found in the turquoise module associated with the OP phenotype. These genes were significantly enriched in thermogenesis-related pathways, including the cAMP signaling pathway that was up-regulated in OR mice. It suggests that these lncRNAs may function in reducing the extent of weight gain. The function of lncRNA is also closely related to their subcellular localization. In the nucleus, lncRNAs regulate gene expression at the epigenetic and transcriptional levels, and in the cytoplasm, they regulate gene expression at the post-transcriptional and translational levels (Zhao et al., 2020). Our results showed that the lncRNA (Gm39490) resided in the nucleus and the 6 lncRNAs (Gm15551, Gm16083, Gm42002, Gm39490, Gm5627, and 4930528G23Rik) were assigned to the cytoplasm. It indicated that 6 lncRNAs and Gm39490 may have different functions, among which 6 lncRNAs in the cytoplasm were more likely to act as a competitive endogenous RNA and interact with mRNAs to participate in expression regulation.

Moreover, we found that seven hub DElncRNAs were closely co-expressed with seven DEmRNAs by constructing the co-expression network. We then verified that the expression trends of these potentially interacting DElncRNAs (4930528G23Rik, Gm39490, Gm5627, Gm15551, Gm16083, Gm36860, Gm42002) and DEmRNAs (C1qb, Lep, Aoc3, Fcer1g, Cmklr1, Ccl5, Dse) in vitro were similar to those in OP group by the RNA-seq analysis. In the brown adipocytes activated with Forskolin or NE, the seven DElncRNAs and seven DEmRNAs were significantly decreased. A previous study suggested that the chemerin-Cmklr1 axis is a physiological negative regulator of thermogenic beige fat, and genetic blockade of adipocytic Cmklr1 was found to protect against diet-induced obesity (Lin et al., 2021). And Ccl5/CCR5 signaling, in part, suppressed energy expenditure and adaptive thermogenesis to exacerbate the development of obesity and insulin resistance (Chan et al., 2022).

C1qb (Gabrielsson et al., 2003), Aoc3 (Agha et al., 2015), and leptin (Machinal-Quelin et al., 2002) were suggested to be biologically relevant to the development of adiposity. Additionally, Fcer1g is positively associated with diet-induced obesity, and the knockout of the Fcer1g gene reduced the development of diet-induced obesity (Van Beek et al., 2015). Thus, we speculate that these DElncRNAs (4930528G23Rik, Gm39490, Gm5627, Gm15551, Gm16083, Gm36860, Gm42002) may play roles in thermogenesis-related pathways by interacting with these DEmRNAs (C1qb, Lep, Aoc3, Fcer1g, Cmklr1, Ccl5, Dse).

In conclusion, we performed the expression profiles and functional analysis of lncRNAs and mRNAs in BAT of HFD-induced OP and OR mice. We found that lncRNA (Gm44502) may play a thermogenic role in obesity resistance by interacting with Adcy8, Adcy6, Ppargc1a, Sos1, and Smarca2. Next, we obtained results consistent with previous studies that DElncRNA (A530050N04Rik) was highly expressed in BAT thermogenesis, supporting its role in thermogenesis. We also found seven DElncRNAs (4930528G23Rik, Gm39490, Gm5627, Gm15551, Gm16083, Gm36860, Gm42002) that may play roles in reducing heat production in obesity susceptibility by interacting with seven DEmRNAs (C1qb, Lep, Aoc3, Fcer1g, Cmklr1, Ccl5, Dse). All of these lncRNAs may be involved in regulating body weight changes by interacting with the mRNAs in BAT. The activation of BAT to expend energy by stimulating thermogenesis has been considered a treatment strategy for metabolic diseases related to obesity (Jung et al., 2019). Therefore, these lncRNAs were potential therapy and diagnosis target for obesity and obesity-related diseases.

Limitations of the study

However, further experimental studies are needed to investigate the functions of these lncRNAs in the OR and OP phenotypes among HFD-induced obesity individuals, as well as to determine the underlying mechanisms involved.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Shankuan Zhu (zsk@zju.edu.cn).

Materials availability

This study did not generate new unique reagents. All reagents used in this study will be made available on request to the lead contact.

Experimental model and subject details

Diet-induced obesity model

Male C57BL/6 mice were kept at the Zhejiang University Institutional Animal Care (Approval no. ZJU20170237). Mice were fed a regular chow diet (13% of calories from fat) or a high-fat diet (Research Diets, #D12492; 60% calories from fat) for 22 weeks beginning at 4 weeks of age. According to the weight gain after 22 weeks and the type of diet fed (De Sa et al., 2015), we divided the mice into three groups: mice fed the regular chow diet (ND, n = 4); mice fed the HFD that gained less weight than the other HFD group (obesity-resistant, OR, n = 3); and mice fed the HFD that gained more weight than the other HFD group (obesity-prone, OP, n = 3). The weight differed significantly among the groups. Bodyweight was measured every week. The BAT depots were rinsed with physiological saline, weighed, and immediately frozen in liquid nitrogen. RNA was extracted from the BAT depots and RNA sequencing was performed using the BGISEQ-500 High-throughput Sequencing System at the Beijing Genomics Institution (BGI).

Cold-induced thermogenesis

For cold-induced thermogenesis, C57BL/6J mice (Shanghai SLAC Laboratory Animal Co., Ltd) fed a regular diet for 8 weeks were moved to 4 °C for 3 h. The mice were sacrificed and the BAT was isolated for further analysis. The experimental protocols were approved by the Institutional Animal Care and Use Committee of Zhejiang University (ZJU20170237).

Adipocyte culture and differentiation

Immortal preadipocyte lines from the BAT of mice (donated by Prof. Zhuoxian Meng (Wang et al., 2014)) were cultured in high-glucose Dulbecco’s modified Eagle medium (DMEM, Gibco) containing 10% fetal bovine serum (Gibco). To induce adipogenic differentiation of these lines to brown fat cells, after reaching confluence, the preadipocytes were cultured in an induction medium containing 20 nM insulin, 1 μM dexamethasone, 0.5 mM 3-isobutyl-1-methylxanthine, 1 nM triiodothyronine (T3), 125 μM indomethacin, and 10% fetal bovine serum (FBS) for 2 days and then in differentiation medium containing 20 nM insulin, 1 nM triiodothyronine (T3), and 10% FBS for an additional 2 days. After that, the differentiation medium was replaced with a fresh differentiation medium every 2 days until day 7. To activate UCP1, after starvation culture in a serum-free medium for 6h, 1 μM norepinephrine bitartrate monohydrate (NE) or 10 μM Forskolin was added to the induction medium for 6 h.

Method details

Hematoxylin and eosin staining

The liver was fixed in 4% paraformaldehyde (pH 7.4) for 24 h and embedded in paraffin. The livers were sliced and stained with hematoxylin and eosin. All sections were observed with a microscope (Olympus).

High-throughput sequencing

Total RNA was extracted from the BAT using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The total RNA from each sample was treated with Ribo-Zero™ Magnetic Kit (Epocentre) to deplete any ribosomal RNA. The retrieved RNA was fragmented by adding First Strand Master Mix (Invitrogen). First-strand cDNA was generated using random primers reverse transcription, followed by second-strand cDNA synthesis in which the RNA template was removed and dTTP was replaced by dUTP. The incorporation of dUTP prevents the second strand of cDNA from being amplified in the subsequent process because the polymerase cannot cross the dUTP site on the template during extension. The double-stranded cDNA product was then ligated with an "A" base and a linker. The ligation product was amplified, and the final 10 cDNA sequencing libraries were obtained after purification (three biological replicates of the OP and OR groups, and four biological replicates of ND). The final library was qualified and quantitated using an Agilent 2100 Bioanalyzer to check the distribution of the size of fragments (quality), and qPCR to quantify the library. The qualified libraries were pair-end sequenced on a HiSeq X Ten platform (BGI-Shenzhen, China).

Quality control, annotation, and expression levels

To ensure the reliability of the subsequent bioinformatics analysis, the raw reads were filtered using SOAPnuke (Cock et al., 2010) with the following conditions: (1) reads that contained the adapter sequence was removed; (2) reads with unknown base N content >5% were removed; and (3) low-quality reads (quality value <10 for >20% of the total number of bases in the reads) were removed. The high-quality clean reads were used for all subsequent analyses. The mouse (Mus musculus) genome assembly (GRCm38.p6) and gene annotation files (GCF_000001635.26_GRCm38.p6) were downloaded from the National Center for Biotechnology Information (NCBI). An index of the reference genome was built using Bowtie2 (Langmead and Salzberg, 2012). The reads from each library were assembled independently using StringTie (Martin and Wang, 2011; Pertea et al., 2015). RSEM (Li and Dewey, 2011) (https://github.com/deweylab/RSEM) was used to calculate the expression levels of genes and transcripts by establishing a model of reads using the maximum likelihood estimation to allocate reads to different transcripts, then using a chain-specific model to distinguish the reads. Gene expression was standardized using FPKM (fragments per kilobase of transcript per million mapped reads), then used to detect differentially expressed genes among the libraries by DEseq2 (Love et al., 2014). Prediction of lncRNA cellular localization used LncLocator (Shanghai Jiao Tong University, http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/) (Cao et al., 2018).

Gene set enrichment analysis

GSEA was performed on BGI’s Dr.Tom system (http://biosys.bgi.com) with the software package provided by the Broad Institute, based on the MSigDB (Molecular Signatures Database). Normalized enrichment score (NES) and false discovery rate (FDR) were used to quantify the enrichment magnitude and statistical significance, respectively (|NES| is greater than 1.0, nominal p-value below 5%, and the false detection rate (FDR) below 25%).

Protein-protein interaction (PPI) network analysis

To illustrate the protein-protein interaction (PPI) information of the screened mRNAs, we mapped the mRNAs to the online tool-Search tool for the retrieval of interacting genes (STRING11) database, and NCBI Reference (https://www.ncbi.nlm.nih.gov/).

Only interactions that enjoyed a required combined score >950 were set as significant. The score range is 0–1000. Subsequently, the candidate key genes were identified based on the PPI network.

Construction of mRNA-lncRNA coexpression network

WGCNA was performed with the “WGCNA” package (Langfelder and Horvath, 2008). A total of 4 non-gray modules were generated. WGCNA was performed on BGI’s Dr.Tom system (http://biosys.bgi.com). The mRNA-lncRNA coexpression was used to clarify the interactions between the mRNAs and lncRNAs. The Pearson’s correlation coefficient (r) was calculated for each mRNA-lncRNA pair and significant pairs (r ≥ 0.9 or r ≥ 0.7) were selected for coexpression construction.

KEGG pathway analysis of DEGs

DElncRNAs and DEmRNAs were detected using the DEGseq software (Wang et al., 2010) (Love et al., 2014). Based on the following criteria: |FC| ≥ 2 and FDR ≤0.001. To clarify the potential functions of the DElncRNAs and their target genes, the gene sequences were searched against the non-redundant protein (NR), Swiss-Prot, Gene Ontology (GO), and Kyoto Encyclopedia of Gene and Genome (KEGG) databases (Kanehisa et al., 2017), using BLAST (Altschul et al., 1990) (E-value < 1.0E−5). Corrected P-values (FDR ≤ 0.05 and Q value ≤ 0.05) were used as the criterion to identify significant enrichment pathways. A heat map was created and enriched using R software (https://www.R-project.org/). The heatmap package in R was used to perform a hierarchical clustering analysis of the DEGs.

QPCR

Total RNA from BAT was extracted with TRIzol reagent (Invitrogen), according to the manufacturer’s instructions. The following PCR conditions were used: 1 cycle at 95°C for 30 s, then 40 cycles at 95 °C for 5 s, and 60 °C for 34 s. The qPCR was performed using a Bio-Rad CFX96 Touch Real-Time PCR Detection System (Bio-Rad Laboratories, Hercules, CA, USA). The primer sequences are presented in Supporting Information Table S16.

Quantification and statistical analysis

Normalization (reads per kilobase pair per million mapped reads, RPKM, and counts per million mapped reads, CPM) analysis was used to control the quality of the sequence data of mRNAs and lncRNAs in our study. The criteria for differentially expressed genes were |log2FC| ≥ 1, FDR ≤0.001. Data are expressed as the mean ± SD. Multiple comparisons among all groups were performed by one-way analysis of variance (ANOVA) using GraphPad Prism 9.0. Comparisons between two groups were performed by Student’s t-test using GraphPad Prism 9.0. p < 0.05 was considered statistically significant.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Brown adipose tissue (10–50 mg) was used for total RNA extraction according to the TRIzol protocol (TRIzol, Life Technologies)	Zhejiang University School of Medicine	N/A
Chemicals, peptides, and recombinant proteins
Paraformaldehyde (4%, pH 7.4)	New Cell & Molecular Biotech Co., Ltd	N1012
High-glucose Dulbecco’s modified Eagle medium	DMEM, Gibco	C11995500BT
Fetal bovine serum	Gibco	10100147
Critical commercial assays
TRIzol Reagent	Invitrogen, Carlsbad, CA, USA	15596026
Deposited data
Total RNA sequencing	This manuscript	GEO: GSE207064
Experimental models: Cell lines
Immortal preadipocyte lines	Prof. Zhuoxian Meng	Wang et al. (2014)
Experimental models: Organisms/strains
Mouse: C57BL/6	Shanghai SLAC Laboratory Animal Co., Ltd	RRID:MGI:2159769
Oligonucleotides
qPCR primers	This paper	See Table S16
Software and algorithms
R	https://www.cran.r-project.org	V4.1.0
DEGseq		RRID:SCR_008480; http://www.bioconductor.org/packages/2.6/bioc/html/DEGseq.html
GraphPad Prism 9.0	GraphPad Software	https://www.graphpad.com/
SOAPnuke	Cock et al. (2010)	RRID:SCR_015025; https://github.com/BGI-flexlab/SOAPnuke
Bowtie2	Langmead and Salzberg (2012)	RRID:SCR_016368; http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
StringTie	Martin and Wang (2011); Pertea et al. (2015)	RRID:SCR_016323; https://ccb.jhu.edu/software/stringtie/
RSEM	https://github.com/deweylab/RSEM	RRID:SCR_013027; http://deweylab.biostat.wisc.edu/rsem/
DEseq2	Love et al. (2014)	RRID:SCR_015687; https://bioconductor.org/packages/release/bioc/html/DESeq2.html
BGI’s Dr.Tom system	Broad Institute	http://biosys.bgi.com
NCBI	https://www.ncbi.nlm.nih.gov/	RRID:SCR_006472;
Weighted Gene Co-expression Network Analysis	Langfelder and Horvath (2008)	RRID:SCR_003302; http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/
NCBI BLAST	Altschul et al. (1990)	RRID:SCR_004870; http://blast.ncbi.nlm.nih.gov/Blast.cgi
Bio-Rad CFX96 Touch Real-Time PCR Detection System	Bio-Rad Laboratories, Hercules, CA, USA	RRID:SCR_018064; https://www.bio-rad.com/en-us/product/cfx96-touch-real-time-pcr-detection-system?ID=LJB1YU15
Other
High Fat Diet (HFD; 60% kcal fat)	Research diets	Cat# D12492
Chow Diet	Xietong Biology	SW9102