Trang TT. Truong1*
Chiara C. Bortolasci1Briana Spolding1Bruna Panizzutti1Zoe SJ. Liu1Srisaiyini Kidnapillai1
Mark Richardson2
Laura Gray1
Craig M. Smith1
Olivia M. Dean1,3
Jee Hyun Kim1,3Michael Berk1,3,4Ken Walder1- 1School of Medicine, IMPACT, Institute for Innovation in Physical and Mental health and Clinical Translation, Deakin University, Geelong, VIC, Australia
- 2Genomics Centre, School of Life and Environmental Sciences, Deakin University, Burwood, VIC, Australia
- 3The Florey Institute of Neuroscience and Mental Health, University of Melbourne, Parkville, VIC, Australia
- 4Orygen, The National Centre of Excellence in Youth Mental Health, Centre for Youth Mental Health, Florey Institute for Neuroscience and Mental Health and the Department of Psychiatry, The University of Melbourne, Melbourne, VIC, Australia
Long non-coding RNAs (lncRNAs) may play a role in psychiatric diseases including bipolar disorder (BD). We investigated mRNA-lncRNA co-expression patterns in neuronal-like cells treated with widely prescribed BD medications. The aim was to unveil insights into the complex mechanisms of BD medications and highlight potential targets for new drug development. Human neuronal-like (NT2-N) cells were treated with either lamotrigine, lithium, quetiapine, valproate or vehicle for 24 h. Genome-wide mRNA expression was quantified for weighted gene co-expression network analysis (WGCNA) to correlate the expression levels of mRNAs with lncRNAs. Functional enrichment analysis and hub lncRNA identification was conducted on key co-expressed modules associated with the drug response. We constructed lncRNA-mRNA co-expression networks and identified key modules underlying these treatments, as well as their enriched biological functions. Processes enriched in key modules included synaptic vesicle cycle, endoplasmic reticulum-related functions and neurodevelopment. Several lncRNAs such as GAS6-AS1 and MIR100HG were highlighted as driver genes of key modules. Our study demonstrates the key role of lncRNAs in the mechanism(s) of action of BD drugs. Several lncRNAs have been suggested as major regulators of medication effects and are worthy of further investigation as novel drug targets to treat BD.
Introduction
Bipolar disorder (BD) is among the top 10 causes of disability globally with potentially devastating consequences for individuals as well as wider society (Association, 2013; James et al., 2018). Although pharmacotherapy is the first-line management for BD, suboptimal outcomes and treatment-resistance are common, and polypharmacy is often required (Garnham et al., 2007; Kukreja et al., 2013; Fornaro et al., 2016). Despite the pressing demand for new therapeutic agents, new drug discovery for BD remains stagnant due to our lack of understanding of both the underlying pathophysiology of the disorder and the mechanism(s) of action of currently available drugs.
Studies have suggested a role for dysregulation of gene expression in the pathophysiology of BD, not only with protein-coding but also with non-coding RNAs (Ponting et al., 2009; Zuo et al., 2016; Luykx et al., 2019). Most of the non-coding transcriptome is made up of long non-coding RNAs (lncRNAs). These comprise >200 nucleotides, and are key regulators of gene expression possibly via epigenetic regulation and chromatin remodelling (Hrdlickova et al., 2014). LncRNAs are highly expressed in the brain, and they may contribute to the development of psychiatric diseases including BD (Qureshi et al., 2010; Zuo et al., 2016; Bella and Campo, 2021). In a genome-wide association study (GWAS), Hou et al. found a significant association between lithium response and four single nucleotide polymorphisms (SNPs) in a region containing two lncRNAs, suggesting a role for these lncRNAs in the biological mechanism(s) of lithium (Hou et al., 2016). Given that little is known regarding how lncRNAs may be involved in the mechanism(s) of action of drugs used to treat BD, comprehensive gene expression analysis that includes lncRNAs will help to elucidate the transcriptional perturbations underlying the therapeutic effects.
When researching brain disorders, reductionist approaches are still dominant: individual candidate receptors or genes are often analysed singly, disregarding their interaction with other molecular factors. To mitigate such limitation, network-based approaches that consider modules of genes as key regulators offer a complementary method to address the complex dysregulation occurring in neuropsychiatric conditions such as BD (Gaiteri et al., 2014). Moreover, it remains unclear how most genes function, which has become more puzzling with the recent identification of novel non-coding genes. A detailed global view of genome-wide transcriptional perturbations can be yielded with the advent of high-throughput sequencing technologies. Genes with functional linkages tend to be co-expressed across various biological states, and co-expression networks yielded from RNA-sequencing data can infer gene function and molecular mechanisms by associating genes of unknown function with biological processes (Gaiteri et al., 2014; Okamura et al., 2015).
In the current study, we investigated the mechanism(s) of action of four widely prescribed BD drugs (i.e., lamotrigine, lithium, quetiapine and valproate) using network-based approaches, with an emphasis on the role of lncRNAs. We hypothesised that mRNA-lncRNA crosstalk and regulatory patterns in neuronal cell models will unveil insights into the complex mechanisms of BD drugs and highlight potential new targets for drug development.
Materials and Methods
The summary flowchart of the methodology was demonstrated in Figure 1.
FIGURE 1. The methodology workflow the current study.
Cell Culture
For the model of human neurons, the study used NT2 human teratocarcinoma cells (CVCL_0034, ATCC, Manassas, VI, United States), differentiated into post-mitotic neuronal cells (NT2-N) after treatment with retinoic acid (Pleasure et al., 1992; Megiorni et al., 2005). Cells were cultured and differentiated as previously described (Bortolasci et al., 2018). Concisely, NT2 cells were maintained in Dulbecco’s modified Eagle’s Medium (DMEM; Life Technologies, Melbourne, Australia), 10% fetal bovine serum (FBS; Thermo Fisher Scientific, Melbourne, Australia) and 1% antibiotic-antimycotic solution (Life Technologies). NT2-N cells were induced from NT2 cell cultures by treating with 10−5 M retinoic acid (Sigma-Aldrich, Sydney, Australia) for 28 days with media refreshed every 2–3 days. For experiments, cells were seeded onto 24-well plates coated with 10 μg/ml poly-d-lysine (Sigma-Aldrich) and 10 μg/ml laminin (Sigma-Aldrich) at 2 × 105 cells/well with further addition of mitotic inhibitors (1 µM cytosine and 10 µM uridine; Sigma-Aldrich) for a total of 7 days, and the media was refreshed every 2–3 days to generate an enriched culture of differentiated neuronal cells (NT2-N). NeuroD, Tau, and GluR expression levels in differentiated cells were measured to verify the effectiveness of the differentiation process (data not shown).
Drug Treatments
NT2-N cells were treated with lamotrigine (50 µM), lithium (2.5 mM), quetiapine (50 µM), or valproate (0.5 mM) for 24 h (4 replicates for each group). All drugs were purchased from Sigma-Aldrich (Sydney, Australia). Vehicle control cells were treated with 0.5% Milli-Q water for lithium or valproate controls, and 0.2% DMSO for lamotrigine or quetiapine controls. The drug doses were chosen based on prior dose response studies in our lab to ensure no impact on cell viability (Kidnapillai et al., 2018). Using these doses, we have previously observed differences in gene expression after 24 h treatment (Bortolasci et al., 2018; Bortolasci et al., 2020b). These drugs have also been demonstrated to alter gene expression levels in vitro after 24 h in other studies (Nahman et al., 2012; Meng et al., 2019).
Genome-Wide Gene Expression Measurement
Total RNA was extracted from cells post-treatment using RNeasy∗ mini kits (Qiagen, Melbourne, Australia), then checked for quality and quantity using an Agilent 2100 Bioanalyzer (Agilent Technologies, Melbourne, Australia) and NanoDrop 1000 (Thermo Fisher Scientific), respectively. RNA-seq libraries were prepared for all samples from 1 µg total RNA using a TruSeq RNA Sample Preparation Kit (Illumina, Victoria, Australia) as per the manufacturer’s instructions. Samples were analysed on a HiSeq 2500 Rapid system (50 bp single end reads; Illumina) to measure genome wide mRNA expression, which yielded an average of about 10 × 106 reads/sample. The single end sequencing offered sufficient quantification of gene expression levels in well-annotated organisms such as Homo sapiens at cheaper cost than the paired-end technology (Conesa et al., 2016).
Genome-Wide Gene Expression Analysis
The raw data yielded in FASTQ format were aligned to reference genomes using the Deakin Genomics Centre RNA-Seq alignment and expression quantification pipeline (https://github.com/m-richardson/RNASeq_pipe). Briefly, the pipeline included raw read quality filtering and adapter trimming (ILLUMINACLIP:2:30:10:4, SLIDINGWINDOW:5:20, AVGQUAL:20 MINLEN:36) with Trimmomatic v35 (Bolger et al., 2014), and alignment to the reference genome using STAR v2.5 in 2-pass mode (Human genome version GRCh38) (Dobin et al., 2012). Raw reads were deposited at the Gene Expression Omnibus (GEO) database under the accession numbers GSE197966.
For differential abundance testing, gene expression quantified at the gene level was compiled into a m x n matrix from individual sample counts. Genes with low expression were omitted (<1 cpm in n samples, where n is the number of samples in the smallest group for comparison), and the data was normalized using the weighted trimmed mean of M-values (TMM) using edgeR (Robinson et al., 2009) in R (R Core Team, 2021). EdgeR was used to determine differential gene expression, and the Benjamini–Hochberg method was used to assess false-discovery rate (FDR) adjusted q-values for multiple tests (Benjamini and Hochberg, 1995).
Weighted Gene Co-Expression Network Analysis
WGCNA is an approach utilizing gene expression data to construct co-expression networks weighted for high correlations (Langfelder and Horvath, 2008) and was used in this study to evaluate correlation between lncRNAs and mRNAs. The RNA-seq data was used as input for the R package WGCNA (Langfelder and Horvath, 2008; R Core Team, 2021), from which a pairwise bi-weight mid-correlation matrix was computed and then transformed into an adjacency matrix. To construct a scale-free network, each absolute mid-correlation value was raised to a soft-thresholding power. Soft-thresholding amplified disparity between strong and weak correlations, leading to the construction of the scale-free network. We chose power 7, which was the lowest power for which the scale-free topology fit index reached 0.8. To minimise the spurious connections, WGCNA utilised the topological overlap measure (TOM) accounting for how large the overlap of each gene pair’s network neighbours. The TOM matrix considered as a similarity measure was then transformed into pairwise dissimilarity measure (calculated by 1-TOM) for the hierarchical clustering of genes. From this, tightly connected genes would be clustered for module assignment (dynamic tree cut algorithm), and unassigned genes with weak connections would not be considered for further analyses. The default value of 0.25 was set as the threshold for cut height to merge possible similar modules. The expression level of each module was represented by an eigengene value. Module membership values were also calculated, which reflect the degree of correlation between genes and modules. Higher absolute values of module memberships mean stronger correlation, while zero values mean no association.
Module Preservation and Module-Trait Analyses
To validate the reproducibility of key modules, module preservation analysis was undertaken with an independent dataset (Bortolasci et al., 2020a). The most common measure of module preservation recommended by the authors of WGCNA is Zsummary, which considers the density and connectivity patterns of module nodes, as well as the overlap among module membership. Higher Zsummary statistics indicate more highly preserved modules. Notably, Zsummary can be heavily dependent on module sizes, which might bias the statistics towards bigger modules (Langfelder et al., 2011). In our network, module sizes were spread in a wide range from 30 to more than 3000, hence we only used a cut off Zsummary >10 as recommended for strong preservation and then ranked the degree of preservation based on medianRank. The lower medianRank one module has, the stronger preservation it tends to exhibit in another dataset. Langfelder et al. found medianRank to be much less dependent on module sizes, and hence it can be a useful measure for comparing relative preservation among multiple modules (Langfelder et al., 2011).
To identify key modules related to the response to drug treatments, association analysis between a module and the trait of each pairwise comparison group (e.g., lithium–water, quetiapine–DMSO) was performed based on module eigengenes. Student asymptotic p-values were calculated for correlation values, and then adjusted for multiple testing with Benjamini–Hochberg false discovery rate correction.
Functional Enrichment Analysis
In order to find potential pathways that might be driven by key modules, functional enrichment analysis was deployed on all genes from each module using the R package ClusterProfiler (Yu et al., 2012) with pathway reference from the Gene Ontology (GO) database (Gene Ontology, 2004) filtered by “Biological Process”. ClusterProfiler is a popular enrichment package, which was extensively utilized in medical studies (Liu et al., 2021; Liu et al., 2022a; Liu et al., 2022b).
Despite their comprehensiveness, GO terms are prone to redundancy, which might complicate the functional interpretation. Hence, we deployed the Enrichment Map from Cytoscape (Merico et al., 2010) to intuitively visualise the enrichment results in a network-based manner. Enrichment Map clusters similar gene sets together (based on how well-connected their genes are), offering a more concise overview of enriched biological functions. In the map, gene sets were illustrated as nodes and would be linked together if they share overlapping genes. The overlapping metric was calculated based on a combination of Jaccard coefficient and overlap coefficient (50% and 50% respectively), and a cut-off of 0.375 was set to define edges forming between nodes.
lncRNA-mRNA Co-Expression Networks, Identification of Hub lncRNAs, Cis- and Trans-Interactions
The lncRNA-mRNA co-expression network of key modules associated with the drug treatments was built using the bi-weight mid-correlation of lncRNA-mRNA pairs. The lncRNA-mRNA pairs whose weights of connection were at least 0.15 were selected for co-expression network construction. The network was imported and visualized in Cytoscape (Shannon et al., 2003). To find potential key regulators of each module, hub lncRNAs were identified by degree of centrality levels calculated using Cytoscape (Shannon et al., 2003).
To find potential cis-targets of lncRNAs in each key module, nearest protein coding genes located 100 kb upstream or downstream from the transcription start site of lncRNAs were identified by BEDTools (Quinlan and Hall, 2010) v2.27 with annotation from human genome version GRCh38. Otherwise, the other mRNAs in the module would be considered to have trans-interaction with lncRNAs in the module. Among these trans-interactions, we evaluated if any of them could be due to the regulation of transcription factors on lncRNAs as co-expression relationships can be interpreted in both directions. Hence, we searched for trans-acting transcription factors found nearby lncRNAs (from −30 to 10 kb away from transcription start sites) using chromatin immunoprecipitation followed by sequencing (ChIP-seq) data from the ChIPBase database (Zhou et al., 2017).
Results
Differential Expression Analysis
The characteristics of RNA-seq results for each drug treatment and differential analysis mainly for mRNAs have been reported in our previous publication (Bortolasci et al., 2020b). For lncRNAs, we found 1044 were expressed at detectable levels in NT2-N cells. Among them, 282 lncRNAs were differentially expressed (FDR <0.05) following treatment of NT2-N cells with one or more drugs used to treat BD (Supplementary Table S1). While differential expression analysis is the initial approach for evaluating genes driving the difference between different phenotypes (in our case, treatments and controls), it can be challenging to study lncRNAs which tend to be expressed at lower levels than mRNAs (Ulitsky, 2016; 2019). The significance of p-values was evaluated based on the whole distribution of all genes detected in our RNA-seq (including mRNAs, lncRNAs), thus any genes with lower fold change would be less likely to be significantly differentially expressed after correction for multiple testing. Therefore, differential expression analysis was complemented with an unsupervised gene co-expression analysis method to discriminate small but coherent patterns of differences.
Weighted Gene Co-Expression Network Analysis
After identifying individual lncRNAs that were differentially expressed, we next sought to identify co-expressed networks of genes from the whole transcriptome, as this approach may reveal additional lncRNAs of interest and other important mechanistic insights. Without any a priori defined groups like in differential expression analysis, weighted gene co-expression network analysis (WGCNA) on the whole transcriptome explored highly correlated and consistent transcriptional patterns across 24 samples of our RNA-seq results (four samples x two vehicles, four samples x four drug treatments). While bigger sample size would be more ideal, 24 samples are within the recommendation of the WGCNA workflow for robust analyses (Langfelder and Horvath, 2022). WGCNA identified 29 modules and the full list of genes in each module is shown in Supplementary Table S2.
To validate the reproducibility of the modules in the current WGCNA network, we performed preservation analysis of topology and connectivity patterns against an independent RNA-seq dataset from our previous study that used different treatments of psychotropic drugs on the same NT2-N cell model (Bortolasci et al., 2020b). Biologically interesting modules are expected to be preserved in other samples with a homogeneous cellular population, which implies sets of genes working in concert for specific biological functions. From the results of the module preservation analysis (illustrated in Supplementary Figure S1), any module with a z-summary score higher than 10 was regarded as strongly preserved, meaning the co-expression relationships of its member genes were maintained across the two datasets. Among the 11 modules found to be strongly preserved in this independent dataset, module five was the most highly preserved module.
Identification of Key Modules Associated With Drug Treatment (Module-Trait Analysis)
Results from the module-trait association analysis for the 11 preserved modules are shown in Figure 2 (each cell in the heatmap has a correlation value and an adjusted p-value). The acquired module-trait relationships enabled the evaluation of which module(s) were most affected by each drug. While the most significant module associated with each drug treatment was of interest, we cannot rule out the possibility that such an association could be accounted for by off-target effects. Hence, we further filtered modules which were coherently associated with at least three drug-vehicle pairs in a similar manner with statistical significance of adjusted p-value < 0.05 (correlation values can be either negative or positive, but they must share the same directionality in at least three traits), as similar effects of multiple drug treatments on one module tend to imply the common therapeutic mechanisms of the drugs rather than distinct side effects of each drug. We identified seven modules meeting such requirements (modules 1, 3, 4, 5, 7, 9, and 10—labelled with a hash symbol in Figure 2) and these modules were regarded as key modules associated with BD drug treatments. In further analyses, we focused on the key modules to find co-expressed relationships of genes belonging to each module, and which drug treatment was most strongly associated with each module.
FIGURE 2. Heatmap of module-trait relationships with corresponding p-values between the detected modules on the y-axis and vehicle-drug traits on the x-axis. The correlation value and Benjamini–Hochberg adjusted p-value for each pair are labelled on the cell in the format: correlation (p-value). Each cell of the heatmap is coloured based on correlation between each module eigenvalue and the trait: blue is a strong positive correlation, red is a strong negative correlation, and white is little to no correlation. For example, regarding the DMSO-Lamotrigine trait, the module 5, with the negative correlation value of −0.66 and significant adjusted p-value of 0.002, tends to have lower overall expression (summarised as eigengene value) in lamotrigine treatment compared with the corresponding DMSO vehicle control. Abbreviation: DMSO_LAM–Lamotrigine treatment versus DMSO vehicle, DMSO_QUE–Quetiapine treatment versus DMSO vehicle, H2O_LIT–Lithium treatment versus water vehicle, H2O_VAL–Valproate treatment versus water vehicle. Modules labelled with hash (#) were coherently regulated with same directionality by at least three drugs. LncRNA-mRNA co-expression networks, identification of hub lncRNAs.
The lncRNA-mRNA co-expression networks were constructed on the six key modules coherently regulated with the same directionality by three or more drugs to shed light on the molecular mechanism of lncRNAs that might be potential targets of drugs used to treat BD. Module nine was eliminated from this step due to the lack of lncRNA presence in the module. Hence, there were five key modules associated with the effects of BD drug treatments (i.e., modules 1, 4, 5, 7, and 10) being constructed for the lncRNA-mRNA co-expression networks. The integrated network of these modules is illustrated in Supplementary Figure S2 and separated sub-networks for single modules are shown in Figure 3.
FIGURE 3. Separated lncRNA-mRNA subnetworks of five key modules. In the network, nodes representing genes and edges representing co-expression connections between them, mRNAs are circle nodes while lncRNAs are diamond shaped. Each subnetwork corresponds to a module: (A) Module 1 (red nodes), (B) Module 4 (pink nodes), (C) Module 5 (orange nodes), (D) Module 7 (yellow nodes), and (E) Module 10 (blue nodes).
Network statistics of lncRNAs identified from these modules are presented in Supplementary Table S3. Genes from each module were ranked by degree of centrality, which measures the number of edges (interactions) each node has. Nodes with greater importance tend to lie on multiple paths (edges) between other nodes and hence finding hub genes with high degrees of centrality can illuminate the greatest influencers of the biological networks (Mason and Verwoerd, 2007). Interestingly, the most connected lncRNAs (hub lncRNAs) identified tended to be differentially expressed in the drug treatments having highest module-trait relationships with these lncRNAs’ modules. The complementary results of differential expression analysis reinforced the role of these lncRNAs in their corresponding module: the eigen value representing a module most correlated with a certain treatment (module-trait relationship in co-expression network analysis), and hence the hub genes that have the greatest contribution to such module regulation tended to have more significant difference induced by such treatment (reflected as log fold change) than other treatments, to an extent that reached the FDR 0.05 cut-off in differential expression analysis for single genes.
Cis- and Trans-Interactions
To explore the characteristics of regulatory relationships between lncRNAs and mRNAs, we searched for cis-targets within each module based on proximity of coding genes to transcription start sites of lncRNAs, as well as genes encoding transcription factors trans-interacting with lncRNAs based on known binding evidence from chromatin immunoprecipitation followed by sequencing (ChIP-seq) data (Supplementary Table S4). In the five lncRNA-mRNA key modules, SETD1B identified as potential cis-targets identified within 100 kb of LINC01089. We identified four lncRNAs that might potentially be regulated themselves by coding genes in the same module that are trans-acting transcription factors found nearby the promoter or enhancer region of those lncRNAs. These transcription factor-lncRNA regulation relationships were found only in several lncRNAs located in two modules, and most trans-interactions of lncRNAs-mRNAs tended to be regulated by lncRNAs rather than the other way around.
Functional Enrichment Analysis
Over-representation analysis was undertaken for genes from five key modules with hub lncRNAs prominently contributing to the module regulation (i.e., modules 1, 4, 5, 7, and 10) to identify the main biological processes driven by each module (detailed results provided in Supplementary Table S5). GO terms from each module with p-values <0.001 were then analysed in Enrichment Map to be grouped into clusters, which highlights the major functional themes enriched in the five key modules (Figure 4).
FIGURE 4. Enrichment Map of five key modules with hub lncRNAs contributing to the modular regulation induced by bipolar disorder drugs. Enriched gene sets are represented as nodes, while edges connect the similar gene sets together. Nodes are coloured by the module they enriched with the coloured patches covering all enriched gene sets by module; one node might be enriched in multiple modules. The bigger the node, the higher the number of genes found in the gene set. The thickness of each edge is proportional to the number of mutual genes between nodes. Gene sets with similar functions are clustered and labelled based on the main theme they belonged to. Abbreviation: ER–endoplasmic reticulum.
The Enrichment Map showed a clear separation of biological functions enriched by the five key modules. Module one appears to be related to epigenetic modification via the enrichment of genes regulating histone acetylation. Module four showed enrichment in cell morphogenesis, RNA metabolism and response to oxidative stress. Module five focused on vesicle docking and fusion of trans-synaptic transmission. Module seven was enriched for peptidyl tyrosine phosphorylation and protein catabolism, which might be involved in the response to endoplasmic reticulum stress. Module 10 may regulate the formation of Purkinje cells and the basement membrane organization.
Discussion
In the current study, multi-stepped approaches using systems biology methodology were deployed to investigate the molecular mechanisms of commonly prescribed BD drugs with a focus on lncRNAs and their potential regulatory connections with mRNAs. LncRNAs are generally expressed at lower levels than protein coding RNAs, but the cascade of regulations induced by lncRNAs may amplify the effects they produce (Ulitsky, 2016; Hu et al., 2018; Encode, 2019). However, they might be overlooked by conventional pair-wise gene expression comparisons, in which the extent of differential expression highlights significant genes for further analyses. Despite their low abundance that could be considered barely changed using conventional differential expression methods, lncRNAs can exert pronounced effects via their distinct mechanisms, e.g., one to ten molecules per cell is likely adequate for enabling cis-acting lncRNAs to affect transcription at a single locus or at several loci through either direct base-pairing with genomic DNA or recruitment of chromatin modifiers (Wu et al., 2021). Several lncRNAs were found to drive vital processes even with subtle changes such as RepA (X-chromosome inactivation) (Zhao et al., 2008), and VELUCT (lung cancer cellular viability) (Seiler et al., 2017). Since genes rarely act alone, considering them in the context of biological networks that they shared interactions with is fundamental to gain better understanding of regulatory entities such as transcription factors and lncRNAs.
Network-based analyses such as co-expression networks prioritised relationships (co-expression) rather than just the states of singular components (mean expression changes). An analysis based only on mean changes in expression could lead to an incorrect conclusion about a particular pathway’s involvement in a condition given that genes may change their partners depending on dynamic biological demands to recruit or disassociate groups of co-regulated genes for a particular task (Mentzen et al., 2009). Pathways with significant mean expression alteration but decreased co-expression could imply a change in functional assignment but do not guarantee whether genes a priori assumed to be in the certain biological pathways are actually dedicated to such pathways since genes could have switched roles and interacted with genes from different pathways (Mentzen et al., 2009; de la Fuente, 2010). The essential need to consider co-expression changes in addition to differential mean expression when comparing gene expression datasets has been emphasised in multiple biological contexts (Mentzen et al., 2009; de la Fuente, 2010; Savino et al., 2020; Weighill et al., 2021). For example, co-expression analyses of cancer gene expression datasets found several transcription factors known to regulate cancer development were identified as highly differentially co-expressed, despite their mean expression levels not having significant changes (de la Fuente, 2010). A proof-of-concept example from Hudson et al. showed co-expression analysis correctly identified a causal gene with a mutation between two bull varieties, whereas the gene had non-significant mean expression changes (Hudson et al., 2009).
In our analyses, co-expression network analysis utilised the whole transcriptome to offer insights on the role of lncRNAs via their connection to other biological entities such as mRNAs (Parikshak et al., 2015). Hence, such a systems biology approach is more ideal to characterise the complexity of BD as well as its treatments (Gaiteri et al., 2014). While WGCNA has been deployed in some psychiatric disorders including BD (Akula et al., 2016; Fromer et al., 2016; Kim et al., 2016; Liu et al., 2019; Zhang et al., 2021), to our knowledge, this is the first-time co-expression network analysis was used for exploring mechanism(s) of action of drugs used to treat BD.
Based on RNA-seq gene expression profiles of four BD drug treatments in NT2-N cells, we applied WGCNA to explore co-expressed genes with a focus on lncRNA-mRNA connectivity and to further identify hub lncRNAs associated with the drug treatments. In many cases, hub genes, particularly intramodular hubs, play a greater role than other network nodes in determining the network’s functionality (Horvath, 2011). Hence, identification of hub lncRNAs might illuminate major lncRNAs influencing the transcriptional regulation of their corresponding modules, which were strongly associated to a certain drug treatment. We also evaluated the type of regulatory interactions between lncRNA and mRNAs, such as potential cis mRNA targets, and trans-acting transcription factor-lncRNA relationships based on current knowledge of transcription factor binding sites. As lncRNAs might be cis- or trans-acting, while the prior can be identified based on the proximity of coding genes to the lncRNAs, pinpointing the latter is challenging due to the complex mechanisms (Statello et al., 2021). Therefore, we identified the potential transcription factors that might regulate lncRNAs, leaving the rest of lncRNA-mRNA trans-interactions more likely to be regulated by lncRNAs.
From 29 modules clustered upon the co-expression pattern of the whole dataset, we identified five key modules for lncRNA-mRNA networks as potential drivers of BD drug treatments. As the roles of lncRNAs are mainly unknown in BD, their functions can be speculated via the connected mRNAs in the modules as per the ‘guilt-by-association’ principle (van Dam et al., 2018). Such association was employed based on how genes co-expressed with one and another across all samples, with the justification that closely regulated genes are more likely to be associated with similar functions. Interestingly, our enrichment analysis showed while there was separation in functional annotations between modules, the enriched functions tended to be complementary, serving several major biological processes such as synaptic vesicle cycle, actin filament organization, endoplasmic reticulum-related functions and neurodevelopment. Findings from previous studies using co-expression network analysis for bipolar disorder have found similar processes significantly enriched by hub genes, e.g., regulation of transcription, postsynaptic density, ribosomal subunit, endocytosis (Akula et al., 2016; Liu et al., 2019), actin filament-based process, axon development (Zhang et al., 2021). Co-expression analysis on RNA-seq of whole blood from BD patients found the co-expressed modules associated to lithium usage enriched endoplasmic reticulum related functions (Krebs et al., 2020). However, it should be acknowledged that these studies mainly focused on mRNAs rather than lncRNAs for their co-expression network construction, hence the involvement of lncRNAs was not evaluated. Nevertheless, the similar enriched biological processes in BD-related phenotypes supported the findings found in our current study.
Our functional enrichment analysis on key modules highlighted their association to the synaptic vesicle cycle, in which module five was enriched for vesicle docking and fusion processes. The synaptic vesicle cycle is highly relevant to BD given that the classic pathophysiological hypothesis has been built upon dysregulation of monoamine transmission (Stahl, 2013). The expression of the hub lncRNA of module 5, GAS6-AS1, was upregulated by valproate–the treatment that module five correlated to the most. While the role of GAS6-AS1 in BD is unclear, it promoted cell proliferation, migration and invasion in several cancer cell lines (Zhang et al., 2019; Li et al., 2020). Interestingly, GAS6-AS1 was also shown to activate the PI3K/AKT pathway (Li et al., 2020), which itself plays a vital role in vesicle trafficking (Bhattacharya et al., 2016; Bilanges et al., 2019).
Targeting the endoplasmic reticulum (ER) and its related processes could be one major function regulated by BD drugs, as inferred from multiple modules in this study: response to ER stress mainly in module seven and partially in module 5. When ER stress occurs, autophagy is activated to restore cellular homeostasis, but this response is compromised in BD (Susanne et al., 2016). There is converging evidence demonstrating that lithium increases autophagy via inositol depletion (Bar-Yosef et al., 2019). An analysis of the lithium response gene network in BD-patient derived lymphoblastoid cell lines also identified ER stress as a major module (Breen et al., 2016). LINC02381 was identified as the hub lncRNA in module 7. LINC02381 has not been researched in BD previously. Evidence from gastric cancer cell lines suggested LINC02381 can reduce Wnt pathway activity and increase apoptosis (Jafarzadeh and Soltani, 2020). Our previous study presented the downregulation of the Wnt pathway in multiple antipsychotic treatments including quetiapine, suggesting this pathway might be a common mechanism induced by different antipsychotics (Panizzutti et al., 2021). In the current study, such an effect could be partially explained by the upregulation of hub lncRNA LINC02381 by quetiapine shown in differential expression analysis, which in turn attenuates the Wnt pathway.
Module 10 appears to contain an integrated cluster of functions related to neurodevelopment via its effect on Purkinje cells and the basement membrane. The basement membrane is part of the extracellular matrix system that plays a critical role in corticogenesis involving Purkinje cells in the cerebral cortex (Franco and Müller, 2011). Damage-associated molecular patterns of extracellular matrix components of the basement membrane induce an immune response suggested to be part of the pathogenesis of BD (Rege and Hodgkinson, 2013). The expression of two hub lncRNAs of module 10, MIR124-2HG and MIR100HG, were found to be differentially expressed by several drug treatments. MIR124-2HG is the host gene of miR-124, which is the most abundant brain-specific miRNA regulating neuronal differentiation during CNS development and adult neurogenesis (Sonntag et al., 2012). Moreover, miR-124 was shown to have a pathophysiological contribution in some neuropsychiatric disorders such as Alzheimer’s disease and autism (Tabarés-Seisdedos and Rubenstein, 2009; Fang et al., 2012). MIR100HG (alternative name Linc-NeD125) was suggested in human neuroblastoma-derived cells to be specifically induced during neuronal differentiation to support cell survival (Bevilacqua et al., 2015). While MIR100HG encodes three miRNAs in its intron (i.e., miR-100, miR-125b and let-7a) (Ottaviani et al., 2018), it might work independently of the hosted microRNAs to promote the required conditions for differentiation by attenuating cell proliferation and exerting its anti-apoptotic function via the activation of BCL-2 (Bevilacqua et al., 2015). These hub genes have yet to be explained in terms of their biological functions in BD, hence more studies are needed to illuminate their involvement in the disorder.
This study is not without limitations. We used the in vitro model of NT2-N cells that have limited capacity to represent the disease state. In addition, the 24-h administration with a single dose for each drug for our NGS analyses limits the overall evaluation on the long-term pharmacological regulation. Finally, the analysis was focused on one dataset due to the limited availability of high-throughput expression data of BD drug responses. More representative results could be obtained from the incorporation of multiple datasets, as well as the addition of other ncRNA species such as miRNAs and siRNAs to obtain more comprehensive insights on the molecular mechanism(s) of BD drugs.
In conclusion, our study demonstrates the potential key role of lncRNAs in the regulatory effect of BD drugs via the associated lncRNA-mRNA co-expression networks. Several major processes were enriched in key modules associated with drug treatments such as synaptic vesicle cycle, cell cycle, endoplasmic reticulum-related functions and neurodevelopment. These results contribute to our understanding of the mechanisms of action of BD drugs and suggest potential novel targets for therapeutic intervention.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Center for Biotechnology Information (NCBI) BioProject database under accession number GSE197966.
Author Contributions
Conceptualisation, KW and MB; methodology, TT, CB, BS, OD, and KW; formal analysis, TT, KW, and JK; investigation, TT, CB, SK, BS, BP, and LG; resources, KW; data curation, MR and CS; writing—original draft preparation, TT; writing—review and editing, CB, SK, BS, BP, ZL, LG, JK, CS, OD, MR, MB, and KW; visualisation: TT, ZL, and JK; funding acquisition, KW.
Funding
This research was funded by National Health and Medical Research Council (NHMRC) Project Grant (1078928) and Centre of Research Excellence (1153607). MB is supported by a NHMRC Senior Principal Research Fellowship (1156072). OMD is an R.D. Wright Biomedical NHMRC Career Development Fellow (1145634) and has received grant support from the Brain and Behavior Foundation, Marion and EH Flack Trust, Simons Autism Foundation, Australian Rotary Health, Stanley Medical Research Institute, Deakin University, Brazilian Society Mobility Program, Lilly, NHMRC, Australasian Society for Bipolar and Depressive Disorders and Sevier. She has also received in-kind support from BioMedica Nutracuticals, NutritionCare and Bioceuticals.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.873271/full#supplementary-material
References
Akula, N., Wendland, J. R., Choi, K. H., and McMahon, F. J. (2016). An Integrative Genomic Study Implicates the Postsynaptic Density in the Pathogenesis of Bipolar Disorder. Neuropsychopharmacology 41 (3), 886–895. doi:10.1038/npp.2015.218
Association, A. P. (2013). Diagnostic and Statistical Manual of Mental Disorders (DSM-5®). American Psychiatric Pub.
Bar-Yosef, T., Damri, O., and Agam, G. (2019). Dual Role of Autophagy in Diseases of the Central Nervous System. Front Cell Neurosci 13, 196. doi:10.3389/fncel.2019.00196
Bella, F., and Campo, S. (2021). Long Non-coding RNAs and Their Involvement in Bipolar Disorders. Gene 796-797, 145803. doi:10.1016/j.gene.2021.145803
Benjamini, Y., and Hochberg, Y. (1995). Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. Ser. B (Methodological) 57 (1), 289–300. doi:10.1111/j.2517-6161.1995.tb02031.x
Bevilacqua, V., Gioia, U., Di Carlo, V., Tortorelli, A. F., Colombo, T., Bozzoni, I., et al. (2015). Identification of Linc-NeD125, a Novel Long Non Coding RNA that Hosts miR-125b-1 and Negatively Controls Proliferation of Human Neuroblastoma Cells. RNA Biol. 12 (12), 1323–1337. doi:10.1080/15476286.2015.1096488
Bhattacharya, S., McElhanon, K. E., Gushchina, L. V., and Weisleder, N. (2016). Role of Phosphatidylinositol-4,5-Bisphosphate 3-kinase Signaling in Vesicular Trafficking. Life Sci. 167, 39–45. doi:10.1016/j.lfs.2016.10.018
Bilanges, B., Posor, Y., and Vanhaesebroeck, B. (2019). PI3K Isoforms in Cell Signalling and Vesicle Trafficking. Nat. Rev. Mol. Cell Biol 20 (9), 515–534. doi:10.1038/s41580-019-0129-z
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a Flexible Trimmer for Illumina Sequence Data. Bioinformatics 30 (15), 2114–2120. doi:10.1093/bioinformatics/btu170
Bortolasci, C. C., Spolding, B., Callaly, E., Martin, S., Panizzutti, B., Kidnapillai, S., et al. (2018). Mechanisms Underpinning the Polypharmacy Effects of Medications in Psychiatry. Int. J. Neuropsychopharmacol. 21 (6), 582–591. doi:10.1093/ijnp/pyy014
Bortolasci, C. C., Spolding, B., Kidnapillai, S., Connor, T., Truong, T. T. T., Liu, Z. S. J., et al. (2020a). Transcriptional Effects of Psychoactive Drugs on Genes Involved in Neurogenesis. Int. J. Mol. Sci. 21 (21). doi:10.3390/ijms21218333
Bortolasci, C. C., Spolding, B., Kidnapillai, S., Richardson, M. F., Vasilijevic, N., Martin, S. D., et al. (2020b). Effects of Psychoactive Drugs on Cellular Bioenergetic Pathways. World J. Biol. Psychiatry 22, 1–15. doi:10.1080/15622975.2020.1755450
Breen, M. S., White, C. H., Shekhtman, T., Lin, K., Looney, D., Woelk, C. H., et al. (2016). Lithium-responsive Genes and Gene Networks in Bipolar Disorder Patient-Derived Lymphoblastoid Cell Lines. Pharmacogenomics J. 16 (5), 446–453. doi:10.1038/tpj.2016.50
Conesa, A., Madrigal, P., Tarazona, S., Gomez-Cabrero, D., Cervera, A., McPherson, A., et al. (2016). A Survey of Best Practices for RNA-Seq Data Analysis. Genome Biol. 17 (1), 13. doi:10.1186/s13059-016-0881-8
de la Fuente, A. (2010). From ‘differential Expression’ to ‘differential Networking’ – Identification of Dysfunctional Regulatory Networks in Diseases. Trends Genet. 26 (7), 326–333. doi:10.1016/j.tig.2010.05.001
Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2012). STAR: Ultrafast Universal RNA-Seq Aligner. Bioinformatics 29 (1), 15–21. doi:10.1093/bioinformatics/bts635
Fang, M., Wang, J., Zhang, X., Geng, Y., Hu, Z., Rudd, J. A., et al. (2012). The miR-124 Regulates the Expression of BACE1/β-Secretase Correlated with Cell Death in Alzheimer's Disease. Toxicol. Lett. 209 (1), 94–105. doi:10.1016/j.toxlet.2011.11.032
Fornaro, M., De Berardis, D., Koshy, A. S., Perna, G., Valchera, A., Vancampfort, D., et al. (2016). Prevalence and Clinical Features Associated with Bipolar Disorder Polypharmacy: a Systematic Review. Neuropsychiatr. Dis. Treat. 12, 719–735. doi:10.2147/NDT.S100846
Franco, S. J., and Müller, U. (2011). Extracellular Matrix Functions during Neuronal Migration and Lamination in the Mammalian central Nervous System. Dev. Neurobiol. 71 (11), 889–900. doi:10.1002/dneu.20946
Fromer, M., Roussos, P., Sieberts, S. K., Johnson, J. S., Kavanagh, D. H., Perumal, T. M., et al. (2016). Gene Expression Elucidates Functional Impact of Polygenic Risk for Schizophrenia. Nat. Neurosci. 19 (11), 1442–1453. doi:10.1038/nn.4399
Gaiteri, C., Ding, Y., French, B., Tseng, G. C., and Sibille, E. (2014). Beyond Modules and Hubs: the Potential of Gene Coexpression Networks for Investigating Molecular Mechanisms of Complex Brain Disorders. Genes Brain Behav. 13 (1), 13–24. doi:10.1111/gbb.12106
Garnham, J., Munro, A., Slaney, C., MacDougall, M., Passmore, M., Duffy, A., et al. (2007). Prophylactic Treatment Response in Bipolar Disorder: Results of a Naturalistic Observation Study. J. Affect Disord. 104 (1), 185–190. doi:10.1016/j.jad.2007.03.003
Gene Ontology Consortium, C. (2004). The Gene Ontology (GO) Database and Informatics Resource. Nucleic Acids Res. 32 (Suppl. l_1), 258D–261D. doi:10.1093/nar/gkh036
Horvath, S. (2011). Weighted Network Analysis: Applications in Genomics and Systems Biology. Springer Science & Business Media.
Hou, L., Heilbronner, U., Degenhardt, F., Adli, M., Akiyama, K., Akula, N., et al. (2016). Genetic Variants Associated with Response to Lithium Treatment in Bipolar Disorder: a Genome-wide Association Study. Lancet 387 (10023), 1085–1093. doi:10.1016/S0140-6736(16)00143-4
Hrdlickova, B., de Almeida, R. C., Borek, Z., and Withoff, S. (2014). Genetic Variation in the Non-coding Genome: Involvement of Micro-RNAs and Long Non-coding RNAs in Disease. Biochim. Biophys. Acta 1842 (10), 1910–1922. doi:10.1016/j.bbadis.2014.03.011
Hu, G., Niu, F., Humburg, B. A., Liao, K., Bendi, S., Callen, S., et al. (2018). Molecular Mechanisms of Long Noncoding RNAs and Their Role in Disease Pathogenesis. Oncotarget 9 (26), 18648–18663. doi:10.18632/oncotarget.24307
Hudson, N. J., Reverter, A., and Dalrymple, B. P. (2009). A Differential Wiring Analysis of Expression Data Correctly Identifies the Gene Containing the Causal Mutation. Plos Comput. Biol. 5 (5), e1000382. doi:10.1371/journal.pcbi.1000382
Jafarzadeh, M., and Soltani, B. M. (2020). Long Noncoding RNA LOC400043 (LINC02381) Inhibits Gastric Cancer Progression through Regulating Wnt Signaling Pathway. Front. Oncol. 10, 562253. doi:10.3389/fonc.2020.562253
James, S. L., Abate, D., Abate, K. H., Abay, S. M., Abbafati, C., Abbasi, N., et al. (2018). Global, Regional, and National Incidence, Prevalence, and Years Lived with Disability for 354 Diseases and Injuries for 195 Countries and Territories, 1990-2017: a Systematic Analysis for the Global Burden of Disease Study 2017. Lancet 392 (10159), 1789–1858. doi:10.1016/S0140-6736(18)32279-7
Kidnapillai, S., Bortolasci, C. C., Udawela, M., Panizzutti, B., Spolding, B., Connor, T., et al. (2018). The Use of a Gene Expression Signature and Connectivity Map to Repurpose Drugs for Bipolar Disorder. World J. Biol. Psychiatry 21, 1–28. doi:10.1080/15622975.2018.1492734
Kim, S., Hwang, Y., Webster, M. J., and Lee, D. (2016). Differential Activation of Immune/inflammatory Response-Related Co-expression Modules in the hippocampus across the Major Psychiatric Disorders. Mol. Psychiatry 21 (3), 376–385. doi:10.1038/mp.2015.79
Krebs, C. E., Ori, A. P. S., Vreeker, A., Wu, T., Cantor, R. M., Boks, M. P. M., et al. (2020). Whole Blood Transcriptome Analysis in Bipolar Disorder Reveals strong Lithium Effect. Psychol. Med. 50 (15), 2575–2586. doi:10.1017/S0033291719002745
Kukreja, S., Kalra, G., Shah, N., and Shrivastava, A. (2013). Polypharmacy in Psychiatry: a Review. Mens Sana Monogr. 11 (1), 82–99. doi:10.4103/0973-1229.104497
Langfelder, P., and Horvath, S. (2008). WGCNA: an R Package for Weighted Correlation Network Analysis. BMC Bioinformatics 9 (1), 559. doi:10.1186/1471-2105-9-559
Langfelder, P., Luo, R., Oldham, M. C., and Horvath, S. (2011). Is My Network Module Preserved and Reproducible? Plos Comput. Biol. 7 (1), e1001057. doi:10.1371/journal.pcbi.1001057
Langfelder, P., and Horvath, S. (2022). WGCNA Package FAQ. [Online]. Available at: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html (Accessed March 10, 2022).
Li, S., Jia, H., Zhang, Z., and Wu, D. (2020). LncRNA GAS6-AS1 Facilitates the Progression of Breast Cancer by Targeting the miR-324-3p/SETD1A axis to Activate the PI3K/AKT Pathway. Eur. J. Cell Biol 99 (8), 151124. doi:10.1016/j.ejcb.2020.151124
Liu, Y., Gu, H. Y., Zhu, J., Niu, Y. M., Zhang, C., and Guo, G. L. (2019). Identification of Hub Genes and Key Pathways Associated with Bipolar Disorder Based on Weighted Gene Co-expression Network Analysis. Front. Physiol. 10 (1081), 1081. doi:10.3389/fphys.2019.01081
Liu, Z., Guo, C., Dang, Q., Wang, L., Liu, L., Weng, S., et al. (2022a). Integrative Analysis from Multi-center Studies Identities a Consensus Machine Learning-Derived lncRNA Signature for Stage II/III Colorectal Cancer. eBioMedicine 75, 103750. doi:10.1016/j.ebiom.2021.103750
Liu, Z., Guo, C., Li, J., Xu, H., Lu, T., Wang, L., et al. (2021). Somatic Mutations in Homologous Recombination Pathway Predict Favourable Prognosis after Immunotherapy across Multiple Cancer Types. Clin. Transl Med. 11 (12), e619. doi:10.1002/ctm2.619
Liu, Z., Liu, L., Weng, S., Guo, C., Dang, Q., Xu, H., et al. (2022b). Machine Learning-Based Integration Develops an Immune-Derived lncRNA Signature for Improving Outcomes in Colorectal Cancer. Nat. Commun. 13 (1), 816. doi:10.1038/s41467-022-28421-6
Luykx, J. J., Giuliani, F., Giuliani, G., and Veldink, J. (2019). Coding and Non-coding RNA Abnormalities in Bipolar Disorder. Genes (Basel) 10 (11). doi:10.3390/genes10110946
Mason, O., and Verwoerd, M. (2007). Graph Theory and Networks in Biology. IET Syst. Biol. 1 (2), 89–119. doi:10.1049/iet-syb:20060038
Megiorni, F., Mora, B., Indovina, P., and Mazzilli, M. C. (2005). Expression of Neuronal Markers during NTera2/cloneD1 Differentiation by Cell Aggregation Method. Neurosci. Lett. 373 (2), 105–109. doi:10.1016/j.neulet.2004.09.070
Meng, Z., Gwag, T., Sui, Y., Park, S. H., Zhou, X., and Zhou, C. (2019). The Atypical Antipsychotic Quetiapine Induces Hyperlipidemia by Activating Intestinal PXR Signaling. JCI Insight 4 (3). doi:10.1172/jci.insight.125657
Mentzen, W. I., Floris, M., and de la Fuente, A. (2009). Dissecting the Dynamics of Dysregulation of Cellular Processes in Mouse Mammary Gland Tumor. BMC Genomics 10 (1), 601. doi:10.1186/1471-2164-10-601
Merico, D., Isserlin, R., Stueker, O., Emili, A., and Bader, G. D. (2010). Enrichment Map: A Network-Based Method for Gene-Set Enrichment Visualization and Interpretation. PLOS ONE 5 (11), e13984. doi:10.1371/journal.pone.0013984
Nahman, S., Belmaker, R. H., and Azab, A. N. (2012). Effects of Lithium on Lipopolysaccharide-Induced Inflammation in Rat Primary Glia Cells. Innate Immun. 18 (3), 447–458. doi:10.1177/1753425911421512
Okamura, Y., Aoki, Y., Obayashi, T., Tadaka, S., Ito, S., Narise, T., et al. (2015). COXPRESdb in 2015: Coexpression Database for Animal Species by DNA-Microarray and RNAseq-Based Expression Data with Multiple Quality Assessment Systems. Nucleic Acids Res. 43 (D1), D82–D86. doi:10.1093/nar/gku1163
Ottaviani, S., Stebbing, J., Frampton, A. E., Zagorac, S., Krell, J., de Giorgio, A., et al. (2018). TGF-β Induces miR-100 and miR-125b but Blocks Let-7a through LIN28B Controlling PDAC Progression. Nat. Commun. 9 (1), 1845. doi:10.1038/s41467-018-03962-x
Panizzutti, B., Bortolasci, C. C., Spolding, B., Kidnapillai, S., Connor, T., Richardson, M. F., et al. (2021). Transcriptional Modulation of the Hippo Signaling Pathway by Drugs Used to Treat Bipolar Disorder and Schizophrenia. Int. J. Mol. Sci. 22 (13). doi:10.3390/ijms22137164
Parikshak, N. N., Gandal, M. J., and Geschwind, D. H. (2015). Systems Biology and Gene Networks in Neurodevelopmental and Neurodegenerative Disorders. Nat. Rev. Genet. 16 (8), 441–458. doi:10.1038/nrg3934
Pleasure, S. J., Page, C., and Lee, V. M. (1992). Pure, Postmitotic, Polarized Human Neurons Derived from NTera 2 Cells Provide a System for Expressing Exogenous Proteins in Terminally Differentiated Neurons. J. Neurosci. 12 (5), 1802–1815. doi:10.1523/JNEUROSCI.12-05-01802.1992
Ponting, C. P., Oliver, P. L., and Reik, W. (2009). Evolution and Functions of Long Noncoding RNAs. Cell 136 (4), 629–641. doi:10.1016/j.cell.2009.02.006
Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a Flexible Suite of Utilities for Comparing Genomic Features. Bioinformatics 26 (6), 841–842. doi:10.1093/bioinformatics/btq033
Qureshi, I. A., Mattick, J. S., and Mehler, M. F. (2010). Long Non-coding RNAs in Nervous System Function and Disease. Brain Res. 1338, 20–35. doi:10.1016/j.brainres.2010.03.110
Rege, S., and Hodgkinson, S. J. (2013). Immune Dysregulation and Autoimmunity in Bipolar Disorder: Synthesis of the Evidence and its Clinical Application. Aust. N. Z. J. Psychiatry 47 (12), 1136–1151. doi:10.1177/0004867413499077
R Core Team (2021). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Available at: https://www.R-project.org/
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2009). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26 (1), 139–140. doi:10.1093/bioinformatics/btp616
Savino, A., Provero, P., and Poli, V. (2020). Differential Co-expression Analyses Allow the Identification of Critical Signalling Pathways Altered during Tumour Transformation and Progression. Int. J. Mol. Sci. 21 (24). doi:10.3390/ijms21249461
Seiler, J., Breinig, M., Caudron-Herger, M., Polycarpou-Schwarz, M., Boutros, M., and Diederichs, S. (2017). The lncRNA VELUCT Strongly Regulates Viability of Lung Cancer Cells Despite its Extremely Low Abundance. Nucleic Acids Res. 45 (9), 5458–5469. doi:10.1093/nar/gkx076
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13 (11), 2498–2504. doi:10.1101/gr.1239303
Sonntag, K. C., Woo, T. U., and Krichevsky, A. M. (2012). Converging miRNA Functions in Diverse Brain Disorders: A Case for miR-124 and miR-126. Exp. Neurol. 235 (2), 427–435. doi:10.1016/j.expneurol.2011.11.035
Stahl, S. M. (2013). Stahl's Essential Psychopharmacology: Neuroscientific Basis and Practical Applications. Cambridge University Press.
Susanne, A. B., Fuchs, R., Lackner, N., Birner, A., Reininghaus, B., Meier-Allard, N., et al. (2016). Endoplasmic Reticulum Stress and Bipolar Disorder - Almost Forgotten Therapeutic Drug Targets in the Unfolded Protein Response Pathway Revisited. Cnsnddt 15 (4), 403–413. doi:10.2174/1871527315666160321104613
Statello, L., Guo, C. J., Chen, L. L., and Huarte, M. (2021). Gene Regulation by Long Non-coding RNAs and its Biological Functions. Nat. Rev. Mol. Cell Biol 22 (2), 96–118. doi:10.1038/s41580-020-00315-9
Tabarés-Seisdedos, R., and Rubenstein, J. L. R. (2009). Chromosome 8p as a Potential Hub for Developmental Neuropsychiatric Disorders: Implications for Schizophrenia, Autism and Cancer. Mol. Psychiatry 14 (6), 563–589. doi:10.1038/mp.2009.2
Ulitsky, I. (2016). Evolution to the rescue: Using Comparative Genomics to Understand Long Non-coding RNAs. Nat. Rev. Genet. 17 (10), 601–614. doi:10.1038/nrg.2016.85
van Dam, S., Võsa, U., van der Graaf, A., Franke, L., and de Magalhães, J. P. (2018). Gene Co-expression Analysis for Functional Classification and Gene-Disease Predictions. Brief Bioinform 19 (4), 575–592. doi:10.1093/bib/bbw139
Weighill, D., Ben Guebila, M., Glass, K., Platig, J., Yeh, J. J., and Quackenbush, J. (2021). Gene Targeting in Disease Networks. Front. Genet. 12. doi:10.3389/fgene.2021.649942
Wu, M., Yang, L.-Z., and Chen, L.-L. (2021). Long Noncoding RNA and Protein Abundance in lncRNPs. Rna 27 (12), 1427–1440. doi:10.1261/rna.078971.121
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS 16 (5), 284–287. doi:10.1089/omi.2011.0118
Zhang, P., Dong, Q., Zhu, H., Li, S., Shi, L., and Chen, X. (2019). Long Non-coding Antisense RNA GAS6-AS1 Supports Gastric Cancer Progression via Increasing GAS6 Expression. Gene 696, 1–9. doi:10.1016/j.gene.2018.12.079
Zhang, Z.-Q., Wu, W.-W., Chen, J.-D., Zhang, G.-Y., Lin, J.-Y., Wu, Y.-K., et al. (2021). Weighted Gene Coexpression Network Analysis Reveals Essential Genes and Pathways in Bipolar Disorder. Front. Psychiatry 12 (320). doi:10.3389/fpsyt.2021.553305
Zhao, J., Sun, B. K., Erwin, J. A., Song, J. J., and Lee, J. T. (2008). Polycomb Proteins Targeted by a Short Repeat RNA to the Mouse X Chromosome. Science 322 (5902), 750–756. doi:10.1126/science.1163045
Zhou, K. R., Liu, S., Sun, W. J., Zheng, L. L., Zhou, H., Yang, J. H., et al. (2017). ChIPBase v2.0: Decoding Transcriptional Regulatory Networks of Non-coding RNAs and Protein-Coding Genes from ChIP-Seq Data. Nucleic Acids Res. 45 (D1), D43–D50. doi:10.1093/nar/gkw965
Keywords: bipolar disorders, co-expression network, WGCNA, mood stabilizers, lncRNAs, treatments, mood disorders
Citation: Truong TT, Bortolasci CC, Spolding B, Panizzutti B, Liu ZS, Kidnapillai S, Richardson M, Gray L, Smith CM, Dean OM, Kim JH, Berk M and Walder K (2022) Co-Expression Networks Unveiled Long Non-Coding RNAs as Molecular Targets of Drugs Used to Treat Bipolar Disorder. Front. Pharmacol. 13:873271. doi: 10.3389/fphar.2022.873271
Received: 10 February 2022; Accepted: 24 March 2022;
Published: 08 April 2022.
Edited by:
Francisco Lopez-Munoz, Camilo José Cela University, SpainReviewed by:
Lifeng Li, First Affiliated Hospital of Zhengzhou University, ChinaLu Jia, Capital Medical University, China
Copyright © 2022 Truong, Bortolasci, Spolding, Panizzutti, Liu, Kidnapillai, Richardson, Gray, Smith, Dean, Kim, Berk and Walder. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Trang TT. Truong, truongtra@deakin.edu.au