Identification of differentially expressed mRNA and the Hub mRNAs modulated by lncRNA Meg3 as a competing endogenous RNA in brown adipose tissue of mice on a high-fat diet

ABSTRACT Obesity is associated with insulin resistance, diabetes, and obesity-related metabolic disorders. Brown adipocytes have emerged as potential targets for the treatment of obesity and obesity-related diseases. However, changes that occur in brown adipose tissue during various stages of high fat diet (HFD)-induced obesity remain poorly understood. The present study aimed to determine the changes occurring in brown adipose tissue during various stages of an HFD by analyzing two microarray expression profiles. A total of 1,337 differentially expressed RNAs (DE RNAs) were identified between the HFD and ND groups, using the limma package in R. The DE RNAs included 1,249 mRNAs, 74 long non coding RNAs (lncRNAs), and 14 pseudogenes. Functional annotation of the DE mRNAs, including GO terms and KEGG pathways were identified using the Database for Annotation, Visualization, and Integrated Discovery. A protein-protein interaction network was constructed using STRING and clusters were obtained through the Molecular Complex Detection plug-in. In the present study, the lncRNA,maternally expressed gene 3 (Meg3), was identified as the DE lncRNA with a significant fold change. The network of Meg3 as a ceRNA was constructed, which demonstrated that Meg3 modulated five hub DE mRNAs via competitive binding to microRNAs.


Introduction
A high-fat-diet (HFD) can induce obesity, which is closely associated with insulin resistance, diabetes, and obesityrelated metabolic diseases [1][2][3]. An increasing number of studies have focused on the impacts of an HFD on white adipose tissue, not only because of its known function as an excess energy reserve, but also because of its function of secreting adipokines during metabolic processes [4]. However, another type of adipose tissue, known as brown adipose tissue (BAT), can dissipate excess energy via nonshivering thermogenesis. During this process, uncoupling protein 1 (UCP1), a unique BAT mitochondrial membrane protein, uncouples respiration from ATP synthesis to produce heat, and stimulates high levels of fatty acid oxidation to release energy [5]. BAT was previously thought to exist only during infancy; however, it has now emerged as a potential therapeutic target to treat obesity and obesityrelated diseases, since BAT has also been found to exist in the dorsal interscapular region of adult humans [6,7].
Conventionally, RNA can be classified into two types according to its protein-coding potential. The first type is coding RNA, which can be translated into proteins that determine distinct phenotypes and biological functions [8] and the second type is non-coding RNA, which lacks protein-coding ability. Non-coding RNA can be further divided into subtypes, including microRNA (miRNA) and long non-coding RNA (lncRNA). LncRNAs were once considered to be transcriptional noise, because they consist of transcripts exceeding 200 nucleotides in length, but with no protein-coding potential. However, in recent years, accumulating evidence has demonstrated that lncRNAs act as competing endogenous RNAs (ceRNAs), which posttranscriptionally modulate mRNAs by competitively binding to microRNAs using complementary base pairing [9,10].
Several lncRNAs, including maternally expressed gene 3 (Meg3), have been identified as ceRNAs implicated in pathophysiological processes of metabolic disorders [11][12][13]. Meg3 has been shown to serve as a ceRNA to regulate the miR-302a-3p-CRTC2 axis and the miR-214-ATF4 axis, and promote hepatic insulin resistance [11,12]. Chen et al. reported that Meg3, acts as a ceRNA by directly binding to miR-145 in cardiomyocytes, thereby upregulating the expression of PDCD4 to induce cardiomyocyte apoptosis under high-glucose conditions [13]. However, changes in Meg3 levels and its corresponding function in BAT under HFD conditions remain ill-defined.
The previous researches pay close attention to the alteration of thermogenic-related mRNAs and lncRNAs during the differentiation of brown preadipocytes or in the BAT exposing to cold temperature [14,15]. To better understand changes in mRNA levels, the functional enrichment of differentially expressed (DE) mRNAs, and the exact function of lncRNAs as ceRNAs in BAT under HFD conditions, two microarray expression profiles were retrieved from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm. nih.gov/geo/), an international public repository of highthroughput microarray data and relevant functional genomic data sets [16]. Data from BAT samples of male C57BL/ 6 mice fed an HFD or a normal diet (ND) were collected from two datasets. A total of 1,337 DE RNAs were identified, including 1,249 mRNAs, 74 lncRNAs, and 14 pseudogenes, between the HFD and ND groups, using the limma package in R. Functional annotation of the DE mRNAs, including gene ontology (GO) terms and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways, was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID). A protein-protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) and clusters were obtained using the Molecular Complex Detection (MCODE) plugin. The lncRNA, Meg3, was identified as the lncRNA with a significant fold change. A network was constructed with Meg3 as a ceRNA, which demonstrated that Meg3 modulated several mRNAs by competitively binding to microRNAs. Five hub DE mRNAs were identified in the Meg3 ceRNA network. These results may deepen our understanding of the BAT function in HFD-induced obesity and obesity-related disorders such as diabetes.

