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⁻⁵ 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 × 10⁵ 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).

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).

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 × 10⁶ 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).

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).

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.

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.

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.

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).

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.

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.

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.

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.

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.

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.

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).

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.