Analysis of differentially expressed long non-coding RNAs in LPS-induced human HMC3 microglial cells

Background Long non-coding RNAs (lncRNAs) are emerging as key modulators of inflammatory gene expression, but their roles in neuroinflammation are poorly understood. Here, we identified the inflammation-related lncRNAs and correlated mRNAs of the lipopolysaccharide (LPS)-treated human microglial cell line HMC3. We explored their potential roles and interactions using bioinformatics tools such as gene ontology (GO), kyoto encyclopedia of genes and genomes (KEGG), and weighted gene co-expression network analysis (WGCNA). Results We identified 5 differentially expressed (DE) lncRNAs, 4 of which (AC083837.1, IRF1-AS1, LINC02605, and MIR3142HG) are novel for microglia. The DElncRNAs with their correlated DEmRNAs (99 total) fell into two network modules that both were enriched with inflammation-related RNAs. However, treatment with the anti-inflammatory agent JQ1, an inhibitor of the bromodomain and extra-terminal (BET) protein BRD4, neutralized the LPS effect in only one module, showing little or even enhancing effect on the other. Conclusions These results provide insight into, and a resource for studying, the regulation of microglia-mediated neuroinflammation and its potential therapy by small-molecule BET inhibitors. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-09083-6.

such as JQ1 and I-BET151 protected mice against LPSinduced sepsis [7] and attenuated the LPS-stimulated expression of pro-inflammatory genes in murine bone marrow-derived macrophages [8]. Similarly, JQ1 attenuated pro-inflammatory chemokine, cytokine, and interferon response genes in the LPS-treated murine glial cell line BV-2 [9], an observation that we confirmed and expanded in the human microglial cell line HMC3 [10].
Thus, it appears that the BET proteins also regulate neuroinflammation.
The present study therefore extends our transcriptomic analysis of the LPS/JQ1-treated HMC3 cells into the category of differentially expressed lncRNAs (DElncRNAs). By re-analyzing the published RNA-seq datasets [10], we identified (i) inflammation-related and BET inhibitorsensitive DElncRNAs, (ii) their correlated DEmRNAs, and (iii) potential functional networks that contain these two classes of transcripts.

Approach
We first identified the LPS-induced DElncRNAs in the dataset GSE155408. Compared to the initial study [10], we added processing methods that are specifically designed to identify DElncRNAs (Fig. 1A); furthermore, we now used a stricter cutoff, which slightly downsized the DEmRNA data set ( Supplementary Fig. 1A). We then identified the DEmRNAs that were correlated with the LPS-induced DElncRNAs. In parallel, we assessed whether or how the LPS-induced DE RNAs (including lncRNAs and mRNAs) were affected by JQ1. Next, we performed a network analysis in order to construct interaction modules incorporating the LPS-induced DElncR-NAs and their correlated DEmRNAs. This was followed by pathway analysis in order to learn about the potential functional significance of these modules. Finally, we constructed a corresponding protein-protein interaction network.

Identification of DElncRNAs and DEmRNAs in LPS-treated HMC3 cells
In the control vs. LPS-treated cells, we identified a total of 5 DElncRNAs (all upregulated; Fig. 1B and C) and 99 DEmRNAs (98 upregulated and 1 downregulated; Supplementary Fig. 1A). Conversely, the DElncRNAs and a majority of the DEmRNAs tended to be downregulated or unaffected by JQ1. When LPS and JQ1 were combined, JQ1 further increased the expression of the two DElncR-NAs that were most increased by LPS (AC083837.1 and LINC02605), but partially (MIR3142HG) or fully (MIR155HG and IRF1-AS1) neutralized the LPS effect on the other DElncRNAs (Fig. 1C). Likewise, JQ1 counteracted LPS or did not alter its effect for the majority of DEmRNAs, but enhanced the LPS effect in some cases. We note that the above DElncRNAs had not yet been identified in the HMC3 cells; the DEmRNAs are often inflammation-and immunity-related (e.g., CCL20, CSF3, CXCL10, TNF, and CXCL8) (Supplementary Fig. 1B and C). The up-and downregulated DElncRNAs and DEmR-NAs are listed in Supplementary Table 1.

Correlations between DElncRNAs and DEmRNAs
The heat map in Fig. 2