Ethical declaration
This research does not directly contain any material obtained from animals or humans. All data used in this study were extracted from public databases.

Microarray data archives
Microarray expression profiles from GSE113315 and GSE116225 datasets were retrieved from the GEO database.
Beginning at 3 weeks of age, the male C57BL/6mice were fed regular chow diet (ND, D12450B, Research Diet, Inc. New Brunswick, NJ, USA) or HFD (D12492, Research Diet, Inc. New Brunswick, NJ, USA) ad libitum for 15 weeks [17,18]. ND serves as a negative control. The mice were kept in enclosures at room temperature with a fixed 12 h light/12 h dark cycle. The expression profiles of GSE113315 were based on the GPL13112 (Illumina HiSeq 2000 [Mus musculus]) platform and the expression profiles of GSE116225 were based on the GPL17021 (Illumina HiSeq 2500 [Mus musculus]) platform. The series matrix files of the two datasets were downloaded from GEO to screen for DE RNAs between BAT samples from the HFD and ND groups.

Microarray data and DEG identification
Following annotation of the matrix files from Ensemble ID to gene symbols using http://grch37.ensembl.org/Mus_ musculus/Info/Index, the sva package in R (version 3.6.0; University of California, Berkeley, CA, USA) was applied to correct background expression values and normalize the data [19]. DE RNAs, including mRNAs and lncRNAs, with the threshold criteria of adjusted p < 0.05 and |log fold change (FC)| > 1 between the HFD and ND group, were screened via the limma package in R [20]. DE RNA biotypes were identified using the online resource at http://grch37. ensembl.org/Mus_musculus/Info/Index. The pheatmap package in R was subsequently used to plot mRNA and lncRNA heatmaps [21].

