Astrocyte fatty acid metabolism as a driver of risk for major depressive disorder

Domain: neuroscience
Published: 2026-04-08
Editorial outcome: major_revision

Abstract
Genome-wide association studies (GWAS) have successfully identified genetic loci associated with major depressive disorder (MDD), yet the complex gene networks underpinning this polygenic risk remain largely uncharacterized. Here we show the neurobiological mechanisms of MDD by analyzing co-expression networks of 94 risk genes in the human prefrontal cortex. By linking these networks to individual symptoms, we identify the FADS1 (fatty acid desaturase 1) network as a central integrator across symptom clusters. We find that the FADS1 network functions primarily in astrocytes to regulate fatty acid metabolism and influence oligodendrocyte-related cell states. Furthermore, we identify FGF2 as a synaptic effector of this pathway and highlight PPARα (peroxisome proliferator-activated receptor alpha) as a putative therapeutic target. These results establish astrocyte fatty acid metabolism as a critical mechanistic contributor to MDD and a promising avenue for treatment.

Introduction
Major depressive disorder (MDD) is a leading contributor to the global burden of disease and is associated with immense socio-economic consequences1,2. Pharmacotherapy is typically the primary treatment strategy for MDD but large clinical trials show high relapse rates with common anti-depressant drugs3. Identifying novel therapies is a priority but is hampered by our currently incomplete mechanistic understanding of MDD.

Genome-wide association studies (GWAS) show conceptual promise to identify MDD mechanisms in an unbiased fashion. Indeed, analysis on GWAS summary statistics have nominated many hundreds of genes that contribute to MDD risk. However, the genetic architecture of MDD is highly polygenic, and as such risk is not driven by individual genes, but through complex gene networks 4,5. Genetically identified risk genes are likely important components of these networks, but in isolation have a minimal contribution to overall risk 5,6.

Co-expression networks are one method to identify the broad networks through which genes act to influence phenotypic variation. They are typically built in an unbiased fashion using RNA-sequencing datasets 7, but an alternative method is to construct seeded co-expression networks, where networks are built around candidate genes such as those identified from post-GWAS analyses. Seeded co-expression networks have been used successfully to interrogate the mechanistic underpinnings of psychiatric disorders 8,9 and previous work has shown this approach can accurately recapitulate CRISPR loss of function experiments 10. In this study, w e hypothesized that generating and characterizing seeded co-expression networks for MDD risk genes would allow us to identify novel mechanisms of MDD for potential therapeutic innovation (Figure 1).

We first identified MDD risk genes by annotating MDD GWAS summary statistics to the gene level using chromatin interaction data from the human dorsolateral prefrontal cortex (DLPFC). The DLPFC is among the most well-studied brain regions in MDD, with evidence from post-mortem, model systems, and imaging studies as well as clinical trials pointing towards an important role of the DLPFC in MDD 11–14. We created seeded co-expression networks for these genes using RNA-seq data from the DLPFC of the GTEx dataset. We validated these networks in internal and external datasets and then associated the most robust networks with the nine individual symptoms of MDD assessed by the patient health questionnaire 9 (PHQ9).

These networks clustered into distinct patterns of association with the individual MDD symptoms and gene ontology (GO) enrichment analysis showing that the clustered networks represented distinct functional mechanisms. The FADS1 (fatty acid desaturase 1) network (i.e. genes positively correlated with FADS1 expression in the DLPFC) was the network most densely connected with other network clusters and emerged as an important candidate for further study. The FADS1 network was associated with fatty acid signalling and primarily expressed in astrocytes. Integration of the FADS1 network with snRNA-seq data from control and MDD samples showed that the FADS1 network could define an astrocyte population depleted in MDD. This depleted astrocyte population was associated with the blood-brain barrier (BBB) and the oligodendrocyte lineage. Further analyses suggested the FADS1 network drove transcriptomic changes in excitatory neurons through changes in astrocyte-excitatory neuron communication, putatively driven by FGF2 (fibroblast growth factor 2). Finally, we identified PPARα (peroxisome proliferator activated receptor alpha) as a promising therapeutic candidate for further investigation in MDD.

Results
Generating seeded co-expression networks to characterize the mechanisms of MDD risk genes
We first identified MDD risk genes by annotating an MDD GWAS 15 to the gene level using Hi-C data from the DLPFC16,17 (Supplementary Figure 1). We identified 153 genes with significant associations following correction for multiple comparisons (hereafter referred to as MDD risk genes, Figure 2A; Supplementary Data 1). The MDD risk genes from HMAGMA were not enriched in any DLPFC co-expression networks identified from any of the eight previous studies we assessed 6–13, all of which used a WGCNA (weighted correlation network analysis) approach (Supplementary Figure 2). Although WGCNA, as shown by many studies, is a powerful approach for mechanistic inference 6–15, considering the broad distribution of MDD risk genes across modules we reasoned a seeded co-expression approach may be more informative for mechanistic inference in this instance. From the 153 MDD risk genes we generated seeded co-expression networks using RNA-seq data from the DLPFC of the GTEx catalogue (“Frontal_cortex_BA9”; 94 MDD risk genes remained after filtering lowly expressed genes).

We tested each MDD risk gene for co-expression with each of the 23,326 genes in the RNA-seq dataset (following filtering, normalisation and regression of covariates; see methods section). Any gene with a correlation coefficient > 0.5 or < -0.5 and an FDR P-value < 0.05 was considered to be positively or negatively co-expressed, respectively. Together all genes positively or negatively co-expressed with a particular seed were considered to form positive or negative co-expression networks. Therefore, we use the word “network” to refer to the set of correlations with an MDD risk gene, and as such does not capture correlations between genes that are not the MDD risk gene (as is done in WGCNA networks for instance).

A positive co-expression network can be considered as the genes with direct functional associations with the MDD risk gene, while the negative co-expression network can be conceptualised as the processes downstream of those positively correlated genes. There was a strong correlation (r = 0.89; Pearson correlation) between the number of genes positively and negatively correlated with a given MDD risk gene (Figure 2B; Supplementary Data 2). Remarkably, 55.4% of all genes tested for co-expression were co-expressed with at least one MDD risk gene (Figure 2C), illustrating a striking degree of complexity.

Processes related to neuronal activity were among the most common gene ontology (GO) terms enriched across all networks, with similar terms enriched in both positive and negative networks, suggesting oppositional regulation of processes by MDD risk gene networks (Figure 2D and 2E).

Of all genes co-expressed genes with an MDD risk gene, 78.9% were co-expressed with more than one seed (Supplementary Figure 3A; Supplementary Data 3). There was also a dose-dependent inverse relation between the number of networks a particular gene was present within and their loss-of-function constraint (LOUEF score) in the gnomAD database 28 (Figure 2F; Kruskal-Wallis p < 2.2e-16). We observed a similar pattern for missense mutation Z-score (Supplementary Figure 3B; Kruskal-Wallis p < 2.2e-16), but a more modest association with synonymous mutation Z-scores that was also statistically significant (Supplementary Figure 3C; Kruskal-Wallis p = 0.027). Psychiatric risk genes are typically highly intolerant to mutations, prompting us to ask if they were present in more networks than the background set of brain-expressed genes. Indeed, MDD risk genes (as defined in this study), SFARI high-confidence autism genes, and schizophrenia risk genes (those associated with rare variants 29) were present in more networks than expected (Figure 2G; Kruskal-Wallis p < 2.2e-16).

We next performed further enrichment analyses using genes divided into three categories: i) those not co-expressed with an MDD risk gene, ii) those present in <=5 networks, and iii) those present in >5 networks. Genes present in >5 networks and thus highly implicated in MDD-related pathways were enriched in genes differentially expressed in rodent models of chronic stress 30 (chronic variable stress-nucleus accumbens p=1.3e-6; chronic variable stress- prefrontal cortex p=1.8e-4; social isolation-nucleus accumbens p=6.6e-17; social isolation- prefrontal cortex p=3.3e-5; Fischer’s Exact Test; Figure 2H). Genes present in >5 networks were also enriched in GWAS of common variants for bipolar disorder 31 (p=2.2e-7; MAGMA) and schizophrenia 32 (p=1.5e-5; MAGMA; Figure 2H). Genes present in <= 5 networks were enriched in the schizophrenia GWAS (p=9.2e-5; MAGMA), but not for other psychiatric disorders 33–37. We found no enrichment using genes not co-expressed with MDD risk genes, suggesting seeded co-expression networks may be an efficient way to investigate pleiotropic mechanisms across psychiatric disorders.

Integration of networks with PHQ9 GWAS identifies symptom-specific MDD molecular associations
We validated the reproducibility of our co-expression networks using internal (adjacent tissue from the GTEx consortium with largely overlapping individuals; Supplementary Figure 4A and Supplementary Figure 5C) and external samples (using data from phase 1 of the CommonMind consortium 38,39; Supplementary Figure 4B and Supplementary Figure 5B). In both validations, we assessed the correlation between each of our primary co-expression networks and the replication dataset. For both validation cohorts, we observed variable replications (internal r: range -0.06 to 0.9, mean = 0.72; external r: range -0.25 to 0.8, mean = 0.43). The replication for individual networks was highly correlated with the number of genes co-expressed with the MDD risk gene (Supplementary Figure 4C; Supplementary Data 5; GTEx internal rho= 0.85 and p < 2.2e-16; CommonMind external rho=0.55, p=1.7e-8; Spearman correlation), which plateaued at 1000 co-expressed genes.

We also note substantial overlap of the genes within these co-expression networks, with 90% and 82% of the internal and external datasets, respectively enriched for genes in their corresponding co-expression networks from the original dataset (Supplementary Figure 5). These data demonstrate networks with a larger number of genes are more reproducible. We therefore constrained our subsequent analyses to those networks with more than 1000 co-expressed genes.

To examine mechanisms proximal to each MDD risk gene, we used GWAS enrichment analysis with our positive co-expression networks in responses to the PHQ-9 included as part of the UK Biobank. We note that our initial characterization of MDD risk genes using the MDD GWAS from Als et al 15, was performed using summary statistics that excluded the UK Biobank sample. This approach avoided the reuse of the UK Biobank in our PHQ-9 symptom-specific analyses. As guided by our previous validations, we used positive co-expression networks with more than 1000 co-expressed genes in this analysis (16 networks).

Hierarchical clustering of the standardized betas from the GWAS enrichment analyses suggested three groups of networks (Figure 3A). Cluster 1 (11 networks) was enriched for GO terms primarily related to neuronal function (Figure 3B). Cluster 2 (2 networks) was enriched for GO terms related to metabolic processes such as lipid metabolism (Figure 2C), and enrichment analyses (Figure 3A) suggested a preferential association with suicidal ideation and feelings of inadequacy. Cluster 3 (3 networks) was enriched for GO terms related to inflammation, myelination and the vasculature (Figure 3D).

We next investigated the inter-network interactions using network analysis, which allowed us to identify densely connected networks (Figure 3E). We generated several network-based metrics including degree, betweenness, closeness, eigenvector centrality and pagerank (Figure 3F ; Supplementary Figure 6). These metrics largely showed similar patterns across networks, with the notable exception of the FADS1 positive network in betweenness, which had by far the highest betweenness of any network (Figure 3F), with dense connections to several networks from the other clusters. We also observed a similar pattern in the internal GTEx validation dataset, with the FADS1 positive network having the largest betweenness of any network. However, in the external CMC validation dataset, the FADS1 positive network did not have the highest betweenness (Supplementary Figure 7). Nevertheless, we considered the FADS1 positive network as an attractive candidate for further investigation, as FADS1 has also been associated with MDD in several other analyses including Mendelian randomization 3, transcriptome-wide association studies 4 and differential expression in human post-mortem tissue 5. The FADS1 positive network was also primarily associated with suicidal ideation, making it the most tractable to investigate through integration with post-mortem MDD datasets. Furthermore, the FADS1 positive network was among the most well-replicated networks in our internal (r = 0.89; enrichment P-value <2.22-16) and CMC external validations (r = 0.61; enrichment P-value <2.22-16). Furthermore, we were able to validate the FADS1 positive network in two other external datasets using DLPFC tissue from control samples (Labonté et al 46 r = 0.59; enrichment P-value <2.22-16 and Chang et al 38, r = 0.52; enrichment P-value <2.22-16).

The FADS1 positive network is astrocyte-specific and altered in MDD
We next aimed to further characterize the functional, cell-specific and MDD-associated characteristics of the FADS1 positive network. Predictably, GO analysis of the network revealed enrichment for processes related to the metabolism of fatty acids (Figure 4A). Encouragingly, the FADS1 positive network was enriched for protein-protein interactions in STRINGdb 42 (Supplementary Figure 8; p < 1.0e-16), suggesting our RNA co-expression is also related to protein function. A down sampling and bootstrapping approach to estimate the number of genes in the FADS1 positive network across the 11 distinct GTEx regions showed a large number of co-expressed genes across all regions, with the largest network in the DLPFC (Supplementary Figure 9A). After the DLPFC, the anterior cingulate cortex, hippocampus and nucleus accumbens held the next largest FADS1 positive networks. Each of these regions is strongly implicated in the neurobiology of MDD 43. The FADS1 positive network was enriched in several WGCNA modules from the human DLPFC across numerous studies (Supplementary Figure 9B), and we found a signature highly overlapping to the FADS1 positive network when we used sparse partial least squares to derive a latent variable of DLPFC gene expression correlated with FADS1 expression (909 of 949 genes overlapped with the FADS1 positive network; Supplementary Figure 9C). Together these data suggest the seeded coexpression approach is reproducible with other distinct methods. The FADS1 positive network was also enriched in mouse (chronic social defeat stress 44) and rat (selectively bred line 45) models of depression-like behaviours (Figure 4B), suggesting it is responsive to the environment.

We then used snRNA-seq data from the DLPFC of control samples 13 to determine the cell specificity of the FADS1 positive network, which showed striking specificity for astrocytes (Figure 4C). We also, however, saw enrichment of the FADS1 positive network in astrocyte genes differentially expressed between non-human primates (gorilla and chimpanzee) and humans (Figure 4D). Together these data suggest the FADS1 positive network is associated with depression-relevant behaviours across species, but there are also likely important human-specific features of the network.

We next asked whether the expression of FADS1 positive network genes was different in human post-mortem MDD samples compared to controls. We thus used snRNA-seq data from the DLPFC of control and MDD samples 13 to investigate astrocyte-specific expression of the FADS1 positive network. We generated a module score and aggregated it to the sample level before comparing across groups, which showed a significant decrease in expression of the FADS1 positive network in MDD (Figure 4E; p = 0.018; Wilcoxon test). We next used a bootstrapping approach with 10,000 iterations to assess the composition of the FADS1 positive network in bulk RNA-seq data from the DLPFC of control and MDD samples 46. We found a striking decrease in the number of co-expressed genes in MDD compared to a sample size-matched control group (Figure 4F; p <2.22-16; Wilcoxon test). We could, reassuringly, reproduce this finding using a microarray dataset 47 (Figure 4G; p <2.22-16; Wilcoxon test). FADS1 expression did not change between control and MDD samples, suggesting a complex relation between FADS1 expression and the number of genes co-expressed with FADS1 (Supplementary figure 9D and 9E).