Identification of DElncRNA-DEmRNA network modules
In order to find potential interactions, we constructed co-expression networks based on the above correlations between the DElncRNAs and DEmRNAs. The resulting nodes and relations fell into two modules, which we dubbed "large turquoise" and "small cyan" (Fig. 3). The large turquoise module was defined by IRF1-AS1 and MIR155HG (the two DElncRNAs that were least upregulated by LPS); the small cyan module was defined by AC083837.1, LINC02605, and MIR3142HG (the DElncR-NAs that were more strongly induced by LPS). Within these modules, MIR155HG has the maximum number of co-expressed genes.
We determined by chi-square test whether there was a significant difference in how the two modules respond to the JQ1 treatment (Supplementary Table 3). The mRNAs of the large turquoise module were significantly related to JQ1, while the mRNAs of the small cyan module were not related to JQ1. Specifically, out of 65 mRNAs in the large turquoise module, the levels of 64 mRNAs were decreased by JQ1. In contrast, out of 29 mRNAs of the small cyan module, 11 were marginally affected and 18 were increased. In more detail, the mRNA levels of the small cyan module that were highly correlated with MIR3142HG were reduced or marginally affected by JQ1, while the mRNAs that were highly correlated with AC083837.1 and LINC02605 showed increased expression or marginal changes. To validate the two modules, which were identified on the basis of RNA-seq data, we repeated the LPS/JQ1 treatments of the HMC3 cells and measured the expression levels of the DElncRNAs and selected (i.e., highly correlated, inflammation-related) DEmRNAs. In the large turquoise module, JQ1 mostly neutralized the LPS effect (22 of 23 tested mRNAs, both tested lncRNAs) ( Fig. 4 and Supplementary Fig. 2), whereas in the small cyan module, most of the LPS-stimulated RNAs were marginally affected or even further increased by the additional presence of JQ1 (15 of 17 mRNAs, 1 of 2 tested lncRNAs) To summarize, while the mRNAs of the large turquoise module mostly were decreased by JQ1, the mRNAs of the small cyan module showed a mixed response pattern.

Functional annotations
Having validated the network modules, we performed a functional classification and pathway enrichment analysis of their respective mRNAs. In both the large turquoise ( Fig. 6A) and small cyan ( Fig. 6C) modules, Gene Ontology (GO) analysis highlighted terms related to inflammatory response such as type I interferon signaling pathway, defense response to virus, and chemokine-mediated signaling pathway. Similarly, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis ( Fig. 6B and D) found a total of 32 KEGG pathways

Protein-protein interaction (PPI) network analysis
The PPI network mathematically calculates the physical interactions between proteins in cells, which allowing a molecular assessment at the molecular and system level. We constructed a PPI network of DEmRNAs, which highly correlated DElncRNAs using the STRING database (Fig. 7). The PPI network contained 63 nodes and 332 edges in the co-expressed mRNAs of the large turquoise module (related to the MIR155HG and IRF1-AS1 lncRNAs) and 29 nodes and 90 edges in the co-expressed mRNA of the small cyan module (related to the AC083837.1, LINC02605, and MIR3142HG lncRNAs). The top 5 mRNAs that had the highest node degrees in the large turquoise module were IL1B (degree = 33), IRF1 (degree = 31), IL6 (degree = 29), CXCL10 (degree = 28), and TNFSF10 (degree = 27). The top 5 mRNAs that had the highest node degrees in the small cyan module were TNF (degree = 21), CXCL8 (degree = 18), NFK-BIA (degree = 18), CXCL1 (degree = 12), and CXCL2 (degree = 12).