GO and pathway enrichment analyses
GO and KEGG pathway analyses were performed using DAVID 6.8 (http://david.ncifcrf.gov) [22]. GO is a commonly used bioinformatics tool that provides comprehensive information on the function of individual gene products, based on defined features. GO analysis involves the identification of biological processes (BP), cellular components (CC), and molecular functions (MF). KEGG is a major database used to understand high-level biological functions and utilities. A gene count >2 and p < 0.05 were set as the threshold values.

PPI network creation and hub gene identification
A PPI network of DE mRNAs was constructed using STRING 11.0 (https://string-db.org/), with a combined score > 0.9 as the cut-off value [23]. Significant modules in the PPI network were identified using MCODE 1.5.1, a Cytoscape software plug-in [24]. The parameter for DE mRNA clustering and scoring was set as follows: MCODE score ≥ 4, degree cut-off = 2, node score cut-off = 0.2, max depth = 100, and k-score = 2.

Target miRNA prediction of Meg3 and corresponding target mRNA prediction
Meg3 was selected to be further investigated for two reasons. Firstly, Meg3 was the DE lncRNA with a significant fold change (Table S1). Secondly, the full sequence of Meg3 was available in the NONCODE database (http://www.noncode.org). The Meg3 sequence was imputed into miRDB (http://www. mirdb.org) to obtain the target miRNAs, according to the principle of complementary base pairing. The top five target scores were chosen as the target prediction results. The corresponding target mRNAs of the predicted miRNAs were identified using the Targetscan database (http://www.targetscan.org). From the list of mRNAs in both the Targetscan results and the DE mRNAs, the top 20 mRNAs were selected according to a comprehensive rank of 'cumulative weighted con-text++ scores'.

Construction of Meg3 ceRNA regulatory network
Cytoscape software (Cytoscape, 3.7.1) was used to construct the network of Meg3 as a ceRNA. The network was visualized using the Meg3-miRNA-mRNA triple competing relationship module.

Statistical analysis
Statistical analyses of DE RNAs were performed using R, with the threshold criteria of adjusted p < 0.05 and | log FC| > 1 between the HFD and ND groups. A gene count > 2 and p < 0.05 were set as the thresholds in the GO and KEGG analyses. The parameters for DE mRNA clustering and scoring in the significant module analyses were set as follows: MCODE score ≥ 4, degree cut-off = 2, node score cut-off = 0.2, max depth = 100, and k-score = 2.

Identification of HFD-induced DE RNAs in BAT
To identify DE mRNAs and DE lncRNAs in BAT between the HFD and ND groups, we retrieved relevant microarray expression profiles from the GSE113315 and GSE116225 datasets in the GEO database. After the consolidation and normalization of microarray data, 1,337 DE RNAs, including 1,203 upregulated RNAs and 134 downregulated RNAs, were screened using the limma package (|logFC| > 1, adjusted p < 0.05) ( Figure 1, Table S2). The DE RNA biotypes were identified and found to include 1,249 DE mRNAs, 74 DE lncRNAs, 14 DE pseudogenes between the HFD and ND groups, as shown in the heatmap (Figures 2 and 3).

GO enrichment analysis of DE mRNAs
To determine the biological features of the DE mRNAs, GO analysis was performed using the DAVID online tools. The BP analysis showed that the DE mRNAs were mostly enriched for cell adhesion, angiogenesis, and inflammatory response terms ( Figure 4). The CC analysis showed that the DE mRNAs were significantly enriched for membrane, cell surface, and extracellular exosome terms ( Figure 4) and the MF analysis showed that the DE mRNAs were significantly enriched for protein binding, integrin binding, actin binding, and ATP binding terms ( Figure 4).

KEGG enrichment analysis of DE mRNAs
To explore the potential mechanism responsible for these DE mRNAs, KEGG pathway analysis was performed using DAVID online tools. The results of the KEGG analysis showed that DE mRNAs were mainly involved in focal adhesion, leukocyte transendothelial migration, Staphylococcus aureus infection, and ECMreceptor interactions ( Figure 5).

PPI network analysis
To identify the most significant clusters of the DE mRNAs, a PPI network of DE mRNAs was constructed using STRING. As shown in Figure 6(a), there were 1,238 nodes and 2,429 edges in the PPI network. There were 79 DE mRNAs in the most significant module (score = 17.838) recognized by MCODE ( Figure 6(b), Table 1).

Construction of the Meg3 network as a ceRNA in BAT under HFD conditions
To better understand the function of DE lncRNAs on coding RNAs, Meg3 was selected to construct a ceRNA network ( Figure 7, Table 2). The overlapping genes between the most significant module of MCODE and mRNAs modulated by Meg3, including Stk10, Leprel1, Itgam, Fam20 c, and Col6a2, were identified as hub genes modulated by Meg3 as a ceRNA in BAT under HFD conditions. Meg3 could compete with Stk10, Leprel1, and Itgam mRNA for binding to mmu-miR -466i-5p miRNA; compete with Fam20 c mRNA for binding to mmu-miR-574-5p; and compete with Col6a2 mRNA for binding to mmu-miR-770-5p, thus upregulating the expression of these five mRNAs in BAT under HFD conditions (Figure 8).

Discussion
Although the function of BAT has been demonstrated in previous studies, the effect of an HFD on the GO annotations and KEGG pathways of mRNA and lncRNA-miRNA-mRNA crosstalk in BAT remain illdefined. In the present study, based on two microarray expression profiles analysed using bioinformatic methods, DE mRNAs and lncRNAs were identified between BAT under HFD and ND conditions. GO analysis indicated that the DE mRNAs were significantly enriched in the BP terms, cell adhesion, angiogenesis, and inflammatory response. Eguchi and colleagues reported that adipocyte adhesion molecule was implicated in adipocyte maturation and the development of obesity and thus, they speculated that cell adhesion may influence the morphology and differentiation of cells via alterations in cell signalling or cytoskeletal organization [25]. However, to the best of our knowledge, no previous studies have demonstrated cell adhesion functions in BAT. Broad-spectrum findings have reported that hypoxia plays a critical role in adipose tissue angiogenesis in response to obesity, which is rapidly induced by HFD [26][27][28]. This is consistent with another BP term enriched in response to hypoxia. It is noteworthy that, within the enrichment in response to hypoxia, uncoupling protein 2 (UCP2) expression was up-regulated, but UCP1 was not. UCP1 is linked to the protection against diet-induced obesity (DIO) as an integral membrane protein unique to BAT mitochondria [29]. However, UCP2 seems to have a different function, since Kim et al. demonstrated that the regulation of microglial UCP2 in vivo is associated with increased levels of inflammatory cytokines and that the deletion of microglial UCP2 prevents DIO [30]. According to their results and the results of the present study, we propose that UCP2 may be an inflammatory inducer in response to hypoxia in BAT under HFD conditions. Further studies are required to clarify the mechanism for the regulation of UCP2 in BAT under HFD conditions, since UCP is a promising therapeutic target for the treatment of obesity and obesity-related diseases.
In GO analysis, the significant DE mRNAs were enriched in the CC terms, membrane, cell surface, and extracellular exosome. Although exosomes derived from normal BAT can improve glucose tolerance [31,32], the functions of exosomes derived from BAT under HFD conditions are ill defined, but these exosomes are an emerging therapeutic target for metabolic diseases.
The GO analysis of MF suggested that the DE mRNAs were most significantly enriched in protein binding, which suggested that the ability of proteins to form bonds with other substances was increased.   Interestingly, KEGG enrichment analysis of DE mRNAs showed that these DE mRNAs were enriched in focal adhesion, which was consistent with the GO analysis of BP, which suggested that cell adhesion may play a critical role in the BAT under HFD conditions. The present study is the first study to demonstrate the upregulation of Meg3 in BAT under HFD conditions and to construct the corresponding network of Meg3 as a ceRNA (Figure 7). Increasing evidence has implicated Meg3 in metabolic disorders. Several studies have demonstrated that Meg3 promotes insulin resistance by serving as a ceRNA of miRNAs, to upregulate mRNA expression and consistently, Meg3 knockdown alleviates insulin resistance in palmitate-treated hepatocytes and in mice fed an HFD [12,33]. However, You et al. demonstrated that the suppression of Meg3 levels decreases insulin secretion from pancreatic islet beta cells [34]. Therefore, the function of Meg3 in the development of metabolic disorders may be tissue specific. It would be of great significance to identify Meg3 as a significant ceRNA in BAT under HFD conditions and to determine the downstream mechanism whereby Meg3 in BAT affects the development of obesity and obesity-related diseases in response to an HFD.
In the present study, three miRNAs, including miR-466i-5p, miR-574-5p and miR-770-5p were validated as the target miRNAs of Meg3, predicted by the bioinformatics analysis. MiRNAs are short RNAs that can be regulated by the sponge role of lncRNA, following by modulating the expression of the downstream genes by targeting their 3′untranslated regions (3′UTR) [35]. An increasing number of studies has demonstrated that miR-574-5p can be sponged by lncRNA, like lncRNA PTCSC3 and lncRNA MFI2-AS1, to modulate growth and metastasis of cancer cells [36,37]. In addition, miR-770-5p has been found to be downregulated by lncRNA TPT1-AS1 and correspondingly upregulated the expression of STMN1, leading to promotion of the proliferation of Glioma cells [38].
Interestingly, the genes targeted by miR-466i-5p, miR-574-5p and miR-770-5p, were identified as hub genes modulated by Meg3 as a ceRNA in BAT under HFD conditions, including Stk10, Leprel1, Itgam, Fam20 c, and Col6a2. Among them, Leprel1, Itgam and Col6a2 are expressed in adipose tissues and involved in the regulation of inflammation, fibrosis, insulin signalling and so on, respectively [39][40][41]. There are no findings linking Stk10 and Fam20 c with adipocytes or adipose tissues, although Stk10      possesses anti-apoptotic property [42,43] and Fam20 c is associated with insulin production in pancreases β cells [44]. Furthermore, to date, there are no previous reports correlating these five genes with Meg3.
Further studies are thus required to validate the correlation between Meg3 and the five hub genes identified in the present study and to determine the mechanisms whereby Meg3 and these five hub genes are involved in metabolic processes and the development of metabolic disorders. Taken together, the present study, for the first time, identifies changes in both mRNA and lncRNA levels contributing to HFD-induced obesity and indicates that lncRNA Meg3 might modulate several mRNAs via binding to microRNAs competitively. However, this study is a purely bioinformatics analysis aiming to help investigators to explore the DE RNAs and DE LncRNAs in brown adipose tissue in mice fed with HFD. The conclusions are required to be verified by in vitro and in vivo experimental settings closer to the clinical reality or in the invidious with obesity or diabetes. Only in this case can we better understand the exact functions of Meg3 and its five downstream genes in human BAT and also uncover mechanisms related to the development and progression of obesity and diabetes.