We then used strict criteria (present in >95% of bootstrapped networks), to define co-expression networks across control and MDD groups, with a high degree of confidence. This approach allowed us to compare FADS1 positive network composition across control and MDD samples. The MDD group had 54 genes co-expressed with FADS1 in more than 95% of bootstrapped networks (hereafter high-confidence co-expressed genes), with only 3 genes common with the control group (Figure 4H). Meanwhile, the control group had 215 high-confidence co-expressed genes, again with only 3 being common with the MDD group (Figure 4I; Supplementary Data 6). The MDD-specific co-expressed genes were enriched for GO terms related to lipoprotein function (Figure 4J; Supplementary Data 7) while, in contrast, the control-specific genes were enriched for GO terms related to myelination and the extracellular matrix (Figure 4K). Astrocyte fatty acid metabolism is essential for correct myelination 48, and decreased myelination in the PFC of patients with MDD has previously been observed 49. Together these data demonstrate the number of genes co-expressed with FADS1 decreases in MDD, which may result in a decreased capacity for myelination.

The FADS1 positive network can define astrocyte states altered in MDD
Astrocytes are a heterogenous cell type with intricate communication patterns with all brain cell types. We aimed to further dissect the impact of the FADS1 positive network on astrocyte function and cell interactions in MDD using a DLPFC snRNA-seq dataset from control and MDD individuals 13. We extracted astrocytes from the dataset using previously defined broad cell type classes 13 and clustered the astrocytes using only expression of the FADS1 positive network, rather than the most variable genes as done in a typical analysis pipeline. This approach allowed us to ask whether the FADS1 positive network could define astrocyte clusters and then ask if these clusters were altered in MDD. Indeed, we could identify three clusters of astrocytes using the FADS1 positive network (Figure 5A-D; Supplementary Figure 10A and 10B; Supplementary Data 8). All three clusters expressed canonical astrocyte genes (Figure 5E). Cluster 0 was enriched for GABRB1 (Gamma-Aminobutyric Acid Type A Receptor Subunit Beta 1) and the extracellular matrix regulator HPSE2 (heparanase 2; Figure 5E). Interestingly the GABA-A receptor, a sub-unit of which is encoded by GABRB1, is the target for the recently approved anti-depressant brexnanolone 50. Cluster 1 was primarily enriched for genes related to blood-brain barrier (BBB) function and GO terms related to cell adhesion, functions which have previously been associated with MDD 51,52 (Figure 5E; Supplementary Figure 9C). Cluster 2 had a notable enrichment for genes and GO terms related to the synapse (Figure 5E; Supplementary Figure 10C), several of which have also previously been associated with MDD 15.

We then assessed the spatial localization of these astrocyte clusters using the human DLPFC spatial transcriptomic dataset from Maynard et al 53. Genes enriched (positive logFC and FDR p-value < 0.05) within each of the three astrocyte clusters (Figure 5A), as well as genes enriched in broad classes of other canonical cell types, were used to build module scores in the spatial transcriptomic data (representative scores of the astrocyte clusters can be seen in Figure 5F). We then correlated the spatial localisations of these scores across 12 distinct sections and calculated the mean spatial correlation between cell type scores (Figure 5G). In line with our previous analysis, there was a strong colocalization between the Cluster 2 score and a score from excitatory neurons, suggesting a role for this subset at excitatory synapses. Cluster 1 was strongly expressed in layer I and the white matter and was highly colocalized with scores from the oligodendrocyte lineage and endothelial cells, once more suggesting a role for this subtype at the BBB. The Cluster 0 score mirrored the correlations of Cluster 1 but with weaker associations suggesting it was more functionally diffuse (Figure 5G).

We next used Milo 54 to identify differentially abundant neighbourhoods using our FADS1 positive defined dimensional reductions. We observed depletion of several neighbourhoods in MDD, which were predominantly localised to Cluster 1. Differential expression analysis between the neighbourhoods depleted in MDD and non-depleted neighbourhood showed a strong signature of BBB genes in the MDD-depleted neighbourhoods (Figure 5I; Supplementary Figure 10D; Supplementary Data 9). We used GSEA to test the enrichment of the FADS1 positive network in these differential-expression statistics and found it highly enriched (Figure 5J). In this analysis there was a polarisation of the genes in the FADS1 positive network to both ends of the differential expression statistics, but particularly those genes enriched in MDD-depleted neighbourhoods. GO term enrichment for the genes differentially expressed between depleted and non-depleted neighbourhoods showed enrichment for synapse-related terms in the non-depleted neighbourhoods (Figure 5K), with enrichment for GO terms related to BBB function within the markers of depleted neighbourhoods (Figure 5L).

We replicated key aspects of these findings in an independent MDD snRNA-seq dataset from Maitra et al 14. For instance, we again identified three clusters of astrocytes using the FADS1 positive network, with enrichment for similar GO terms (Supplementary Figure 11). We also identified MDD-depleted astrocytes using Milo and these depleted neighbourhoods were also enriched for BBB-related genes (Supplementary Figure 12). Correlation of the differential expression statistics for the astrocyte-depleted and non-depleted signatures between the two datasets showed a strong correlation (Supplementary Figure 13; R = 0.6 and p < 2.2e-16; Pearson correlation), suggesting we are capturing a reproducible pattern of astrocyte depletion that is defined by the FADS1 positive network.

Taken together, these data demonstrate the FADS1 positive network is involved in diverse astrocyte functions within distinct astrocyte clusters. Further, there is a reproducible depletion of an astrocyte population associated with the BBB and the oligodendrocyte lineage in MDD.

FGF2 is a putative regulator of these altered astrocyte states on downstream synaptic function
We next asked what mechanisms might act downstream of the FADS1 positive network. We reasoned that a parsimonious extension of our analysis, which has focused on genes positively correlated with FADS1, would be to investigate genes negatively correlated with FADS1 as the genes likely downstream of the FADS1 network. The FADS1 negative network numbered 388 genes, which were primarily expressed in excitatory neurons (Figure 6A), enriched in GO terms related to synaptic function (Figure 6B) and enriched at an uncorrected p-value threshold in several psychiatric-related GWAS (including MDD and suicide attempt; Figure 6C). We did not observe preferential expression of the FADS1 negative network in excitatory neuron subtypes, suggesting the FADS1 negative network represents a general pattern of expression common to DLPFC excitatory neurons (Supplementary Figure 14). This is perhaps not surprising considering the potential role of fatty acids in cellular metabolism.

Astrocyte and neuronal functions are intimately linked, therefore, our characterization of the FADS1 positive network led us to hypothesize the FADS1 positive and FADS1 negative networks may be linked through cell-cell communication. To test this, we leveraged multiple published datasets to define the subcellular enrichment of the FADS1 positive and FADS1 negative networks. We first used data from a mouse study that characterized the sub-cellular localisation of astrocyte proteins using an in vivo proximity labelling approach 55. The FADS1 positive network was enriched for several subcellular astrocyte compartments, with the strongest enrichment seen in proteins proximate to GLT1, which is localised to the astrocyte portion of the tripartite synapse and is involved in glutamate uptake (Figure 6D).

We then used synapse-annotated genes from the synGO database 56 and proteins present at the human post-synaptic density 57 to determine if the FADS1 negative was also enriched at the synapse. Indeed, we saw a marked enrichment of the FADS1 negative network in both synaptic gene sets (Figure 6E). Together these data suggest the FADS1 negative network acts within astrocytes and affects the expression of synapse-related genes in excitatory neurons.

We then used NicheNet 58 to identify putative astrocyte ligands (that were in the FADS1 positive network) with the regulatory potential to affect excitatory neuron genes in the FADS1 negative network. This analysis nominated FGF2 as a putative candidate to drive astrocyte-excitatory neuronal signalling in the context of FADS1 co-expressed genes (Figure 6F). In rodents, astrocyte FGF2 is differentially modulated by both stress and anti-depressants, and FGF2 has been associated with astrocyte-neuron communication 59–61. This finding also agrees with a large body of literature using both human and animal models to describe a role for FGF signalling in MDD 62–64.

We then wondered if genes within the global FGF signalling pathway were differentially expressed in astrocytes of MDD samples. We used the GO term “FGF signalling pathway” to define genes involved in FGF signalling and a pseudobulk analysis with the voom-limma method to test for differential expression within our previously defined astrocyte clusters (from Figure 5A). We observed a striking downregulation of FGF pathway genes in astrocyte Cluster 1 (which was characterized by a BBB signature and neighbourhoods depleted in MDD; Supplementary figure 14D; Supplementary Data 10). Out of 79 genes tested for differential expression in Cluster 1, 40 passed correction for multiple comparisons and 75 had a negative fold change. These data show profound differential expression of FGF signalling components within astrocytes associated with the BBB and the oligodendrocyte lineage in MDD and suggest a role for FGF2 in driving downstream changes in synapse-related gene expression in excitatory neurons.

PPARα as a putative therapeutic modulator in MDD
We finally sought to identify druggable therapeutic targets from our analyses (Figure 7A). Our analyses identified two prominent and reproducible signatures that were disrupted in MDD, which represented promising signatures to therapeutically target; the 212 genes that were only co-expressed in control samples (and not MDD; Figure 4I and 4K), and the genes that were enriched in MDD-depleted astrocyte neighbourhoods (Figure 5I). The DrugMatrix database is among the largest databases assaying the transcriptomic effects of various drugs and therefore, represented an ideal database in which to test for enrichment of our gene signatures 65,66. This analysis revealed 622 and 481 compounds associated with co-expressed genes lost in MDD and the markers of MDD-depleted astrocytes, respectively (Figure 7B and Supplementary Data 11 and 12). We note these analyses define agents that modulate our signatures and may therefore be pro- or anti-depressant. However, we observed enrichment for several drugs used to treat MDD (e.g. paroxetine, fluoxetine and citalopram), other psychiatric disorders (e.g. quetiapine, fluphenazine, olanzapine, clozapine, lorazepam and diazepam) and drugs that target previously postulated MDD-related mechanisms, such as inflammation (e.g. mefenamic acid), glucocorticoid signalling (e.g. hydrocortisone and dexamethasone) and sex-hormone signalling (e.g. diethylstilbestrol). We also saw enrichment for multiple agonists of PPARα (including bezafibrate, clofibrate, fenofibrate and gemfibrozil), several of which appeared among the most enriched drugs. Together the enrichment we see for well-validated and characterized MDD-related drugs lends credence to more novel targets, such as PPARα.

The DrugMatrix database contains signatures across multiple cell lines and concentrations. We therefore reasoned that an association across multiple conditions would provide stronger evidence for an association. Among the most frequently associated drugs were two that target PPARα (Fenofibrate and Bezafibrate). Moreover, each of the two drugs were associated with both gene signatures (Figure 7C). PPARα is a nuclear receptor that regulates the expression of a range of lipid-modulating genes 67 and so represented an intriguing candidate. Encouragingly PPARA was also a member of the FADS1 positive network (Supplementary Figure 15A , r = 0.69; Pearson correlation), but other PPAR family members PPARD and PPARG were not co-expressed with FADS1 (Supplementary Figure 14A and 14C; r = 0.02 and -0.12; Pearson correlation).

We finally used the Open Targets platform to search for genetic associations with PPARA. Surprisingly there was an MDD GWAS significant locus 88Kb from PPARA (Figure 7D; Supplementary Data 13) 68, making it the second closest gene to the locus (WNT7B was the closest at 84Kb). We could not identify any brain expression QTLs for PPARA, despite its robust expression, which negated the possibility of confirming the association using colocalization or Mendelian randomization. There was however also evidence for PPARA as the causal gene in multiple other GWAS, including several pertaining to co-morbidities or risk factors for MDD including addiction, inflammation and body weight (Figure 7D). Together these data support the further study of PPARα as a potential therapeutic target in MDD.

Methods
The research published within this manuscript was conducted on previously collected data and therefore not subject to ethical approval related to data collection. All relevant ethical guidelines were stringently adhered to through the analysis process and access to protected data in the case of GTEx were through approval (accession number phs000424.v9.p2).