Discussion
In this study, we identified 5 DElncRNAs in LPS-treated HMC3 human microglial cells. 4 of these transcripts (AC083837.1, IRF1-AS1, LINC02605, and MIR3142HG) have not previously been invoked in the microglial inflammatory response. We further found that the BET inhibitor JQ1 can exert both positive and negative effects on the LPS-induced DElncRNAs. Importantly, we were able to correlate the novel DElncRNAs with inflammation-related DEmRNAs, in line with the idea that the DElncRNAs modulate the microglial immune response. More specifically, we found that the LPSinduced DElncRNAs and their correlated mRNAs fall into two modules that we dubbed "large turquoise" (defined by IRF1-AS1 and MIR155HG) and "small cyan" (defined by AC083837.1, LINC02605, and MIR3142HG). Interestingly, JQ1 affected their patterns in opposite ways: Inflammation processes were decreased by JQ1 in the large turquoise module, but marginally affected or even significantly increased in the small cyan module. We performed qRT-PCR to verify the expression of lncRNAs and mRNAs of the divided modules based on the Pearson correlation coefficient (r) of DElncRNAs-DEmRNAs. We found some genes with conflicting DElncRNAs-DEmR-NAs correlations as a result of qRT-PCR, which appeared as a probability of false-positive in Pearson's correlation analysis. Pearson correlation analysis is the commonly used analysis method to measure the strength of association in a linear relationship between two datasets based on the mapped read count value. Although the Pearson correlation analysis result is reliable, it can have false positives due to the influence of several variables such as the magnitude of the slope at which the points are clustered and the heteroscedasticity [23], and thus results of RNAseq analysis and qRT-PCR may be different. Therefore, verification of qRT-PCR of RNA-seq data is necessary.
As a result of co-expression network analysis, PPI network analysis, and functional annotation analysis, we found that genes and biological terms related to immune response, inflammatory responses, cytokine, and chemokine activity were enriched in both modules, and TNFα and NF-κB signaling-related genes were also found in both modules. Interestingly, in terms of inflammatory response, the large turquoise module was enriched with interferon signaling-related functions, and the small cyan module was enriched with the CXC chemokine in particular. Regarding cell proliferation, survival, and metabolism, the large turquoise module was associated with negative regulation of metabolic processes, and positive regulation of apoptotic processes such as TRAIL signaling and necroptosis. On the other hand, the small cyan module showed positive regulation of metabolic processes and more enriched MAPK signaling regulation including activated tak1 mediates p38 mapk activation compared to the large turquoise module.
As summarized above, TNFα-related genes (including also some TNF family members themselves) and NF-κBrelated genes (including also NF-κB components) are present in both modules. As such, this is in line with the known roles of the TNFs and NF-κB as global regulators of inflammation, i.e. as regulators of multiple aspects of inflammation. However, the response of the TNFα-or NF-κB-related genes to JQ1 is different between the two modules, such that the TNFα-or NF-κB-related genes are downregulated along with the IFN-related genes in the large turquoise module but upregulated or unchanged along with the CXC chemokine genes in the small cyan module. Hence, different aspects of inflammation regulated by TNFs or NF-κB respond differently to JQ1, hence are under distinct epigenetic controls. Given that the TNFs and NF-κB are each related to more than one signaling pathway [24], it will be of interest to relate the large turquoise vs. small cyan modules to those specific pathways.
A complementary message emerges from the findings that we made with the more specialized inflammatory regulators, namely the IFNs and CXC chemokines. The direct biological effects of IFNs (mediated mainly by JAK/STAT signaling) are related especially, albeit not exclusively, to the modulation of immune system [25]; the direct biological effect of CXC chemokines is especially, albeit not exclusively, to promote neutrophil migration [26].
Here, we found that along with the TNFα/NF-κB genes of the large turquoise module, the IFN (mainly type I and II)-related genes (including IFNB1), which were limited to that same module, were also inhibited by JQ1. This agrees with the literature, since the direct inhibition of IFN response by BET inhibitors has been well documented [27]. Furthermore, we found that along with the TNFα/NF-κB-related genes of the small cyan module, the CXC chemokine genes (encoding CXCL1/2/3/5/8, all known to bind the chemokine receptor CXCR2), which were limited to that same module, were not inhibited by JQ1. We note in this respect that the GO term "activated tak1 mediated p38 mapk signaling activation" was also enriched in the small cyan module, in line with the known role of tak1/p38 signaling in CXC chemokine signaling [28]. Thus, not the type of cytokine per se (TNF vs. IFN vs. CXC chemokine) determines the outcome of JQ1 treatment, but the belonging to the large turquoise vs. small cyan module.
The selectivity of the gene response to I-BET has been reported to be related to the epigenetic status of the responding gene and the mechanism of BET recruitment. Nicodeme E et al. reported that some LPS-induced cytokines and chemokines, such as Tnf, Ccl2-5 and Cxcl1/2, have a highly selective effect on I-BET in bone marrow-derived macrophages. Further study on the mechanism of selectivity for BET inhibitors is needed, and considering the previous reports, selectivity for BET inhibitors may be related to the BET recruitment pathway or histone acetylation level for each gene, and the binding of the BET protein to the gene promoter or super-enhancer.
Of note, the DElncRNAs that we identified in the LPSstimulated HMC3 cells have been observed in other inflammation-related contexts before. IRF1-AS1 (Lnc-SLC22A5-6) was found to act as a positive modulator of the IFN response in esophageal squamous cell carcinoma (ESCC), functioning as a tumor suppressor by regulating cell proliferation and apoptosis [29]. MIR155HG was increased in human M1 (inflammatory-type) macrophages that were derived from monocytes by treatment with LPS and IFN-gamma [30]. In human B cell lymphoma cells, LPS induced the nuclear translocation of an NF-κB p50/p65 heterodimer that could bind to the MIR155HG promoter, suggesting that MIR155HG was a direct NF-κB target gene [31]. MIR3142HG, the host gene for miR-3142 and miR-146a, was reported to regulate the IL-1β-induced inflammatory response in idiopathic pulmonary fibrosis (IPF) [32] and to be highly expressed in LPS-exposed human pulmonary microvascular endothelial cells (HPMECs) [33]. It is therefore of interest that our study showed that the LPS-induced increase of MIR3142HG was reduced by JQ1. However, the expression levels of MIR3142HG-correlated genes are controversial. In our study, JQ1 was not or only marginally effective on the MIR3142HG-correlated genes. Finally, LINC02605 (IL7-AS; Lnc-ZC2HC1A-1) was reported to regulate the immune response in a cell-type specific manner [34]. IL7-AS is a positive regulator of the IL1β-induced inflammatory response in human A549 epithelial cells, but a negative regulator in LPS-stimulated human THP-1 monocytes and mouse RAW 264.7 macrophages. Additionally, knockdown of IL7-AS increased the IL-6 release in IPF-derived fibroblasts, indicating that it is a negative regulator [32]. Our co-expression network analysis showed that while LINC02605 was induced by LPS, neither LINC02605 nor its correlated genes were significantly affected by JQ1. Interestingly, this was observed for all three DElncRNAs and their correlated genes in the small cyan module. These results are consistent with our previous transcriptome analyses of human and mouse microglial cells and mouse bone marrowderived macrophages [8][9][10].

Conclusion
Here, we present a comprehensive analysis of inflammation-related DElncRNA and DEmRNA expression profiles and functional networks in the human microglial cell line HMC3. We identified 5 DElncRNAs (including 4 novel ones in microglia) and 99 DEmRNAs. We constructed DElncRNA-DEmRNA co-expression networks, which fell into two separate modules, and investigated their functions and pathways, which -for both modules-turned out as largely known to be inflammationrelated. We determined that although considered as an anti-inflammatory agent, the BET inhibitor JQ1 regulates the two modules differently, showing an almost uniform anti-inflammatory effect on one module, but little or even enhancing effect on the other. This interesting result will have to be complemented and validated with methods such as knockdown and overexpression; it is also of interest whether it can be extended to other models of neuroinflammation. Altogether, the RNA expression modules that we identified here provide a resource for further studies of human microglial neuroinflammation through both computational analysis and functional approaches.

Identification of differentially expressed lncRNAs and mRNAs
For this study, we used the RNA-seq data of our previous paper (GSE155408). To identify DElncRNAs, a comprehensive reference list of known lncRNAs was included in the processing of the RNA-seq data [10,35,36]. Briefly, FASTQ data were quality controlled and trimmed with Trimmomatic (version 0.36) [37]. The FASTQ files were aligned using STAR (version 2.7.8) [38] alignment software with the GENCODE Homo sapiens reference sequence GRCh38 (Release 27). DElncRNAs and DEmRNAs were normalized to sequencing depth and RNA composition using the median method with default parameters of DESeq2 [39]. Differential expression analysis of lncRNA and mRNA was conducted using the DESeq2 R package. The DElncRNAs and DEmR-NAs were selected with a cutoff of | log 2 fold change (log 2 FC) | ≥ 1.2, | log 2 FC | ≤ -1.2, and adjusted p-value (padj) ≤ 0.05 in LPS-treated HMC3 cells.