Annotation of MDD GWAS using Hi-C data with H-MAGMA
For the MDD GWAS annotation we used the largest publicly available summary statistics (Als et al), excluding the UK Biobank and the 23 and me sample. We excluded the UK Biobank sample to eliminate the potential for bias introduced by the reuse of datasets (we use GWAS summary statistics from the UK Biobank PHQ-9 in Figure 3 to associate co-expression networks with specific MDD symptoms). The 23 and me sample was not made publicly available and as such could not be included 15. The MDD summary statistics were downloaded from iPSYCH as described in the original publication 15 (https://ipsych.dk/en/research/downloads/). We then annotated these summary statistics to the gene level using Hi-C data from the DLPFC 16, which was obtained from https://zenodo.org/record/6382668#.Y9PSG3bMJPY, as previously described through the H-MAGMA framework 17,100. We used SNP annotations from the 1000 genomes European dataset and gene annotations from the NCBI website build 38. We chose H-MAGMA over MAGMA because the region-specific approach to annotating the GWAS summary statistics aligned nicely with our region-specific approach to seeded co-expression networks

Bulk RNA sequencing quality control and normalisation
GTEx v8 data were used because of their well validated RNA extraction and sequencing methodology, large sample size and well documented technical phenotype data. Raw gene-level read counts, transcript-level expected counts, sample attributes, and sequencing metrics were downloaded from the GTEx portal and dbGaP (accession number phs000424.v9.p2). Only samples and genes passing GTEx QC were retained. Analyses were restricted to Brain – Frontal Cortex (BA9). Samples were excluded if marked ineligible or if the donor had the following: amyotrophic lateral sclerosis, Alzheimer’s or other dementia, encephalitis, influenza, multiple sclerosis, Parkinson’s disease, stroke or cerebrovascular disorder, systemic lupus, bacterial or viral sepsis, active cancer, schizophrenia, major depressive disorder, substance use (heroin or prescription drug abuse), pneumonia, unexplained seizures, or positive blood cultures. After filtering, 111 samples remained. Genes with >0.1 counts per million (CPM) in >25% of samples were retained (23,326 genes). Counts were normalized using the TMM method in edgeR v3.38.1, log₂-transformed, and filtered for low expression. Outlier samples were identified by connectivity Z-scores from biweight midcorrelations (|Z| > 3) and removed (n = 4). Thirty-two sequencing metrics (e.g., mapping rate, read quality, duplication rate, insert size) were standardized and summarized by PCA; the first ten components were used as covariates. Additional covariates included RNA integrity number (RIN), Hardy scale, refrigeration interval, post-mortem interval category, sex, and age. Expression residuals were obtained by ordinary least squares regression of normalized counts on all covariates. The resulting residualized expression matrix was used for all downstream co-expression analyses.

Labonté et al: Data from Labonté et al 46 were downloaded from GSE102556 in FPKM form. Genes were filtered to those with an FPKM (Fragments Per Kilobase of transcript per Million mapped reads) > 0.5 in 25% of samples, leaving 22893 genes in the dataset. The effects of age, age2, pH, sex, post-mortem interval, RIN (RNA integrity number) and RIN2 were then regressed from the data.

Chang et al: Microarray data were used as in a previous meta-analysis 38. Data were normalised using the “rma” function in the affy v1.72.0 package 102. Outlying samples were identified by a connectivity Z-score greater than three standard deviations and were removed. The following covariates were then regressed from the data using the “ComBat” function in the sva v3.42.0 package 103; age, race, sex, post-mortem interval, pH, RIN, RNA degradation (quantified with the “AffyRNAdeg” function) and batch. Probe IDs were converted to gene names using the biomaRt package in R.

Generation of seeded co-expression networks
Co-expression networks were generated using a Pearson correlation between the MDD risk gene and all genes expressed in the dataset. As per previous studies 9, a gene was said to be positively co-expressed with an R > 0.5 and an FDR p-value < 0.05, a gene was said to be negatively co-expressed with an R < -0.5 and an FDR < 0.05. The FDR was calculated for each network independently and based on the 23,326 correlations conducted. Bootstrapping was carried out using 10,000 iterations while down sampling to the smallest sample size among compared groups.

Mutational intolerance analysis
Mutation metrics were downloaded from the gnomAD database v2.1.1 28. We used the LOEUF (loss-of-function observed/expected upper bound fraction) score as a conservative measure of mutational intolerance, along with missense Z-score and synonymous Z-score as provided. All genes expressed in the GTEx dataset that were tested for co-expression were used as background in the analysis.

Enrichment analyses
Fisher’s Exact Test was conducted to test enrichment between various gene sets of interest as implemented in the GeneOverlap v1.30.0 package 104. The number of genes from which the test gene set was generated, were used as the background in each analysis.

Gene set enrichment analysis (GSEA)
GSEA was conducted using the R package fgsea v1.20.0 105 using differential gene expression statistics ranked by log fold change. Where necessary rodent gene lists were converted to their human orthologues using the R biomaRt package v2.50.3 106. The normalisation function of the fgsea package has been reported as unstable with gene sets larger than 500, therefore, we used the 500 genes with the strongest correlation co-efficient with FADS1 in the analysis. Of note, a similar analysis with the full complement of FADS1 correlated genes produced similar results.

Internal and external validation
Internal validation was conducted using the GTEx “Brain_cortex” data, which is considered a technical replicate of the “Brain_frontal_cortex_BA9” dataset that forms the basis of the majority of our analyses. The “Brain_cortex” data underwent dissection and processing using a different protocol and at a different site to the “Brain_frontal_cortex_BA9” dataset, but both datasets contain overlapping individuals. These data were downloaded and analysed in an identical manner using the same covariates as the “Brain_frontal_cortex_BA9” analysis.

External validation was done using the phase 1 data from the DLPFC of the common mind consortium as described previously 38. Data were filtered and normalized in an identical way to the GTEx data, leaving 26,599 genes expressed at a CPM > 0.1 in more than 25% of samples. These genes were normalized using the TMM method and log normalisation. Outliers with a connectivity score more than 3 standard deviations were removed. The effect of covariates was regressed from the data (covariates: Age, Age2, Batch, pH, Sex, PMI_hrs, RIN, RIN2 and the first 10 sequencing principal components). Co-expression networks were then generated in an identical way to the GTEx data.

GWAS enrichment
We conducted GWAS enrichment using MAGMA v1.10 100, with GWAS annotated to the gene level using H-MAGMA 17. We used SNP annotations from the 1000 genomes European dataset and gene annotations from the NCBI website build 38. We annotated GWAS using Hi-C data from the DLPFC 16, which was obtained from https://zenodo.org/record/6382668#.Y9PSGS3bMJPY.

For the PHQ9 analysis, the beta coefficients from the regression (produced by the gene-set analysis function in MAGMA) were scaled for each network and hierarchically clustered using the pheatmap v1.0.12 package 107. The dendrograms were then cut to form three clusters of networks.

For details of GWAS summary statistics used and their sources see Supplementary Data 14.

Gene ontology (GO) analysis
We conducted GO analysis using the clusterProfiler v4.2.2 package 108, with all genes expressed in the relevant dataset used as background in the analysis. We ran enrichment for biological processes, cellular components, and molecular functions in a single analysis. For Figure 3B-3D, significant terms were imported into Cytoscape v3.8.0 109, where similar terms were clustered and a word cloud was used to suggest words common to all terms within a given cluster. We then manually refined these terms for brevity and clarity.

Network analysis
Network analysis on co-expression networks was conducted with igraph v1.5.1 110. To generate the adjacency matrix, we first identified the co-expression networks with more than 1000 positively co-expressed genes (defined as the number of genes with an FDR<0.05 and R>0.5). Then for these networks we correlated the coefficients of all genes, with those from the other networks. Specifically, the correlation matrix was generated from an input matrix containing 23,326 rows (the number of genes against which co-expression was tested) and 14 columns (the number of networks with more than 1000 genes positively co-expressed). The adjacency matrix was then constructed from this correlation matrix by thresholding the presence or absence of a correlation. Edge thickness was set to represent the weight between nodes. The layout was automatically determined using the Fruchterman-Reingold layout algorithm. The betweenness was then calculated on this adjacency network using the “betweenness” function within igraph.

Protein-protein interactions
Estimates of protein-protein interactions were made using STRINGdb 42 using standard settings. Protein-protein interaction networks were then exported to Cytoscape v3.8.0 109 for plotting.

Sparse partial least squares (sPLS)
A sparse partial least squares (sPLS) analysis was performed using the mixOmics package (v6.18.1) to derive a latent variable capturing the expression pattern of FADS1 in the GTEx Brain Frontal Cortex (BA9) dataset. Expression data were preprocessed as previously described. The predictor matrix (X) included all genes except FADS1, and the response vector (y) contained normalized FADS1 expression values. An sPLS model with one component and 1,183 selected genes (keepX = 1183) was fit to mirror the number of genes coexpressed with FADS1 (FDR < 0.05, R > 0.5) identified in our seeded coexpression network. The latent variable was sign-adjusted to ensure a positive correlation with FADS1 expression. The latent was highly correlated with FADS1 expression r = 0.86 and p-value < 2.2e-16.

Single cell expression of networks
We downloaded the Nagy et al raw data from GSE144136 and used the cell type assignations from the original study, but grouped them to form broad cell classes. We excluded nuclei expressing fewer than 700 and greater than 7000 genes, before normalization with the “NormalizeData” function. Data were then split and integrated by batch using Seurat v4.3.0 111 and the standard Seurat pipeline before using the “AddModuleScore” function to identify cell type specific expression patterns.

To compare module scores between control and MDD, module scores were aggregated by sample and condition before comparison with a Wilcoxon test.

Analysis of snRNA-seq data
We used the Nagy et al and Maitra et al data as they are the largest published snRNA-seq datasets for MDD. To cluster astrocytes based on variation in expression of genes within the FADS1 positive network we started with the Nagy et al raw data. We subset this data to nuclei annotated as astrocytes and removed two batches as they had too few nuclei to perform integration (batches 1M and 2M). We then subset the count matrix to genes within the FADS1 positive network and ran the standard Seurat integration and analysis pipeline as previously described. We used the first 10 principal components and a resolution of 0.2 to generate clusters and a UMAP dimensional reduction of the nuclei. We then extracted the metadata with the FADS1 positive based clusters and combined it with the full astrocyte count matrix to perform downstream analysis, such as identification of cluster marker genes, differential expression and differential abundance.

Maitra et al dataset: We downloaded the Maitra et al 14 data from GSE213982 and analysed the data in an identical pipeline to the Nagy et al dataset. We used the first 10 principal components and a resolution of 0.1 to generate clusters and a UMAP dimensional reduction of the nuclei.

Differential abundance and gene expression testing with miloR
We conducted differential abundance testing between control and MDD cells using the miloR v1.2.0 package 54. We started with our FADS1 positive clustered object containing the full count matrix and converted it to a Single Cell Experiment object and subsequently to a Milo object. We used the following parameters for the analysis; k = 20, d = 20 and prop = 0.2. We then tested neighbourhoods for differential abundance using batch and age as covariates. Differentially abundant neighbourhoods were defined as those with a spatial FDR < 0.1. The differentially abundant neighbourhoods with a negative log fold change (i.e. decreased in MDD), were tested for marker genes by comparison with all other neighbourhoods using log normalized counts with the “findNhoodGroupMarkers” function in miloR.

Spatial analysis of module scores
We downloaded 10X Visium data from Maynard et al (samples: Br2743, Br2720, Br3942, Br6432, Br6471, Br6522, Br8325, Br8492, Br8667) 53. We normalized data using the Seurat pipeline and the “SCTransform” function. We created a module score using the “AddModuleScore” function. This score was created using all genes enriched in broad cell types or the FADS1 positive defined astrocyte clusters. Marker genes were defined as those with a positive logFC and an FDR P-value < 0.05. Correlations between cell types were generated independently for all slices. We then averaged the correlation coefficients of all slices to obtain the final correlation matrix for plotting.

Differential gene expression in the FGF signalling pathway
To test for differential gene expression within the FGF pathway we used genes annotated to the GO term “fibroblast growth factor receptor signalling pathway” (GO:0008543). The astrocyte count matrix were subset to these genes and differential analysis run using pseudobulk counts with the limma-voom method implemented within the Libra R package, which has been shown to be among the best performers in single cell differential expression testing 112.

Drug analysis
We conducted drug enrichment analysis using the DrugMatrix databases and the iLINCS software 66. We used the online interface provided and standard settings to identify patterns of association (Supplementary Data 11 and 12). Results were returned in the form of absolute Z-score. We aggregated the significant results by perturbagen to identify targets associated with multiple concentrations or cell lines.

Open Targets
We used the Open Targets Genetics platform 113 to search for GWAS hits associated with PPARA. We conducted our search in December 2022. We manually annotated the results into broad groups to facilitate interpretation. We provide the data using the Open Targets L2G (locus-to-gene) score, which is analogous to a probability score based on a machine learning pipeline trained to identify likely causal genes from significant GWAS hits. The L2G score can be interpreted as probability the gene is the target of the causal locus under the assumption the GWAS is similar to those included in the training set.

Statistical analysis
All analysis unless otherwise stated were carried out using R v4.1.1 114 and Rstudio v1.4.1717 115. GWAS enrichment was conducted using the MAGMA package. All package versions are specified in their relevant methods section. Plots were generated using either the packages mentioned throughout the methods section or in ggplot2 v3.4.0 116. Plots were combined into figures using Inkscape v1.0. The Benjamini-Hochberg method (False Discovery Rate; FDR) was used to correct for multiple comparisons where appropriate, except for annotation of the MDD GWAS where a Bonferroni correction was applied. The relevant statistical tests used are described throughout the results section and the figure legends.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "this is a huge amount of data presented combining human work from different sources, crosslinking this to animal work and\neventually also testing therapeutic options. the authors researched every corner possible and need to be commended on the\ntremendous work. the manuscript is extremely dense yet clear and logical. Figures are supplementary to the work and while\nthere are many, illustrate the findings and the text nicely.\nhowever, when it comes to the methods it was really hard to follow what the authors had done and why. i got lost in all the\nabbreviations used and would like to see them written out. also, some background or at least some more info of why studies\nfor comparison were chosen and why methods were used. i recognise a good number of empty brackets where i believe\nreferences are missing for the methods and tests described. some additional sentences would enlighten the reader as to the\nwhy and what.\ni wasn't able to install and run the code", "reviewer_2": "Fitzgerald and colleagues conducted an integrative network analysis using major depressive disorder (MDD) risk genes as\nseed genes to explore the neurobiological mechanisms underlying the disorder. They identified the FASD1 network that is\nassociated with fatty acid metabolism as a central connector across network clusters. Additionally, the network is mainly\nexpressed in astrocytes. They further identified a putative regulator of the network and eventually showed a potential\ntherapeutic target for MDD. While several interesting findings were presented, I have some concerns that need to be\naddressed.\n1. The identification of the gene coexpression network positively correlated with FADS1 as a key connector across network\nclusters, based on the highest betweenness, is a significant finding in this study. Please clarify why utilizing the highest\nbetweenness as a criterion for identifying key connectors is appropriate for exploring the mechanisms underlying MDD?\n2. The expression level of FADS1 did not show significant changes between control and MDD groups in this study. If\nFADS1 does not exhibit differential expression, it is important for the authors to explain how it can still play a pivotal role\nwithin the network clusters and in the mechanisms underlying the development of MDD.\n3. The authors attempt to replicate the initially generated networks using both internal and external samples and evaluate\ncorrelations between the seeded networks and replication datasets. However, more thorough replication is necessary to\nvalidate the results of the network analysis. I recommend that the authors conduct seed network analysis using the\nreplication datasets and demonstrate the statistical significance of the overlap in module genes between the networks from\nthe initial study and those generated from the replication datasets. Additionally, it should be clarified whether a consensus\nnetwork positively correlated with FADS1 from the replication datasets also exhibits the highest betweenness centrality.\n4. The authors hypothesize that a network positively correlated with the seed genes is likely to represent direct functional\nassociations with those genes, while a network negatively correlated with the genes may relate to downstream processes.\nPlease elaborate on the rationale behind this hypothesis and provide evidence or theoretical frameworks that support this\nidea.", "reviewer_3": "In this article, the authors explore the functional relationship of MDD risk genes and their potential mechanistic impact in the\netiopathogenesis of the disorder by applying a seeded co-expression network approach to RNA-seq data of healthy\ncontrols.\nI am positively impressed by the amount of genomic and transcriptomic analyses implemented to obtain converging\nevidence from multiple data sources with the aim of disentangling the complex architecture of MDD. The article reads well.\nThroughout the manuscript, the authors offer a clear explanation of the hypotheses behind the analyses implemented to test\nthem. When interpreting the results, they propose an intriguing biological mechanism for MDD implicating fatty acid\nmetabolism in astrocytes and a potential therapeutic target.\nI think that some sections of the manuscript could benefit from additional detail and clarification. Additionally, I have strong\ndoubts about the technical foundation of the network analysis. I offer the following suggestions:\n• The foundation of the article is based on the fact that the authors apply H-MAGMA to identify MDD risk genes and since\nthey do not find enrichment for these genes in any of the two WGNA-based co-expression networks considered (Hartl et\nal.,2021; Punzi et al. 2022), they propose the seeded co-expression approach as alternative. Although the authors briefly\nmention in the limitation section of the discussion the problem relative to annotating GWAS SNPs to genes, I believe a\nfurther investigation of the other methods could be useful and informative. For example, why consider the H-MAGMA genes\nas MDD risk genes instead of the ones prioritized by the GWAS from Als. et al. (the one used for the H-MAGMA analysis)? Is\nthere an overlap between the genes identified? What is the rationale behind preferring H-MAGMA to other methods?\nConsidering that there are multiple mapping strategies, such as physical sequence proximity (MAGMA), gene expression\nregulation (eQTLs), or long-range interactions by measuring chromatin interaction profiles (H-MAGMA), it could be more\ninformative to either derive a consensus of the different methodologies or, even better, the union of genes obtained by the\ndifferent methods.\n• On the same line, I reckon more robust evidence is needed to explain why the “seeded\nco-expression approach may be more informative for mechanistic inference”. Besides the only two ones analyzed, there are several WGCNA-based co-\nexpression networks on the DLPFC that should be investigated to obtain a broader picture of how the MDD risk genes are\nfound or not found co-expressed in this brain region (Fromer et al.,2016; Gandal et al., 2018; Pergola et al., 2017,2019,2024;\nRadulescu et al., 2020). Additionally, networks available nowadays include a relevant age dimension that is not considered\nby Hartl and Punzi, who analyzed adult transcriptomics. Consider Li et al., 2019, Walker et al., 2019, Werling et al., 2020,\nPergola et al., 2023. There should be a greater effort to reconcile previously published co-expression evidence with this new\nwork. The seeded approach has historically been somewhat intriguing, though very old and, in many ways, much more\nlimited than WGCNA. The validity of seeding a single gene or a few genes rests on the quality of the assay, which is\nobviously lower for a single measure than for transcriptome-wide measures. And if the same evidence cannot be retrieved\nby any published network out there, this casts doubt on the possibility that this unique approach is free of artifacts. A recent\nstudy has recapitulated CRISPR evidence—one criterion that authors cite in support of their approach—very effectively\nusing published WGCNA networks (Borcuk et al., 2024). To investigate the context of specific risk genes, the question of\nwhy this should be more informative than analyzing the WGCNA modules where these genes are found co-expressed still\nlingers and I honestly cannot agree with such a statement. I point to this as a crucial comparative analysis, because the\ninterpretation that “all genes positively or negatively co-expressed with a particular seed were considered to form positive or\nnegative co-expression networks” and that “ a positive co-expression network can be considered as the genes with direct\nfunctional associations with the MDD risk gene, while the negative co-expression network can be conceptualised as the\nprocesses downstream of those positively correlated genes” can be misleading and should be taken cautiously, since the\nseeded approach does not give any information of the relationship between the genes identified co-expressed with a seed\ngene (which is what a network implies). This is even disregarding the fact that unsigned networks incorporate both positive\nand negative relationships, so this is not an exclusive attribute of the seeded approach. As pointed out by the authors,\nalmost 80% of the co-expressed genes are co-expressed with more than one seed. This means that the seeded approach\ngives you a close-up of a specific seed gene but does not capture the molecular environment surrounding those genes. It\ncould very well be that one gene is negatively correlated with multiple seed genes as a result of a transcription factor\nregulating the seed genes or genes positively correlated with the seed genes that belong to a specific biological pathway.\nUnfortunately, this common and plausible biological scenario remains undetectable by the seeded approach implemented,\nwhile it is well captured by WGCNA or other co-expression network analyses. Finally, if each new study picks their favorite\ngenes, the applicability of the findings is questionable, especially if the same findings prove impossible to retrieve using\nunbiased, data-driven approaches.\nMinor suggestions:\n• In the generation of seeded co-expression networks the number of multiple comparisons used in the FDR correction is not\nmentioned. Is it 94 or 94 by 23,326?\n• In the validation analyses for the co-expression networks, it is not clearly defined what are the terms of the “correlation\nbetween each of our primary co-expression networks and the replication dataset”. I am interpreting that for each seed gene\nyou correlated the correlation coefficients of the co-expressed genes found in the primary and replication analysis. Is that the\ncase? If you performed the exact same analysis in the replication datasets (R > |0.5| and an FDR p-value < 0.05), what is the\noverlap of the positively and negatively correlated genes between the primary and replication datasets? Do you find the\nexact same genes?\n• In the Network analysis on the co-expression networks conducted with igraph, it is not clear how the weighted adjacency\nmatrix for all networks was constructed. What is represented in each cell of the adjacency matrix?"}, "author_response": "Response to Reviewer’s comments\nAstrocyte fatty acid metabolism as a driver of risk for major depressive disorder\nEamon Fitzgerald, Nicholas O'Toole, Irina Pokhvisneva, Eric J Nestler, Gustavo Turecki,\nCorina Nagy, Michael J Meaney\nWe would first like to thank all the Reviewers for their positive and constructive\ncomments on our manuscript.\nWe believe the strength of our manuscript lies in its contextualization of genetic risk for\nMDD through a network-based approach. We generate and interrogate MDD risk\nnetworks using multiple complementary and orthogonal computational methods,\nintegrating diverse datasets and data modalities. Notably, these distinct approaches\nconverge on astrocyte fatty acid metabolism as a novel mechanism underlying MDD.\nThe additional analyses recommended by the reviewers, including further technical\nvalidations, have further reinforced our conclusions and enhanced the overall\nrobustness of our study.\nPlease find a detailed response to each of the comments below. We believe our\nmanuscript has been greatly improved through this process and we very much\nappreciate the Reviewers’ time and effort.\nThroughout this response:\n The Reviewer’s original comments appear in black\n Our response appears in blue\n Text quoted from the revised manuscript appears in red\nReviewer #1\nThis is a huge amount of data presented combining human work from different sources,\ncrosslinking this to animal work and eventually also testing therapeutic options. The\nauthors researched every corner possible and need to be commended on the\ntremendous work. the manuscript is extremely dense yet clear and logical. Figures are\nsupplementary to the work and while there are many, illustrate the findings and the text\nnicely.\nhowever, when it comes to the methods it was really hard to follow what the authors had\ndone and why. i got lost in all the abbreviations used and would like to see them written\nout. also, some background or at least some more info of why studies for comparison\nwere chosen and why methods were used. I recognise a good number of empty\nbrackets where i believe references are missing for the methods and tests described.\nsome additional sentences would enlighten the reader as to the why and what.\nWe thank the Reviewer for their kind words. We thoroughly revised the methods section\nto increase clarity, describe motivation for the choice of the various datasets and\nremove excess abbreviations. Regarding the empty parentheses noticed by the\nReviewer, we had initially used these parentheses as a signal to the reader that we\nwere referring to a coding function. For instance, when we referenced the function to\nnormalize single cell RNA-seq data using the R package Seurat, we referred to the\n“NormalizeData() function”. To run such a function in R, you would then populate the ()\nwith the name of the Seurat object for normalization e.g.\n“NormalizeData(SeuratObject)”. We now realize this caused confusion and have\nswitched instead to using “” when referring to a function’s name. Please see below for a\nmore explicit example.\nWe also now clearly describe our rationale for including datasets in our analysis. In brief\nour criteria typically amounted to datasets that were publicly available, had a sufficiently\nlarge sample size and bore relevance to our study (e.g. tissue from the dorsolateral\nprefrontal cortex).\nPlease see below for several examples of changes made in the revised manuscript to\naddress these points.\nLine 478-480 of the methods: Example explication of acronyms\n“Data were normalized using DESeq2 v1.34.0 with the TMM (trimmed mean of M\nvalues) method and log normalization. Genes were filtered to those with >0.1 CPM\n(counts per million) in >25% of samples”\nLine 599 of the methods: Example removal of empty brackets and denoting coding\nfunctions with “”\n“We created a module score using the “AddModuleScore” function”\nLine 454-459 of the methods: Example 1 of increased clarity and explanation for dataset\nchoice\n“For the MDD GWAS annotation we used the largest publicly available summary\nstatistics (Als et al), excluding the UK Biobank and the 23 and me sample. We excluded\nthe UK Biobank sample to eliminate the potential for bias introduced by the reuse of\ndatasets (we use GWAS summary statistics from the UK Biobank PHQ-9 in Figure 3 to\nassociate co-expression networks with specific MDD symptoms). The 23 and me\nsample was not made publicly available and as such could not be included.”\nLine 466-468 of the methods: Example 2 of increased clarity and explanation of dataset\nchoice\n“GTEx v8 data were used because of their well validated RNA extraction and\nsequencing methodology, large sample size and well documented technical phenotype\ndata.”\nReviewer #1 (Remarks on code availability):\ni wasn't able to install and run the code\nWe thank the Reviewer for noting this point. We have double-checked the code and\ndata, which should be available for download (this has been confirmed by several of our\ncolleagues). We also now provide a detailed description of the files available in the\nrepository and the parts of the manuscript to which they pertain. We provide the link\nhere for ease of access: https://zenodo.org/records/10181579.\nReviewer #2\nFitzgerald and colleagues conducted an integrative network analysis using major\ndepressive disorder (MDD) risk genes as seed genes to explore the neurobiological\nmechanisms underlying the disorder. They identified the FASD1 network that is\nassociated with fatty acid metabolism as a central connector across network clusters.\nAdditionally, the network is mainly expressed in astrocytes. They further identified a\nputative regulator of the network and eventually showed a potential therapeutic target\nfor MDD. While several interesting findings were presented, I have some concerns that\nneed to be addressed.\n1. The identification of the gene coexpression network positively correlated with FADS1\nas a key connector across network clusters, based on the highest betweenness, is a\nsignificant finding in this study. Please clarify why utilizing the highest betweenness as a\ncriterion for identifying key connectors is appropriate for exploring the mechanisms\nunderlying MDD?\nThe Reviewer raises an important point. Indeed, there are many measures for\nquantifying networks and our initial presentation of these metrics now seems too\nnarrow. In the revised manuscript we provided several additional network-based\nmeasures (degree, closeness, eigenvector centrality and pagerank). We provided these\ndata below in Reviewer Figure 1 and in Supplementary Figure 5 of the revised\nmanuscript. We also revised the manuscript to better explain our choice for highlighting\nbetweenness and the FADS1 network. We provide this text below.\npositive\nReviewer Figure 1 Additional network-based metrics pertaining to the graph in Figure\n3E of the revised manuscript. The colours within each plot refer to the cluster identities\nderived from the clustering of PHQ9 symptoms in Figure 3.\nLine 176-179 of the results section: Description of network-based measures and reason\nfor focusing on betweenness and the FADS1 network\npositive\n“These metrics largely showed similar patterns across networks, with the notable\nexception of the FADS1 network in betweenness, which had by far the highest\npositive\nbetweenness of any network (Figure 3F), with dense connections to several networks\nfrom the other clusters.”\n2. The expression level of FADS1 did not show significant changes between control and\nMDD groups in this study. If FADS1 does not exhibit differential expression, it is\nimportant for the authors to explain how it can still play a pivotal role within the network\nclusters and in the mechanisms underlying the development of MDD.\nWe expanded our discussion in our revised manuscript to address this issue. The\nprotein encoded by FADS1 is a fatty acid desaturase, which catalyses the conversion of\ndihomo-γ-linolenic acid and eicosatetraenoic acid into arachidonic acid (AA) and\neicosapentaenoic acid (EPA) 1,2. Both AA and EPA are active molecules with a variety of\nascribed functions, ranging from metabolism to inflammation 1–3. Our results within the\ncontext of the relevant literature suggest AA and EPA are likely central players in the\nassociations we describe with the FADS1 network.\npositive\nAlthough FADS1 expression itself does not change significantly between control and\nMDD groups, altered function of the FADS1 encoded protein may arrive from differential\nsubstrate availability and consequentially differential production of AA and EPA. We\nsuggest two possible routes for altered substrate availability: altered astrocyte fatty acid\nmetabolism or dietary availability. We now expound upon this point in the discussion of\nour revised manuscript.\nThe primary consideration here is that our emphasis is on a FADS1 gene network rather\nthan on a specific single gene. We see downregulation of the FADS1 network\npositive\n(rather than the FADS1 gene itself) in MDD and altered abundance of FADS1\npositive\nnetwork defined astrocyte states in MDD. These findings align with the systems biology\nperspective of considering the network itself (rather than an individual gene or protein)\nas the unit in which disease risk manifests 4,5. Therefore, changes within the dynamics\nof the network could alter the context of FADS1 encoded protein function and thereby\nalter its relation with MDD phenotypes.\nWe believe the Reviewer’s point also touches on a key limitation of our study, notably\nwe are using RNA-seq data to infer functional pathway activity occurring in a different\nmodality (i.e. lipids). We have revised the discussion section to further emphasize this\nlimitation of our study.\nLine 398-405 of the discussion: Discussion of potential FADS1 expression-independent\ncontributions to MDD.\n“We did not observe differential expression of FADS1 itself but instead observed\ndownregulation of the whole FADS1 gene network, which aligns with the systems\npositive\nbiology perspective of considering the network itself (rather than an individual gene or\nprotein) as the unit in which disease risk manifests 4,5. Together with the previously\npublished literature, this data implies that astrocyte fatty acid substrate availability, may\nmediate the effect of FADS1 on the FADS1 network and downstream MDD\npositive\nphenotypes. The capacity for diet-supplied fatty acids to modulate these pathways and\ndisease states is an outstanding question and one that is of central importance.”\nLine 437-440 of the discussion: Acknowledgment our study is limited to inferring lipid\npathway activity via RNA-seq data.\n“Second, we also relied on transcriptomic methods throughout this study. However,\nthese methods are limited in their ability to infer events that may not be well represented\nat the transcriptome level, such as lipid metabolism. Emerging single-cell metabolomic\nmethods promise to build upon our results.”\n3. The authors attempt to replicate the initially generated networks using both internal\nand external samples and evaluate correlations between the seeded networks and\nreplication datasets. However, more thorough replication is necessary to validate the\nresults of the network analysis. I recommend that the authors conduct seed network\nanalysis using the replication datasets and demonstrate the statistical significance of the\noverlap in module genes between the networks from the initial study and those\ngenerated from the replication datasets. Additionally, it should be clarified whether a\nconsensus network positively correlated with FADS1 from the replication datasets also\nexhibits the highest betweenness centrality.\nWe agree with the Reviewer, further validations would strengthen our study. In line with\nthe Reviewer’s comments, we conducted several analyses, which we outline below and\nwhich are now included in the revised manuscript.\n1. We tested the enrichment of our co-expression networks in their corresponding\nnetworks among the validation datasets. We provided these data below in\nReviewer Figure 2 and in Supplementary Figure 4 of the revised manuscript. We\nobserved a similar pattern as in our original validation, i.e. some networks\nreplicate better than others and the strength of the replication is dependent of\nthe number of co-expressed genes.\n2. We generated graphs and graph-based metrics for both the internal and external\ndatasets as suggested by the Reviewer. In our internal validation dataset, we\nalso found that the FADS1 network has, by a large margin, the highest\npositive\nbetweenness- replicating our original data. However, we did not see this in our\nexternal validation dataset. We believe reasons for this may include biological\n(e.g. age, sex) or technical (e.g. precise dissection sites, dissection protocols,\nRNA extraction and sequencing protocols) differences between datasets. We\ntempered our language around the betweenness finding in the results section of\nthe revised manuscript, noting it only validated in the internal dataset. However,\nwe believe the FADS1 network warranted further investigation despite the\npositive\nlimited validation data. First, it was among the most well-replicated networks in\nour internal (r = 0.89; enrichment FDR P-value < 2.2 × 10⁻¹⁶) and external\nvalidations (r = 0.61, enrichment FDR P-value < 2.2 × 10⁻¹⁶). Second, many\nprevious studies have implicated FADS1 in MDD 6–8. We present these new data\nbelow in Reviewer Figure 3 and in Supplementary Figure 6 of the revised\nmanuscript.\nReviewer Figure 2 Enrichment analyses for the original networks in those derived from\nthe validation datasets. In A) the percentage of networks significantly enriched (or not\nenriched) for genes in their corresponding networks in both the internal and external\nvalidation datasets is shown. In B) and C) the correlation between the odds ratio\n(derived from a Fisher’s exact test) and the number of co-expressed genes in a network\nis presented for the external (B) and internal (C) validations.\nLines 155-156 of the results: Discussion of enrichment validations\n“We also note substantial overlap of the genes within these co-expression networks,\nwith 90% and 82% of the internal and external datasets, respectively, enriched for\ngenes in their corresponding co-expression networks from the original dataset\n(Supplementary Figure 4).”\nReviewer Figure 3 Betweenness metric for graphs derived from the internal validation\n(A) and external (B) datasets. The color scheme is consistent across both plots.\nLine 179-188 of the results section: Discussion of betweenness validations\n“We also observed a similar pattern in the internal GTEx validation dataset, with the\nFADS1 network having the largest betweenness of any network. However, in the\npositive\nexternal CMC validation dataset, the FADS1 network did not have the highest\npositive\nbetweenness (Supplementary Figure 6). Nevertheless, we considered the FADS1\npositive\nnetwork as an attractive candidate for further investigation, as FADS1 has also been\nassociated with MDD in several other analyses including Mendelian randomization 6,\ntranscriptome-wide association studies 7 and differential expression in human post-\nmortem tissue 8. The FADS1 network was also primarily associated with suicidal\npositive\nideation, making it the most tractable to investigate through integration with post-\nmortem MDD datasets. Furthermore, the FADS1 network was among the most\npositive\nwell-replicated networks in our internal (r = 0.89; enrichment FDR P-value < 2.2 × 10⁻¹⁶)\nand external validations (r = 0.61, enrichment FDR P-value < 2.2 × 10⁻¹⁶).”\n4. The authors hypothesize that a network positively correlated with the seed genes is\nlikely to represent direct functional associations with those genes, while a network\nnegatively correlated with the genes may relate to downstream processes. Please\nelaborate on the rationale behind this hypothesis and provide evidence or theoretical\nframeworks that support this idea.\nWe thank the Reviewer for flagging these sentences. As noted by the Reviewer, in the\noriginal manuscript we suggested negatively correlated genes represent downstream\nprocesses. This statement is potentially misleading. For instance, take a simple\nexample where a transcription factor activates the expression of one gene and\nrepresses the expression of a different gene. The transcription factor would be positively\ncorrelated with the first gene and negatively correlated with the second gene. Therefore,\nto suggest the downregulated gene is downstream of the transcription factor, but the\nupregulated gene is not, would be unjustified.\nWe therefore revised our manuscript to discuss the nature of these correlations with\nmore nuance, acknowledging the points noted above.\nLine 102 -108 of the results: Explanation of seeded co-expression networks\n“Any gene with a correlation coefficient > 0.5 or < -0.5 and an FDR P-value < 0.05 was\nconsidered to be positively or negatively co-expressed, respectively. Together all genes\npositively or negatively co-expressed with a particular seed were considered to form\npositive or negative co-expression networks. Therefore, we use the word “network” to\nrefer to the set of correlations with an MDD risk gene, and as such does not capture\ncorrelations between genes that are not the MDD risk gene (as is done in WGCNA\nnetworks for instance).”\nReviewer #3 (Remarks to the Author):\nIn this article, the authors explore the functional relationship of MDD risk genes and\ntheir potential mechanistic impact in the etiopathogenesis of the disorder by applying a\nseeded co-expression network approach to RNA-seq data of healthy controls.\nI am positively impressed by the amount of genomic and transcriptomic analyses\nimplemented to obtain converging evidence from multiple data sources with the aim of\ndisentangling the complex architecture of MDD. The article reads well.\nThroughout the manuscript, the authors offer a clear explanation of the hypotheses behind the analyses\nimplemented to test them. When interpreting the results, they propose an intriguing\nbiological mechanism for MDD implicating fatty acid metabolism in astrocytes and a\npotential therapeutic target.\nI think that some sections of the manuscript could benefit from additional detail and\nclarification. Additionally, I have strong doubts about the technical foundation of the\nnetwork analysis. I offer the following suggestions:\n• The foundation of the article is based on the fact that the authors apply H-MAGMA to\nidentify MDD risk genes and since they do not find enrichment for these genes in any of\nthe two WGNA-based co-expression networks considered (Hartl et al.,2021; Punzi et al.\n2022), they propose the seeded co-expression approach as alternative. Although the\nauthors briefly mention in the limitation section of the discussion the problem relative to\nannotating GWAS SNPs to genes, I believe a further investigation of the other methods\ncould be useful and informative. For example, why consider the H-MAGMA genes as\nMDD risk genes instead of the ones prioritized by the GWAS from Als. et al. (the one\nused for the H-MAGMA analysis)? Is there an overlap between the genes identified?\nWhat is the rationale behind preferring H-MAGMA to other methods? Considering that\nthere are multiple mapping strategies, such as physical sequence proximity (MAGMA),\ngene expression regulation (eQTLs), or long-range interactions by measuring chromatin\ninteraction profiles (H-MAGMA), it could be more informative to either derive a\nconsensus of the different methodologies or, even better, the union of genes obtained\nby the different methods.\nThe Reviewer raises a valid point that merits comment in the manuscript (please see\nbelow). There are alternative ways to annotate GWAS summary statistics to the gene\nlevel. We did not use the gene level annotations provided in the original Als et al\nmanuscript because we used a subset of their data, which excluded the UK Biobank\n(and the 23 and me sample; which was not made publicly available). We used this\nsubset to eliminate concerns about reuse of datasets and the potential for introducing\nbias therein, stemming from our use of the UK Biobank to run enrichment for MDD\nsymptoms (derived from the UK Biobank PHQ9). For a similar reason we did not use a\nTWAS approach (as many large eQTL studies include GTEx, where we generate our\nco-expression networks).\nWe then chose H-MAGMA over MAGMA because the region-specific approach to\nannotating the GWAS summary statistics aligned nicely with our region-specific\napproach to seeded co-expression networks.\nLine 440-443 of the discussion: Acknowledgment of the limitations of GWAS annotation.\n“Third, there are many methods to annotate GWAS to the gene level with no accepted\nbest approach. A more holistic understanding of the genetic architecture of disease will\nallow for more accurate identification of causal genes from GWAS loci and facilitate their\ninterpretation through studies such as ours.”\n• On the same line, I reckon more robust evidence is needed to explain why the “seeded\nco-expression approach may be more informative for mechanistic inference”. Besides the only two ones analyzed, there are several WGCNA-based co-expression networks\non the DLPFC that should be investigated to obtain a broader picture of how the MDD\nrisk genes are found or not found co-expressed in this brain region (Fromer et al.,2016;\nGandal et al., 2018; Pergola et al., 2017,2019,2024; Radulescu et al., 2020).\nAdditionally, networks available nowadays include a relevant age dimension that is not\nconsidered by Hartl and Punzi, who analyzed adult transcriptomics. Consider Li et al.,\n2019, Walker et al., 2019, Werling et al., 2020, Pergola et al., 2023. There should be a\ngreater effort to reconcile previously published co-expression evidence with this new\nwork. The seeded approach has historically been somewhat intriguing, though very old\nand, in many ways, much more limited than WGCNA. The validity of seeding a single\ngene or a few genes rests on the quality of the assay, which is obviously lower for a\nsingle measure than for transcriptome-wide measures. And if the same evidence cannot\nbe retrieved by any published network out there, this casts doubt on the possibility that\nthis unique approach is free of artifacts. A recent study has recapitulated CRISPR\nevidence—one criterion that authors cite in support of their approach—very effectively\nusing published WGCNA networks (Borcuk et al., 2024). To investigate the context of\nspecific risk genes, the question of why this should be more informative than analyzing\nthe WGCNA modules where these genes are found co-expressed still lingers and I\nhonestly cannot agree with such a statement. I point to this as a crucial comparative\nanalysis, because the interpretation that “all genes positively or negatively co-expressed\nwith a particular seed were considered to form positive or negative co-expression\nnetworks” and that “ a positive co-expression network can be considered as the genes\nwith direct functional associations with the MDD risk gene, while the negative co-\nexpression network can be conceptualised as the processes downstream of those\npositively correlated genes” can be misleading and should be taken cautiously, since\nthe seeded approach does not give any information of the relationship between the\ngenes identified co-expressed with a seed gene (which is what a network implies). This\nis even disregarding the fact that unsigned networks incorporate both positive and\nnegative relationships, so this is not an exclusive attribute of the seeded approach. As\npointed out by the authors, almost 80% of the co-expressed genes are co-expressed\nwith more than one seed. This means that the seeded approach gives you a close-up of\na specific seed gene but does not capture the molecular environment surrounding those\ngenes. It could very well be that one gene is negatively correlated with multiple seed\ngenes as a result of a transcription factor regulating the seed genes or genes positively\ncorrelated with the seed genes that belong to a specific biological pathway.\nUnfortunately, this common and plausible biological scenario remains undetectable by\nthe seeded approach implemented, while it is well captured by WGCNA or other co-\nexpression network analyses. Finally, if each new study picks their favorite genes, the\napplicability of the findings is questionable, especially if the same findings prove\nimpossible to retrieve using unbiased, data-driven approaches.\nThe Reviewer makes several points each of which we respond to below. We will make\nthe general point that we are big proponents of WGCNA (and indeed have used the\napproach in many of our publications). We are certainly not suggesting seeded co-\nexpression is a better approach than WGCNA, only that in our manuscript using a\nseeded approach allowed us to provide mechanistic insight into MDD. We added more\nnuance to our discussion of WGCNA throughout the manuscript and its ability to provide\nmechanistic inference, as exemplified by the Borcuk et al paper.\n As suggested by the Reviewer we increased the diversity of WGCNA modules in\nwhich we test enrichment for MDD risk genes. We added WGCNA modules from\nsix additional studies, referenced by the Reviewer. We provide these data below\nin Reviewer Figure 4 and in Supplementary Figure 1 of the revised manuscript.\nWe again found no enrichment among any of these published WGCNA modules\nfor the MDD risk genes.\n Motivated by the Reviewer’s surprise at a lack of association between our data\nand published WGCNA modules, we conducted a further analysis in which we\ntested enrichment of the FADS1 network within these WGCNA modules.\npositive\nThis analysis tested the hypothesis that it is our seeded co-expression networks,\nrather than the full swathe of MDD risk genes, that hold a relation to the WGCNA\nmodules. Indeed, we found a very strong enrichment of the FADS1 network\npositive\namong WGCNA modules generated in every dataset. Interestingly our\nFADS1 network was enriched for multiple WGCNA modules in each dataset,\npositive\nsuggesting it is capturing a related but distinct signature to these modules.\nTherefore, these analyses show there is a relation between our seeded co-\nexpression networks and published WGCNA modules. We present this analysis\nbelow in Reviewer Figure 5.\n Without question the methodology described in Borcuk et al is very powerful and\nindeed based on this impression we have several preliminary projects applying\nthese methods in the lab. We agree with the Reviewer seeded co-expression\nnetworks are an “old” method, but we believe they still allow for mechanistic\ninsight and ultimately it is mechanistic insight for which we are searching. We\nadded several sentences to the discussion expanding on these points, and in\nparticular noting how methods such as Borcuk et al will likely facilitate further\nmechanistic insight in the future.\n The word “network” does not have a strict definition, and upon reflection we\nrealize that it would be easy for the reader to misinterpret our use of “network”.\nWe thank the Reviewer for this point. We revised our manuscript to\noperationalize the word “network” and thereby avoid any of the confusion\nreferenced by the Reviewer.\n In the previous version of the manuscript, we explain co-expression networks in\nthe following terms: “A positive co-expression network can be considered as the\ngenes with direct functional associations with the MDD risk gene, while the\nnegative co-expression network can be conceptualised as the processes\ndownstream of those positively correlated genes”. Prompted by the Reviewer’s\ncomments we deleted these sentences from the revised manuscript. In\nretrospect we agree with the Reviewer, these sentences were misleading and not\nhelpful for the reader. Instead, we describe our methodology emphasizing our\nnetworks consist only of straightforward correlations.\nReviewer Figure 4 Enrichment analysis for the MDD risk genes identified via HMAGMA\nand WGCNA modules generated in the DLPFC of eight different studies, including those\nof different ages.\nReviewer Figure 5 Enrichment analysis of the FADS1 network in the same\npositive\nWGCNA networks as those in Reviewer figure 5. The dashed line represents the\nBonferroni corrected P-value threshold (the same threshold as in Reviewer Figure 5).\nDue to a very high overlap, some points (i.e. those laying along the upper line of the\nplot) had a P-value of 0.\nLine 106-108 of the results: Operationalizing the word “network”\n“Therefore, we use the word “network” to refer to the set of correlations with an MDD\nrisk gene. We also note this definition does not capture correlations between genes that\nare not the MDD risk gene (as is captured in WGCNA networks for instance)”\nLine 93-98 of the results: Softening and caveating our language around WGCNA,\nincluding the references noted by the Reviewer\n“The MDD risk genes from HMAGMA were not enriched in any DLPFC co-expression\nnetworks identified from any of the eight previous studies we assessed 9–16, all of which\nused a WGCNA (weighted correlation network analysis) approach (Supplementary\nFigure 1). Although WGCNA, as shown by many studies, is a powerful approach for\nmechanistic inference 9–18, considering the broad distribution of MDD risk genes across\nmodules we reasoned a seeded co-expression approach may be more informative for\nmechanistic inference in this instance.”\nLine 445-447 of the discussion: Acknowledgment that WGCNA approaches are great\ntools for mechanistic inference, in particular the approach of Borcuk et al\n“Although we use seeded co-expression networks throughout our study, alternative\napproaches using WGCNA networks, like the elegant approach from Borcuk et al, are\nalso powerful tools for mechanistic inference”\nMinor suggestions:\n• In the generation of seeded co-expression networks the number of multiple\ncomparisons used in the FDR correction is not mentioned. Is it 94 or 94 by 23,326?\nThe FDR corrects for 23,326 comparisons for each network. We have added this to the\nmethods section\nLines 499- 500 of the methods\n“The FDR was calculated for each network independently and based on the 23,326\ncorrelations conducted.”\n• In the validation analyses for the co-expression networks, it is not clearly defined what\nare the terms of the “correlation between each of our primary co-expression networks\nand the replication dataset”. I am interpreting that for each seed gene you correlated the\ncorrelation coefficients of the co-expressed genes found in the primary and replication\nanalysis. Is that the case? If you performed the exact same analysis in the replication\ndatasets (R > |0.5| and an FDR p-value < 0.05), what is the overlap of the positively and\nnegatively correlated genes between the primary and replication datasets? Do you find\nthe exact same genes?\nYes indeed, the Reviewer is correct both in their interpretation of our initial validation\nand that the suggested analysis would be a valuable addition to the manuscript. We\npresent the new analysis in Reviewer Figure 2 (which we reproduce below for ease of\nreference) and in Supplementary Fig 4 of the revised manuscript. As seen in Reviewer\nFigure 2, most co-expression networks (90% of the internal validation and 82% of the\nexternal validation) are highly enriched in their corresponding networks in the validation\ndatasets. We also see a similar pattern to our previous validation efforts, where\nnetworks with more co-expressed genes had a stronger validation.\nReviewer Figure 2 Enrichment analyses for the original networks in those derived from\nvalidation datasets. In A) the percentage of networks significantly enriched (or not\nenriched) for genes in their corresponding networks in both the internal and external\nvalidation datasets. In B) and C) the correlation between the odds ratio\n(derived from a Fisher’s exact test) and the number of co-expressed genes in a network\nis presented for the external (B) and internal (C) validations.\nLines 155-157 of the results: Discussion of the new enrichment analysis\n“We also note substantial overlap in the genes within these co-expression networks,\nwith 90% and 82% of the internal and external datasets, respectively, enriched for\ngenes in their corresponding co-expression networks from the original dataset\n(Supplementary Figure 4).”\n• In the Network analysis on the co-expression networks conducted with igraph, it is not\nclear how the weighted adjacency matrix for all networks was constructed. What is\nrepresented in each cell of the adjacency matrix?\nThe adjacency matrix represents the presence or absence of a correlation between the\nfull co-expression networks. To generate the adjacency matrix, we first identified the co-\nexpression networks with more than 1000 positively co-expressed genes (defined as\nthe number of genes with an FDR<0.05 and R>0.5). For these networks we then\ncorrelated the coefficients of all genes, with those from the other networks. Specifically,\nthe correlation matrix was generated from an input matrix containing 23,326 rows (the\nnumber of genes against which co-expression was tested) and 14 columns (the number\nof networks with more than 1000 genes positively co-expressed). The adjacency matrix\nwas then constructed from this correlation matrix by thresholding the presence or\nabsence of a correlation.\nWe thank the Reviewer for this point and have added these words verbatim to the\nmethods section.\nLine 549-556 of the methods section: Description of the adjacency matrix\n“To generate the adjacency matrix, we first identified the co-expression networks with\nmore than 1000 positively co-expressed genes (defined as the number of genes with an\nFDR<0.05 and R>0.5). Then for these networks we correlated the coefficients of all\ngenes, with those from the other networks. Specifically, the correlation matrix was\ngenerated from an input matrix containing 23,326 rows (the number of genes against\nwhich co-expression was tested) and 14 columns (the number of networks with more\nthan 1000 genes positively co-expressed). The adjacency matrix was then constructed\nfrom this correlation matrix by thresholding the presence or absence of a correlation.”\nReferences\n1. Sambra, V., Echeverria, F., Valenzuela, A., Chouinard-Watkins, R. & Valenzuela,\nR. Docosahexaenoic and Arachidonic Acids as Neuroprotective Nutrients\nthroughout the Life Cycle. Nutrients 13, 986 (2021).\n2. Tosi, F., Sartori, F., Guarini, P., Olivieri, O. & Martinelli, N. Delta-5 and Delta-6\nDesaturases: Crucial Enzymes in Polyunsaturated Fatty Acid-Related Pathways\nwith Pleiotropic Influences in Health and Disease. 61–81 (2014) doi:10.1007/978-\n3-319-07320-0_7.\n3. Dyall, S. C. Long-chain omega-3 fatty acids and the brain: A review of the\nindependent and shared effects of EPA, DPA and DHA. Front Aging Neurosci 7,\n125532 (2015).\n4. Parikshak, N. N., Gandal, M. J. & Geschwind, D. H. Systems biology and gene\nnetworks in neurodevelopmental and neurodegenerative disorders. Nature\nReviews Genetics 2015 16:8 16, 441–458 (2015).\n5. Akil, H. et al. Treatment resistant depression: A multi-scale, systems biology\napproach. Neurosci Biobehav Rev 84, 272–288 (2018).\n6. Davyson, E. et al. Metabolomic investigation of major depressive disorder\nidentifies a potentially causal association with polyunsaturated fatty acids. Biol\nPsychiatry 0, (2023).\n7. Levey, D. F. et al. Bi-ancestral depression GWAS in the Million Veteran Program\nand meta-analysis in >1.2 million individuals highlight new therapeutic directions.\nNature Neuroscience 2021 24:7 24, 954–963 (2021).\n8. Lalovic, A., Klempan, T., Sequeira, A., Luheshi, G. & Turecki, G. Altered\nexpression of lipid metabolism and immune response genes in the frontal cortex\nof suicide completers. J Affect Disord 120, 24–31 (2010).\n9. Pergola, G. et al. Prefrontal Coexpression of Schizophrenia Risk Genes Is\nAssociated With Treatment Response in Patients. Biol Psychiatry 86, 45–55\n(2019).\n10. Pergola, G. et al. Consensus molecular environment of schizophrenia risk genes\nin coexpression networks shifting across age and brain regions. Sci Adv 9, (2023).\n11. Fromer, M. et al. Gene expression elucidates functional impact of polygenic risk\nfor schizophrenia. Nature Neuroscience 2016 19:11 19, 1442–1453 (2016).\n12. Radulescu, E. et al. Identification and prioritization of gene sets associated with\nschizophrenia risk by co-expression network analysis in human brain. Mol\nPsychiatry 25, 791–804 (2020).\n13. Walker, R. L. et al. Genetic Control of Expression and Splicing in Developing\nHuman Brain Informs Disease Mechanisms. Cell 179, 750-771.e22 (2019).\n14. Gandal, M. J. et al. Shared molecular neuropathology across major psychiatric\ndisorders parallels polygenic overlap. Science (1979) 359, 693–697 (2018).\n15. Punzi, G. et al. Genetics and Brain Transcriptomics of Completed Suicide.\nhttps://doi.org/10.1176/appi.ajp.2021.21030299 179, 226–241 (2022).\n16. Hartl, C. L. et al. Coexpression network architecture reveals the brain-wide and\nmultiregional basis of disease susceptibility. Nature Neuroscience 2021 24:9 24,\n1313–1323 (2021).\n17. Borcuk, C. et al. Network-wide risk convergence in gene co-expression identifies\nreproducible genetic hubs of schizophrenia risk. Neuron 112, 3551-3566.e6\n(2024).\n18. Pergola, G. et al. DRD2 co-expression network and a related polygenic index\npredict imaging, behavioral and clinical phenotypes linked to schizophrenia.\nTranslational Psychiatry 2017 7:1 7, e1006–e1006 (2017)."}, {"round": 2, "reviewer_comments": {"reviewer_2": "Thank you for taking the time to thoroughly revise the manuscript. After reviewing the updated version, I’m glad to see that\nthe main issues raised in the initial review have been properly addressed. The reanalysis and reinterpretation of the data are\nwell-executed, and the authors have provided comprehensive answers to all my questions and comments.\nThe revisions have noticeably strengthened the manuscript. The added clarifications and adjustments to the methodology\nimprove both the transparency and the overall rigor of the study. The conclusions now feel more solid, as they are better\nsupported by the presented data and analysis.\nIn summary, the authors have responded thoughtfully and effectively to the feedback. I don’t have any further major\nconcerns.", "reviewer_3": "The authors made an effort to address my comments thoughtfully. The reply to the point about H-MAGMA is adequate and\nreasonable. I would suggest writing in the methods section the reasoning behind your choice of using H-MAGMA just as\nclearly as you have written it in the rebuttal. Likewise, the analyses on additional WGCNA networks are interesting, and the\nfailure to identify MDD convergence in WGCNA networks should be reported at least in the supplement. The failure in\nexternal replication in response to Rev #2 points is also important and should be reported.\nAs results currently stand, the Authors have stumbled upon an interesting finding regarding FADS1, but the route taken is\narguably invalid and not replicable. I am deeply sympathetic to the efforts made by the authors to justify their approach, and\nthey have considered the points raised quite thoroughly and transparently. However, I cannot turn a blind eye to the\nmethodological validity.\nIn case of results discrepancy due to methods, we need to understand why results differ. Using betweenness to identify\nnetworks in a thresholded setting (here, Pearson's R > 0.5 and FDR<.05) is a valid way to identify networks, but the authors\nshow that it does not replicate in CMC. They also show that WGCNA networks are enriched for FADS1 partners, so the\nquestion is: what is the difference between WGCNA methods and the seeded approach used here?\nI think the main difference is that this approach does not use an exponent to bring the similarity matrix into scale invariance.\nThe Authors are well aware that they can easily verify this in their own work by correlating their similarity values with log-\ntransformed values. Publications by Horvath and Langfelder elaborate on why WGCNA enforces scale invariance. Briefly,\nthe reason is to obtain a scale-free topology, characterized by a small number of highly connected genes (\"hubs\") and many\ngenes with fewer connections. Specifically, scale invariance ensures the following:\n1. Robustness and Stability. Networks with scale-free topology are robust to random perturbations, e.g., measurement errors.\nNote the lack of replication in CMC.\n2. Biological Realism. Biological networks, including gene co-expression networks, typically exhibit scale-free properties.\n3. Identification of Key Hub Genes. The Authors intend to show that FADS1 is a hub.\n4. Facilitation of Network Comparisons. Scale-free topology simplifies comparisons across datasets or biological conditions\n(e.g., healthy vs. diseased states) because similar structural features (hub genes and connectivity patterns) become clearly\nidentifiable. Note the lack of consistency between published networks and the network here reported, as well as the lack of\nreplication in CMC.\nThe dependence of results on the number of genes included is also bad news, as it suggests that the connections between\na small number of genes are too unreliable. The end result is that the outcome reported here is unlikely to be replicated.\nAgain, I commend the Authors for the work they have put into this finding, which is prospectively interesting, assuming that a\nrigorous methodology can be found to support it."}, "author_response": "Response to Reviewer’s comments\nAstrocyte fatty acid metabolism as a driver of risk for\nmajor depressive disorder\nEamon Fitzgerald, Nicholas O'Toole, Irina Pokhvisneva, Eric J Nestler, Gustavo Turecki,\nCorina Nagy, Michael J Meaney\nPlease find a detailed response to each of the Reviewer’s comments below. We believe\nour manuscript has been greatly improved through this process and we very much\nappreciate the Reviewers’ time and effort.\nThroughout this response:\n The Reviewer’s original comments appear in black and in quotation marks.\n Our response appears in blue\n Text quoted from the revised manuscript appears in red\nReviewer #2 (Remarks to the Author):\n“The reanalysis and reinterpretation of the data are well-executed, and the authors\nhave provided comprehensive answers to all my questions and comments.”\n“The revisions have noticeably strengthened the manuscript. The added clarifications\nand adjustments to the methodology improve both the transparency and the overall\nrigor of the study. The conclusions now feel more solid, as they are better supported by\nthe presented data and analysis.”\n“In summary, the authors have responded thoughtfully and effectively to the feedback. I\ndon’t have any further major concerns.”\nThank you for taking the time to thoroughly revise the manuscript. We greatly\nappreciate the extensive and very constructive commentary on our manuscript. The\npaper has been strengthened considerably as a function of the comments.\nReviewer #3 (Remarks to the Author):\nWe thank the Reviewer once again for their engagement with our manuscript.\n1. “The dependence of results on the number of genes included is also bad news,\nas it suggests that the connections between a small number of genes are too\nunreliable. The end result is that the outcome reported here is unlikely to be\nreplicated. Again, I commend the Authors for the work they have put into this\nfinding, which is prospectively interesting, assuming that a rigorous methodology\ncan be found to support it.”\nWe agree with the Reviewer that the validation of results is of paramount importance.\nThus, in the original version of our manuscript, we validated the FADS1 co-expression\nnetwork in two independent cohorts (see Supplementary Figure 3). NOTEWe\ndemonstrated highly significant concordance between the two validation datasets and\nour initial discovery sample using two validation metrics. For the first metric, we\ncorrelated all of the correlation coefficients from the co-expression network with those\nfrom the validation datasets (r = 0.81 and 0.69; P < 2.2 × 10⁻¹⁶ for both validations). For\nthe second metric, we showed strong enrichment of thresholded genes (i.e., FDR < 0.05\nand R < 0.5) from the validation datasets in the discovery dataset (again, P < 2.2 × 10⁻¹⁶\nfor both validations). These data can be found on Line 187-189 of the manuscript.\nIn this revision, we now present two additional validations. We use datasets from the\ndorsolateral prefrontal cortex of control tissue from a bulk RNA-seq (Labonté et al) and\na microarray (Chang et al) study. In the previous round of revisions, the Reviewer\nsuggested that enrichment of the FADS1 co-expression network in a replicate co-\nexpression network was their preferred metric of validation, which is therefore the one\nwe provide below and in the newly revised manuscript. In total, we now validate the\nFADS1 co-expression network in four distinct cohorts.\n Labonté et al (RNA-seq): P < 2.2 × 10⁻¹⁶, OR = 18.95\n Chang et al (microarray): P < 2.2 × 10⁻¹⁶, OR = 150.69\nThese highly significant replication metrics confirm that the FADS1 co-expression\nnetwork is robust across different gene-expression platforms and sample collections. By\nintegrating multiple orthogonal validations (bulk reanalysis, single-cell and spatial RNA-\nseq, genetic enrichment, and protein–protein interaction data) we provide a rigorous,\nmulti-modal foundation for our findings.\nLines 189-191 of the results:\n“Furthermore, we were able to validate the FADS1 network in two other external\ndatasets using DLPFC tissue from control samples (Labonté et al 46 r = 0.59;\nenrichment P-value <2.22-16 and Chang et al 38, r = 0.52; enrichment P-value <2.22-\n16).”\n2. “The authors made an effort to address my comments thoughtfully. The reply to\nthe point about H-MAGMA is adequate and reasonable. I would suggest writing\nin the methods section the reasoning behind your choice of using H-MAGMA just\nas clearly as you have written it in the rebuttal. Likewise, the analyses on\nadditional WGCNA networks are interesting, and the failure to identify MDD\nconvergence in WGCNA networks should be reported at least in the supplement.\nThe failure in external replication in response to Rev #2 points is also important\nand should be reported.”\nWe thank the Reviewer and have added the requested explanation of H-MAGMA to the\nrevised manuscript and the additional WGCNA analysis to the supplement. Our external\nreplication data was included in the manuscript during the previous round of revision,\nbut we provide it again below for ease of access.\nLine 469-471 of the methods:\n“We chose H-MAGMA over MAGMA because the region-specific approach to\nannotating the GWAS summary statistics aligned nicely with our region-specific\napproach to seeded co-expression networks”\nLine 202-204 of the results:\n“The FADS1 network was enriched in several WGCNA modules from the human\npositive\nDLPFC across numerous studies suggesting at least some congruence between the\ndifferent methods (Supplementary Figure 8B)”\nLine 179-182 of the results:\n“We also observed a similar pattern in the internal GTEx validation dataset, with the\nFADS1 network having the largest betweenness of any network. However, in the\npositive\nexternal CMC validation dataset, the FADS1 network did not have the highest\npositive\nbetweenness (Supplementary Figure 6).”\n“As results currently stand, the Authors have stumbled upon an interesting finding\nregarding FADS1, but the route taken is arguably invalid and not replicable. I am deeply\nsympathetic to the efforts made by the authors to justify their approach, and they have\nconsidered the points raised quite thoroughly and transparently. However, I cannot turn\nan eye to the methodological validity.\nIn case of results discrepancy due to methods, we need to understand why results\ndiffer. Using betweenness to identify networks in a thresholded setting (here, Pearson's\nR > 0.5 and FDR<.05) is a valid way to identify networks, but the authors show that it\ndoes not replicate in CMC. They also show that WGCNA networks are enriched for\nFADS1 partners, so the question is: what is the difference between WGCNA methods\nand the seeded approach used here?”\nOur goal is not to perform a direct head-to-head comparison of WGCNA versus seeded\nco-expression, but rather to leverage a hypothesis-driven framework that starts with\nMDD risk genes (e.g. FADS1) and uncovers their mechanistic partners. Simply, we used\nseeded co-expression analysis in our study for two reasons:\n1. MDD risk genes were not enriched in WGCNA networks\n2. We wanted to infer mechanisms proximal to MDD risk genes, thereby looking at\nthe WGCNA modules in which MDD risk genes were present would amount to\nlooking at all modules and all genes.\nBy anchoring the network on FADS1, we infer gene–gene relationships most likely to\nreflect functional associations of the seed gene, which we then validate across multiple\ndatasets. Through bootstrapping and orthogonal validations, we demonstrate that these\nFADS1 co-expression networks are robust, reproducible, and biologically meaningful.\n“I think the main difference is that this approach does not use an exponent to bring the\nsimilarity matrix into scale invariance. The Authors are well aware that they can easily\nverify this in their own work by correlating their similarity values with log-\ntransformed values. Publications by Horvath and Langfelder elaborate on why WGCNA\nenforces scale invariance. Briefly, the reason is to obtain a scale-free topology,\ncharacterized by a small number of highly connected genes (\"hubs\") and many\ngenes with fewer connections. Specifically, scale invariance ensures the following:\n1. Robustness and Stability. Networks with scale-free topology are robust to random\nperturbations, e.g., measurement errors. Note the lack of replication in CMC.\n2. Biological Realism. Biological networks, including gene co-expression networks,\ntypically exhibit scale-free properties.\n3. Identification of Key Hub Genes. The Authors intend to show that FADS1 is a hub.\n4. Facilitation of Network Comparisons. Scale-free topology simplifies comparisons\nacross datasets or biological conditions (e.g., healthy vs. diseased states) because\nsimilar structural features (hub genes and connectivity patterns) become clearly\nidentifiable. Note the lack of consistency between published networks and the network\nhere reported, as well as the lack of replication in CMC. The dependence of results on\nthe number of genes included is also bad news, as it suggests that the connections\nbetween a small number of genes are too unreliable. The end result is that the outcome\nreported here is unlikely to be replicated. Again, I commend the Authors for the work\nthey have put into this finding, which is prospectively interesting, assuming that a\nrigorous methodology can be found to support it.”\nWe agree with the Reviewer, sound results that reproduce are of the utmost importance.\nAs mentioned earlier we have now shown the FADS1 co-expression network validates\nrobustly in four independent cohorts. We also show that the FADS1 network is highly\nenriched in protein–protein interactions (Supplementary Figure 7), which acts as\nanother form of validation. We have additionally validated our single-nucleus RNA-seq\nfindings in an independent cohort (Supplementary Figures 10 and 11). These single-\nnucleus RNA-seq analyses build upon the FADS1 co-expression network, and therefore\ntheir reproducibility again speaks to the integrity of our findings. We must therefore,\nrespectfully, disagree with the Reviewer when they say our results are “unlikely to be\nreplicated.”\nWe fully acknowledge the merits of WGCNA in transcriptomic research. Indeed, we\nhave previously published transcriptomic analyses centered on the use of WGCNA\n(Fitzgerald et al., Nature Communications, 2023; Arcego et al., Biological Psychiatry,\n2023). However, in the current study, seeded co-expression networks offered a more\ndirect path to mechanistic insights into major depressive disorder (MDD). Using a\nseeded approach enabled us to investigate functionally coherent gene networks\ncentered on biologically relevant seed genes, such as FADS1.\nWe also note that we are not alone in applying this methodology. Seeded co-expression\nnetworks have been used effectively in several high-profile publications, including:\nWillsey et al. (Cell, 2013): Coexpression Networks Implicate Human Midfetal Deep\nCortical Projection Neurons in the Pathogenesis of Autism; Kim et al. (Nature\nNeuroscience, 2021): Brain gene co-expression networks link complement signaling\nwith convergent synaptic pathology in schizophrenia; Liao et al. (Cell Genomics, 2023):\nConvergent coexpression of autism-associated genes suggests some novel risk genes\nmay not be detectable in large-scale genetic studies."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors have made commendable efforts in their revisions, including validating the FADS1 co-expression network\nacross four independent cohorts and providing additional support from protein–protein interactions and single-nucleus RNA-\nseq.\nThese steps clearly strengthen the manuscript.\nNonetheless, I find Reviewer #3’s concerns to be fair and justified.\nThe key issue is not WGCNA versus seeded networks alone, but the need to mitigate correlation bias, demonstrate\nsensitivity and specificity, and show convergence across unbiased methods.\nThe lack of replication in the high-quality CMC dataset underscores this point.\nIn my view, further addressing these methodological issues is both feasible and necessary to reinforce the robustness of the\nstudy’s conclusions.", "reviewer_3": "The authors have basically reiterated their points regarding the technique they used, which seems to be the only one\navailable to obtain the desired results. Correlation alone is not a valid parameter for replication, as shown in the following\nexample: suppose a dataset shows a positive correlation bias (as is usually the case, see Freytag et al., 2015); now, you can\nisolate the desired correlation and easily claim it replicates in other datasets, but that is not necessarily the case due to a\nbiological reality, as instead derives from inflated correlation across the board. Sensitivity and specificity are key. This is the\nreason why other authors have taken a number of steps to reduce artifacts, including ensuring scale invariance, adopting\npreservation statistics, controlling for gene length and GC content within permutations, using positive and negative gene\nsets when removing unwanted variation, and so on. The authors are fully aware of the literature I am citing, therefore, I see\nno need to point to individual papers.\nThere is NO replication in CMC, which is the best quality and best powered RNAseq dataset considered (microarray\ncorrelation bias is notoriously even higher).\nThe attempt to deflect criticism with the argument that others have followed a given path is remarkably weak, because peer\nreview is concerned with this paper, and if a paper published in Cell back in 2013 did as the authors do now, and over the\nlast 12 years these seminal techniques have been replaced by more reliable ones, that is simply the advancement of\nscience. Laboratory protocols and standards evolve; therefore, I do not see that as a valid response.\nThis is not an issue of WGCNA vs something else, as claimed by the authors, but of removing artifacts and demonstrating\nmethod convergence. If the authors used, say, Sparse Decomposition of Arrays (Hore et al., 2016), which is not a graph\ntechnique, or several other VALID and UNBIASED techniques to reach their desired result, then I would rethink my position,\nas perhaps WGCNA is not always the best, given that the convergence of separate methods onto the same scientific result\nis an essential aspect of true findings.\nSimply insisting on using a notoriously biased approach will not convince me that the approach is unbiased, sorry. What I\nhave learned from this paper so far is that WGCNA modules are not enriched for MDD hits, which, for me, calls into question\nwhether MDD, as clinically defined, is a diagnostic entity biologically comparable to schizophrenia, where the evidence that\nGWAS results converge into biological processes is overwhelming. If that is the case, should we even bother to study the\ntranscriptomics of such a heterogeneous and environmentally driven group of conditions without questioning the clinical\nlabel?", "reviewer_4": "I do not share the concerns that reviewer 3 has with the network analysis. Network analysis is always subjective, and the\nauthors do a good effort to replicate, plus they explain the pros and cons of their choice. You can’t do any better: the rest is a\nmatter of taste.\nI am more worried about the relation they find between gene expression and MDD (which is key for this paper). The\nassociation between gene expression and MDD status can be a symptom of, or a cause for depression, or causality can be\nbidirectional. Starting with GWAS genes identified in GWAS helps a bit (because DNA is causal) but they go from these\ngenes to a network of genes from which we have no clue about their causal involvement with MDD. Also, their finding is an\nincreased networked structure/connection in controls compared to cases. This can also be due to lifestyle variables like\nsmoking or BMI, which is increased in MDD, and not taken into account in the analysis. These two variable have an\nenormous impact on gene expression (and probably their correlations). I therefore find the conclusion “PPAR is a putative\ntherapeutic target in MDD” way to strong!\nAnyway, a comment in the discussion about this weakness would suffice (as far as I am concerned). Overall the authors do\nan enormous overwhelming amount of analysis but the methods lack details or may be not done appropriately (like the\ncovariates used in the association analysis with MDD (smoking? BMI?)"}, "author_response": "Response to Reviewer’s comments\nAstrocyte fatty acid metabolism as a driver of risk for\nmajor depressive disorder\nEamon Fitzgerald, Nicholas O'Toole, Irina Pokhvisneva, Eric J Nestler, Gustavo Turecki,\nCorina Nagy, Michael J Meaney\nPlease find a detailed response to each of the Reviewer’s comments below. We believe\nour manuscript has been greatly improved through this process and we very much\nappreciate the Reviewers’ time and effort.\nThroughout this response:\n The Reviewer’s original comments appear in black and in quotation marks.\n Our response appears in blue\n Text quoted from the revised manuscript appears in red\nREVIEWER COMMENTS\nReviewer #2 (Remarks to the Author):\nThe authors have made commendable efforts in their revisions, including validating the\nFADS1 co-expression network across four independent cohorts and providing additional\nsupport from protein–protein interactions and single-nucleus RNA-seq.\nThese steps clearly strengthen the manuscript.\nNonetheless, I find Reviewer #3’s concerns to be fair and justified.\nThe key issue is not WGCNA versus seeded networks alone, but the need to mitigate\ncorrelation bias, demonstrate sensitivity and specificity, and show convergence across\nunbiased methods.\nThe lack of replication in the high-quality CMC dataset underscores this point.\nIn my view, further addressing these methodological issues is both feasible and\nnecessary to reinforce the robustness of the study’s conclusions.\nWe thank the Reviewer for their constructive comments and understand their concerns\nregarding potential methodological bias introduced by the seeded co-expression\nnetwork approach.\nWe first note data that were included in a response to Reviewer 3 during a previous\nround of review (and which had previously been included as Supplementary Figure 8B\nof the revised manuscript) that may also be relevant here (see Reviewer Figure 1\nbelow). These data show enrichment of the FADS1 network in WGCNA modules from\neight independent studies that used WGCNA in human post-mortem DLPFC tissue.\nWe believe this analysis demonstrates that the seeded coexpression approach captures a\ncomparable signature to that identified by WGCNA, while also enabling the generation\nof networks centered on specific genes of interest, which is a feature particularly\nrelevant to the goals of our analysis.\nWe have also conducted new analyses to directly address the Reviewer’s comment\nregarding the “need to mitigate correlation bias, demonstrate sensitivity and specificity,\nand show convergence across unbiased methods”. We used sparse partial least\nsquares (sPLS) analysis to derive a latent variable of gene expression associated with a\nsingle gene (i.e. FADS1). The sPLS approach has the distinct advantage of addressing\na similar objective to the seeded coexpression approach, while relying on distinct\nstatistical assumptions. We compared the genes identified by sPLS with those obtained\nfrom the seeded coexpression approach and observed a striking degree of overlap:\n96% of genes (909 out of 949) identified by sPLS were also captured within the seeded\ncoexpression network. These data are presented below in Reviewer Figure 2 and in\nSupplementary figure 8C of the revised manuscript.\nIn summary, these results demonstrate the seeded co-expression approach exhibits\nstrong convergence with multiple independent methodologies, each grounded in\ndifferent statistical frameworks.\nReviewer Figure 1 The genes with the FADS1 positive coexpression network were\ntested for enrichment amount WGCNA modules created in eight different published\nstudies of the postmortem human DLPFC. We observed very strong enrichment of the\nFADS1 network in modules of each study. The dashed line indicates the threshold for\nmultiple comparisons and each study has at least one network exceeding this threshold.\nReviewer Figure 2 Venn diagram showing the overlap in genes showing a positive\nassociation with FADS1 in the DLFPC. The seeded co-expression (white) is the network\npresented throughout the manuscript. The sPLS (sparse partial least squares) is an\napproach building a latent variable of DLPFC genes aiming to explain FADS1\nexpression, the result being a latent variable associated of 949 genes that is associated\nwith FADS1 expression (r = 0.86 and p-value < 2.2e-16 for the Pearson correlation\nbetween the latent variable and FADS1 expression). Of these 949 genes, 909 are\nshared with the seeded co-expression network for FADS1, which is highly statistically\nsignificant by Fisher’s Exact Test (shown underneath the Venn Diagram).\nLines 204 – 209 of the results section\n“The FADS1 network was enriched in several WGCNA modules from the human\npositive\nDLPFC across numerous studies (Supplementary Fig. 8B). We also found a signature\nhighly overlapping to the FADS1 network when we used sparse partial least\npositive\nsquares to derive a latent variable of DLPFC gene expression correlated with FADS1\nexpression (909 of 949 genes overlapped with the FADS1 network; Supplementary\npositive\nFig. 8C). Together these data suggest the seeded co-expression approach is\nreproducible with other distinct methods.”\nReviewer #3 (Remarks to the Author):\nThe authors have basically reiterated their points regarding the technique they used,\nwhich seems to be the only one available to obtain the desired results. Correlation alone\nis not a valid parameter for replication, as shown in the following example: suppose a\ndataset shows a positive correlation bias (as is usually the case, see Freytag et al.,\n2015); now, you can isolate the desired correlation and easily claim it replicates in other\ndatasets, but that is not necessarily the case due to a biological reality, as instead\nderives from inflated correlation across the board. Sensitivity and specificity are key.\nThis is the reason why other authors have taken a number of steps to reduce artifacts,\nincluding ensuring scale invariance, adopting preservation statistics, controlling for gene\nlength and GC content within permutations, using positive and negative gene sets when\nremoving unwanted variation, and so on. The authors are fully aware of the literature I\nam citing, therefore, I see no need to point to individual papers.\nThere is NO replication in CMC, which is the best quality and best powered RNAseq\ndataset considered (microarray correlation bias is notoriously even higher).\nThe attempt to deflect criticism with the argument that others have followed a given path\nis remarkably weak, because peer review is concerned with this paper, and if a paper\npublished in Cell back in 2013 did as the authors do now, and over the last 12 years\nthese seminal techniques have been replaced by more reliable ones, that is simply the\nadvancement of science. Laboratory protocols and standards evolve; therefore, I do not\nsee that as a valid response.\nThis is not an issue of WGCNA vs something else, as claimed by the authors, but of\nremoving artifacts and demonstrating method convergence. If the authors used, say,\nSparse Decomposition of Arrays (Hore et al., 2016), which is not a graph technique, or\nseveral other VALID and UNBIASED techniques to reach their desired result, then I\nwould rethink my position, as perhaps WGCNA is not always the best, given that the\nconvergence of separate methods onto the same scientific result is an essential aspect\nof true findings.\nSimply insisting on using a notoriously biased approach will not convince me that the\napproach is unbiased, sorry. What I have learned from this paper so far is that WGCNA\nmodules are not enriched for MDD hits, which, for me, calls into question whether MDD,\nas clinically defined, is a diagnostic entity biologically comparable to schizophrenia,\nwhere the evidence that GWAS results converge into biological processes is\noverwhelming. If that is the case, should we even bother to study the transcriptomics of\nsuch a heterogeneous and environmentally driven group of conditions without\nquestioning the clinical label?\nWhat I have learned from this paper so far is that WGCNA modules are not enriched for\nMDD hits, which, for me, calls into question whether MDD, as clinically defined, is a\ndiagnostic entity biologically comparable to schizophrenia, where the evidence that\nGWAS results converge into biological processes is overwhelming. If that is the case,\nshould we even bother to study the transcriptomics of such a heterogeneous and\nenvironmentally driven group of conditions without questioning the clinical label?\nWe first note data included in a previous response (shown in Supplementary Figure 8B\nof the revised manuscript and reproduced above in Reviewer Figure 1). These data\ndemonstrate enrichment of the FADS1 network within networks from eight independent\nstudies that used WGCNA to generate modules from human post-mortem DLPFC\ntissue-studies originally suggested by the Reviewer. The indicated analysis shows that\nthe seeded coexpression approach captures a signature highly comparable to that\nobtained within WGCNA modules from eight different studies (we also note the\nReviewer suggested these studies).\nTo further address the Reviewer’s concern, we conducted an additional analysis to test\nthe convergence between our seeded coexpression methodology and networks derived\nusing an orthogonal approach. We selected sparse partial least squares (sPLS) because\nit allows the construction of a latent variable representing gene expression associated\nwith a single gene (e.g., FADS1), thereby achieving a similar objective to seeded\ncoexpression while relying on distinct statistical assumptions. The overlap between\ngenes identified by the two approaches is shown above in Reviewer Figure 2. Notably,\n96% of genes (909 out of 949) identified by sPLS were also captured in our seeded\ncoexpression network, demonstrating a high degree of convergence.\nIn summary, these findings show that the seeded coexpression approach exhibits\nstrong agreement with multiple independent methodologies that employ distinct\nstatistical frameworks.\nReviewer #4 (Remarks to the Author):\nI do not share the concerns that reviewer 3 has with the network analysis. Network\nanalysis is always subjective, and the authors do a good effort to replicate, plus they\nexplain the pros and cons of their choice. You can’t do any better: the rest is a matter of\ntaste.\nI am more worried about the relation they find between gene expression and MDD\n(which is key for this paper). The association between gene expression and MDD status\ncan be a symptom of, or a cause for depression, or causality can be bidirectional.\nStarting with GWAS genes identified in GWAS helps a bit (because DNA is causal) but\nthey go from these genes to a network of genes from which we have no clue about their\ncausal involvement with MDD. Also, their finding is an increased networked\nstructure/connection in controls compared to cases. This can also be due to lifestyle\nvariables like smoking or BMI, which is increased in MDD, and not taken into account in\nthe analysis. These two variable have an enormous impact on gene expression (and\nprobably their correlations). I therefore find the conclusion “PPAR is a putative\ntherapeutic target in MDD” way to strong!\nAnyway, a comment in the discussion about this weakness would suffice (as far as I am\nconcerned). Overall the authors do an enormous overwhelming amount of analysis but\nthe methods lack details or may be not done appropriately (like the covariates used in\nthe association analysis with MDD (smoking? BMI?)\nWe thank the Reviewer for their thoughtful comments and the time spent evaluating our\nmanuscript.\nWe have revised the Methods to specify that we primarily adjusted for technical\ncovariates, choosing to minimally correct for non-technical variables (such as smoking\nand BMI) because they may have causal relations with depressive phenotypes, and\nregressing them could have removed salient biological signal from the data. We agree\nthat unmeasured factors, such as smoking and BMI, may influence gene expression\nand potentially contribute to the observed differences. To address this in the manuscript,\nand in accordance with the suggestion of the reviewer, we have adopted more cautious\nwording around causality and clarified in the Discussion that unmeasured variables\ncould influence our findings.\nWe also expanded our Methods section to be more comprehensive, with an example\nprovided below.\nWe appreciate the Reviewer’s overall positive assessment of our paper, in particular the\nnote acknowledging replication, and thank them again for their time and effort.\nIntroduction lines 85-86\n“Finally, we identified PPAR (peroxisome proliferator activated receptor alpha) as a\npromising therapeutic candidate for further investigation in MDD.”\nDiscussion lines 415- 416\n“Our data suggest that PPARα represents a possible avenue for future therapeutic\ninvestigation in MDD.”\nDiscussion lines 452-454\n“..our findings are largely derived from observational studies and, despite grounding in\ngenetics, unmeasured variables that correlate with MDD status could influence our\nfindings”\nMethods lines 479-497\n“GTEx v8 data were used because of their well validated RNA extraction and\nsequencing methodology, large sample size and well documented technical phenotype\ndata. Raw gene-level read counts, transcript-level expected counts, sample attributes,\nand sequencing metrics were downloaded from the GTEx portal and dbGaP (accession\nnumber phs000424.v9.p2). Only samples and genes passing GTEx QC were retained.\nAnalyses were restricted to Brain – Frontal Cortex (BA9). Samples were excluded if\nmarked ineligible or if the donor had the following: amyotrophic lateral sclerosis,\nAlzheimer’s or other dementia, encephalitis, influenza, multiple sclerosis, Parkinson’s\ndisease, stroke or cerebrovascular disorder, systemic lupus, bacterial or viral sepsis,\nactive cancer, schizophrenia, major depressive disorder, substance use (heroin or\nprescription drug abuse), pneumonia, unexplained seizures, or positive blood cultures.\nAfter filtering, 111 samples remained. Genes with >0.1 counts per million (CPM) in\n>25% of samples were retained (23,326 genes). Counts were normalized using the\nTMM method in edgeR v3.38.1, log₂-transformed, and filtered for low expression.\nOutlier samples were identified by connectivity Z-scores from biweight midcorrelations\n(|Z| > 3) and removed (n = 4). Thirty-two sequencing metrics (e.g., mapping rate, read\nquality, duplication rate, insert size) were standardized and summarized by PCA; the\nfirst ten components were used as covariates. Additional covariates included RNA\nintegrity number (RIN), Hardy scale, refrigeration interval, post-mortem interval\ncategory, sex, and age. We chose to primarily adjust for technical covariates, and to\nminimally correct for non-technical variables (such as smoking and BMI) because they\nmay have causal relationships with depressive phenotypes and regressing them could\nhave removed salient biological signal from the data. Expression residuals were\nobtained by ordinary least squares regression of normalized counts on all covariates.\nThe resulting residualized expression matrix was used for all downstream co-expression\nanalyses.”"}, {"round": 4, "reviewer_comments": {"reviewer_2": "The revised manuscript incorporates several helpful additions, including evidence of FADS1-positive network enrichment\nacross multiple WGCNA-derived DLPFC studies and an sPLS analysis showing strong overlap with the seeded co-\nexpression network. These updates strengthen the overall context of the work and offer useful supporting evidence for\naspects of the proposed model. However, several central methodological concerns raised in the previous review remain\ninsufficiently addressed. Most importantly, the manuscript still does not discuss the lack of replication of the seeded co-\nexpression network in the high-quality CMC dataset, which is essential for evaluating the robustness of the findings. In\naddition, concerns regarding several other methodological issues have only been partially addressed, and these points\nremain not fully resolved. Providing clearer clarification of these issues and acknowledging them more explicitly as\nlimitations would substantially improve the reliability and interpretability of the study’s conclusions.", "reviewer_4": "No further comments."}, "author_response": "Response to Reviewer’s comments\nAstrocyte fatty acid metabolism as a driver of risk for\nmajor depressive disorder\nEamon Fitzgerald, Nicholas O'Toole, Irina Pokhvisneva, Eric J Nestler, Gustavo Turecki,\nCorina Nagy, Michael J Meaney\nPlease find a detailed response to each of the Reviewer’s comments below. We believe\nour manuscript has been greatly improved through this process and we very much\nappreciate the Reviewers’ time and effort.\nThroughout this response:\n The Reviewer’s original comments appear in black and in quotation marks.\n Our response appears in blue\n Text quoted from the revised manuscript appears in red\nREVIEWER COMMENTS\nReviewer #2 (Remarks to the Author):\nThe revised manuscript incorporates several helpful additions, including evidence of\nFADS1-positive network enrichment across multiple WGCNA-derived DLPFC studies\nand an sPLS analysis showing strong overlap with the seeded co-expression network.\nThese updates strengthen the overall context of the work and offer useful supporting\nevidence for aspects of the proposed model. However, several central methodological\nconcerns raised in the previous review remain insufficiently addressed. Most\nimportantly, the manuscript still does not discuss the lack of replication of the seeded\nco-expression network in the high-quality CMC dataset, which is essential for evaluating\nthe robustness of the findings. In addition, concerns regarding several other\nmethodological issues have only been partially addressed, and these points remain not\nfully resolved. Providing clearer clarification of these issues and acknowledging them\nmore explicitly as limitations would substantially improve the reliability and\ninterpretability of the study’s conclusions.\nWe thank the Reviewer once more for their time in reviewing our manuscript.\nWe would like to start with a clarification, the FADS1 network does replicate in the CMC\ndataset, as it does in a total of four independent datasets. We also note that it replicates\nacross 2 different metrics (Pearson’s correlation and Fisher’s exact test). We present\nthese results below for the Reviewer and also on lines 187-192 of the revised\nmanuscript.\nWhat did not replicate in the CMC dataset was the network analysis. We note the\nReviewer’s concern with this lack of replication and we have added text acknowledging\nthis to the discussion section, which can be found below as well as lines 456-457 of the\nrevised manuscript.\nWe also offer a comprehensive limitations section in the discussion to acknowledge the\nissues to which the reviewer is referring.\nLines 187-192 of the results, noting the replication of the FADS1 cohort:\n“the FADS1 network was among the most well-replicated networks in our internal (r\npositive\n= 0.89; enrichment P-value <2.22-16) and CMC external validations (r = 0.61;\nenrichment P-value <2.22-16). Furthermore, we were able to validate the FADS1\nnetwork in two other external datasets Furthermore, we were able to validate the\nFADS1 network in two other external datasets using DLPFC tissue from control samples\n(Labonté et al 46 r = 0.59; enrichment P-value <2.22-16 and Chang et al 38, r = 0.52;\nenrichment P-value <2.22-16).”\nLines 456-457 of the discussion:\n“we note the network analysis we conducted did not replicate in the CMC dataset and\nthus must be interpreted with caution”\nLines 444-459 of the discussion:\n“The current study also comes with limitations. First, the generation of robust co-\nexpression networks requires large sample sizes and study designs, studies such as\nours, will benefit from increasing MDD cohort sizes. Second, we also relied on\ntranscriptomic methods throughout this study. However, these methods are limited in\ntheir ability to infer events that may not be well represented at the transcriptome level,\nsuch as lipid metabolism. Emerging single-cell metabolomic methods promise to build\nupon our results. Third, there are many methods to annotate GWAS to the gene level\nwith no accepted best approach. A more holistic understanding of the genetic\narchitecture of disease will allow for more accurate identification of causal genes from\nGWAS loci and facilitate their interpretation through studies such as ours. Fourth, our\nfindings are largely derived from observational studies and, despite grounding in\ngenetics, unmeasured variables that correlate with MDD status could influence our\nfindings. Finally, the findings emerging from these analyses are intended to identify\nplausible novel mechanisms and therapeutic targets, but require validation.\nFurthermore, we note the network analysis we conducted did not replicate in the CMC\ndataset and thus must be interpreted with this in mind. Although we use seeded co-\nexpression networks throughout our study, alternative approaches using WGCNA\nnetworks, like the elegant approach from Borcuk et al, are also powerful tools for\nmechanistic inference.”"}]

Decision letter