Weighted gene co-expression network analysis (WGCNA)
First, the Pearson correlation coefficient (r) values were calculated to assess the similarity of the expression patterns of transcripts. Then, a scale-free network was obtained by weighting the correlation coefficient between transcripts with soft-thresholding power. A module is defined as a cluster of densely interconnected transcripts in terms of co-expression. We considered a |r|≥ 0.75 as a meaningful value. Cytoscape MCODE plug-in (Version 3.4.0, available online: http:// www. cytos cape. org/) [40] was applied for visualization of the co-expression networks.

Functional annotation and canonical pathway analysis
Database for Annotation, Visualization, and Integrated Discovery (DAVID, version 6.8) software (http:// david. abcc. ncifc rf. gov/ home. jsp) was used to analyze the biological functions in the datasets [41]. DAVID uses a modified Fisher's exact p-value to examine gene ontology (GO) enrichment. A false discovery rate (FDR) ≤ 0.05 was used as the criterion for GO term analysis. The KEGG Orthology Based Annotation System (KOBAS, version 3.0) software (http:// kobas. cbi. pku. edu. cn/) [42] was used to analyze the enriched KEGG pathways [43] in the datasets. FDR ≤ 0.05 was used as the criterion for KEGG pathway enrichment analysis.

Cell culture and treatment
HMC3 human microglial cells were purchased from the Korean Cell Line Bank (Seoul, Korea). The cells were cultured in minimum essential medium (MEM) supplemented with 10% fetal bovine serum (FBS), 100 IU/ ml penicillin, and 10 μg/ml streptomycin and were maintained in a humidified incubator at 37 °C with 95% air/5% CO 2 . The cells were treated with 100 ng/ml LPS (Sigma-Aldrich, St. Louis, MO, USA) and/or 500 nM JQ1 for 4 h under standard culture conditions. The LPS and JQ1 were dissolved in dimethyl sulfoxide (DMSO; Sigma-Aldrich, St. Louis, MO, USA).

Quantitative RT-PCR
Total RNA extractions and cDNA preparation were performed according to the manufacture's instruction (Takara, Shiga, Japan). Quantitative Reverse Transcription PCR (qRT-PCR) was performed using an ABI 7500 real-time PCR system (Applied Biosystems Inc., Foster City, CA, USA). The critical threshold (△CT) value was normalized by the expression of an internal control, glyceraldehyde-3-phosphate dehydrogenase (GAPDH). Finally, the results were also analyzed using the comparative critical threshold (△△CT) method. The primers were designed using Primer Bank (http:// pga. mgh. harva rd. edu/ prime rbank/ index. html) and are listed in Supplementary Table 4.

Protein-protein interaction (PPI) network analysis
The Search Tool for the Retrieval of Interacting Genes (STRING, http:// string. embl. de/) [44] was used to construct the PPI network for DEmRNAs (minimum required interaction score > 0.7). The interaction relationships of the proteins encoded by DEmRNAs were searched by STRING online software, and the combined score > 0.7 was used as the cut-off criterion. Cytoscape MCODE plug-in (Version 3.4.0, available online: http:// www. cytos cape. org/) was applied for visualization of the protein-protein interaction. Additionally, the network analyzer was used to compute the basic properties of the PPI network, including average clustering co-efficient distribution, closeness centrality, average neighborhood connectivity, node degree distribution, shortest path length distribution, and topological coefficients.

Statistical analysis
All data are expressed as the mean ± standard deviation of the mean (SD). The chi-square test was performed to confirm whether there was a statistically significant relationship between categorical variables by JQ1 treatment. A p-value or padj ≤ 0.05 was considered significant. The statistical analyses were performed using IBM SPSS Statistics ver. 26.0 (IBM Corporation, Armonk, NY, USA). All qRT-PCR data were tested using one-way ANOVA followed by Tukey's honestly significant difference (HSD) post hoc test. Differences for which p ≤ 0.05 were considered significant.

Availability of data and materials
The RNA-seq data used and analyzed during the current study are included in this published article and its supplementary files. The open RNA-seq data (Accession number: GSE155408) is available in the NCBI database (https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE15 5408).