Network preservation reveals shared and unique biological processes associated with chronic alcohol abuse in NAc and PFC.

Chronic alcohol abuse has been linked to the disruption of executive function and allostatic conditioning of reward response dysregulation in the mesocorticolimbic pathway (MCL). Here, we analyzed genome-wide mRNA and miRNA expression from matched cases with alcohol dependence (AD) and controls (n = 35) via gene network analysis to identify unique and shared biological processes dysregulated in the prefrontal cortex (PFC) and nucleus accumbens (NAc). We further investigated potential mRNA/miRNA interactions at the network and individual gene expression levels to identify the neurobiological mechanisms underlying AD in the brain. By using genotyped and imputed SNP data, we identified expression quantitative trait loci (eQTL) uncovering potential genetic regulatory elements for gene networks associated with AD. At a Bonferroni corrected p≤0.05, we identified significant mRNA (NAc = 6; PFC = 3) and miRNA (NAc = 3; PFC = 2) AD modules. The gene-set enrichment analyses revealed modules preserved between PFC and NAc to be enriched for immune response processes, whereas genes involved in cellular morphogenesis/localization and cilia-based cell projection were enriched in NAc modules only. At a Bonferroni corrected p≤0.05, we identified significant mRNA/miRNA network module correlations (NAc = 6; PFC = 4), which at an individual transcript level implicated miR-449a/b as potential regulators for cellular morphogenesis/localization in NAc. Finally, we identified eQTLs (NAc: mRNA = 37, miRNA = 9; PFC: mRNA = 17, miRNA = 16) which potentially mediate alcohol's effect in a brain region-specific manner. Our study highlights the neurotoxic effects of chronic alcohol abuse as well as brain region specific molecular changes that may impact the development of alcohol addiction.

Introduction

Alcohol use disorder (AUD) is a debilitating psychiatric illness with negative health, economic, and social consequences for nearly 15.1 million affected adults worldwide [1]. AUD risk is dependent upon both genetic and environmental factors, with a heritability of 0.49 [2]. The neurobiological framework for understanding how benign, recreational alcohol use leads to AUD follows various hypotheses [3-5], with the most commonly accepted being the cyclical model of addiction [6]. This hypothesis provides valuable insight into the functional specialization of different brain regions that underlie behavioral maladaptations associated with AUD [7]. However, the genetic architecture and molecular mechanisms contributing to alcohol-facilitated neuroadaptations remain widely unknown.

Postmortem brain studies provide the unique opportunity to interrogate neurobiological changes associated with addiction across brain regions and neural pathways [8, 9]. Among these, the mesocorticolimbic system (MCL), which connects the ventral tegmental area (VTA) to the prefrontal cortex (PFC), and nucleus accumbens (NAc), has proven especially sensitive to alcohol-associated neuroadaptations [10-12]. Recent postmortem brain studies of AUD have focused on examining gene and microRNA (miRNA) expression as the biological intermediate between genetic variation and molecular function [13-19]. Studying mRNA and miRNA interactions may also reveal functional relationships that mediate the differential expression of risk AUD genes based on the role miRNAs play in the destabilization and degradation of their target genes [20]. While single gene expression differences are continuously explored, network approaches, such as weighted gene co-expression network analysis (WGCNA), allows genes with correlated expression, and therefore likely related functions, to cluster into modules that then can be analyzed to identify dysregulated biological processes and molecular pathways associated with AUD [21]. Others and we have successfully implemented this method to identify gene networks associated with AUD within the MCL and other brain regions [16, 18]. While postmortem brain expression differences alone are insufficient to infer a causal relationship between AUD and neurobiological function, the integration of genetic information via expression quantitative trait loci (eQTL) analysis can help elucidate the regulatory mechanisms by which genetic variants associated with AUD impact gene expression [22].

Thus, in this study, we seek to expand upon previous research by jointly analyzing two key MCL areas, the NAc and PFC, to identify unique and shared neurobiological processes associated with alcohol dependence (AD). To achieve this, we utilize a case/control study design to identify genes and co-expressed gene networks associated with AD. We then performed a network preservation analysis to determine how well significant modules and their respective biological processes are conserved between the PFC and NAc of chronic alcohol abusers. Within the significant modules, we identified the most connected genes (termed hubs), which were then integrated with miRNA expression data analyzed using the same methodological framework. Finally, we assessed the genetic factors that might impact the functions of risk AD genes via expression quantitative trait loci (eQTL). The miRNA and eQTL analyses were performed in order to identify the regulatory mechanisms by which gene networks identified in PFC and NAc contribute to alcohol addiction.

Materials and methods

Tissue processing and RNA extraction

Postmortem brain tissue from 41 AD cases and 41 controls was provided by the Australian Brain Donor Programs of New South Wales Tissue Resource Centre (NSW TRC) under the support of The University of Sydney, National Health and Medical Research Council of Australia, Schizophrenia Research Institute, National Institute of Alcohol Abuse and Alcoholism, and the New South Wales Department of Health [8]. Samples were excluded based on: (1) history of infectious disease, (2) circumstances surrounding death, (3) substantial brain damage, and (4) post-mortem interval > 48 hours. Total RNA was isolated from PFC (the superior frontal gyrus) and NAc tissue using the mirVANA-PARIS kit (Life Technologies, Carlsbad, CA) following the manufacturer’s suggested protocol. RNA concentrations and integrity (RIN) were assessed via Quant-iT Broad Range RNA Assay kit (Life Technologies) and Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA) respectively. Samples were matched for RIN, age, sex (all male), ethnicity, brain pH, and PMI as part of a previous study [18] yielding a total of 18 case-control matched pairs (n = 36). Due to our matching, the RINs in PFC were slightly lower (mean = 4.5, ±2.04) compared to NAc (mean = 6.9, ±0.84). Previous reports, however, have demonstrated that in post-mortem brain studies reliable results are readily obtained even with RINs ≤4 [23]. For demographic information see S1 Table.

Gene expression microarray and data normalization

Gene expression was assayed using Affymetrix GeneChip Human Genome U133A 2.0 (HG-U133A 2.0) on 22,214 probe sets spanning ~ 18,400 mRNA transcripts, and the Affymetrix GeneChip miRNA 3.0 microarray interrogating the expression of 1733 mature miRNAs as previously described [24]. None of the mRNA or miRNA probes were excluded based on quality control criteria outlined in previous studies [18]. Raw probe data were GCRMA background corrected, log2 transformed, and quantile normalized using Partek Genomics Suite v6.23 (PGS; Partek Inc., St. Louis, MO) to obtain relative gene expression values. A principal component analysis was used to identify potential outlier samples. Only one case sample was removed from the analyses, leaving 18 controls and 17 cases (n = 35) for both brain regions. It has become widely accepted to verify a subset of microarray-generated gene expression changes via an independent platform such as qPCR. Considering limited tissue availability and our extensive use of the Affymetrix platform in the past, we did not include microarray validation in this study which is similar to what other groups have done in the past [25]. We have previously ‘validated’ the same array and platform in independent qPCR experiments with a concordance between microarray and qPCR platforms exceeding 80% in the past [18].

Analysis of differential gene expression

The relationship between AD case status and gene expression in PFC and NAc was analyzed via bidirectional stepwise regression for each gene. This approach is better suited to adjust for the confounding effect of covariates within each transcript’s regression model than the robust linear regression approach employed previously in the analyses of NAc [18]. Regression coefficients were calculated in RStudio (ver. 1.66) using the Stats package (ver. 3.5.1). We further observed that brain pH, RIN, and neuropathology were the most influential covariates in the analyses of NAc expression data, while RIN and smoking history were the two most important covariates in the PFC expression analysis. Finally, we assessed proportion of variance explained by each covariate via the variancePartition package (ver. 1.20) [26]. For more details on our bidirectional stepwise regression, see S1 File.

Network analyses

WGCNA was performed using the WGCNA package in RStudio (ver. 1.66). All nominally significant genes (p≤0.05) were used to generate a signed similarity matrix via pair-wise Pearson correlations. The nominal significance was chosen to (1) include genes with smaller effect sizes, albeit true positive signals, (2) exclude genes with low disease variance, i.e., likely not associated with AD and (3) to provide a sufficient number of genes for the network analysis. WGCNA was performed as outlined previously by us and others [16, 18] and in S1 File. Module eigengenes (MEs), serving as a single aggregate expression value for each of the modules, were correlated to AD case-status and available demographic/biological covariates. To validate WGCNA module clustering, we performed a bootstrap based resampling of 100 iterations with replacement. Next, using WGCNA with the clusterRepro (ver. 0.9) package in RStudio, we identified the level of module preservation between the PFC and NAc by comparing adjacency matrices and calculating the composite preservation statistic (Zsummary). A Zsummary >10 indicates strong evidence for network preservation, Zsummary <10 >2 indicates weak evidence of network preservation and Zsummary <2 indicates no module preservation, as outlined previously [27].

Gene set enrichment analysis

Gene set enrichment was performed using ShinyGo (ver. 0.61) gene annotation database [28]. Gene lists from the significant AD modules from NAc and PFC were enriched using GO biological processes consisting of 15,796 gene sets from the Ensembl BioMart release 96; all p-values for significantly enriched gene sets are FDR adjusted (FDR of 5%). We further performed cell type enrichment using the “userListEnrichment” option within the WGCNA package in R (ver. 1.66) as previously described [18]. Statistical significance of brain-list enrichment was determined via a hypergeometric test; all p-values were adjusted at FDR of 5%.

Hub gene prioritization

Hub genes were defined based on the strength of intramodular connectedness, (also referred as module membership (MM)) calculated from the absolute value of the Pearson’s correlation coefficient between module eigengene and expression values. Hub genes were prioritized for downstream analysis based on MM of r≥0.80 and a significant gene correlation with AD (at p ≤0.05).

eQTL analysis and GWAS/GTEx enrichment

DNA from the postmortem brain sample was processed and genotyped as part of a larger GWAS study [18]. Genotypes with excessive missingness (greater than 20%) and monomorphic for homozygous major and minor alleles were removed. We then selected only, local, cis-eQTLs, defined as SNPs 500kb from the start/stop positions for each hub gene. Such selected SNPs were pruned with Plink v1.9 to exclude variants in linkage disequilibrium (R2 ≥0.7). For eQTL detection, SNP effect on hub gene expression was analyzed via MatrixEQTL package (ver. 2.2) in R using a linear regression model adjusting for covariates. To identify potential disease risk eQTLs, we further tested for an interaction (SNP x AD) term between genotype and AD status using the “modelLINEAR_CROSS” argument. A significant genotype/disease interaction for a SNP/gene pair would indicate that the effect of genotype on expression is significantly different in AD cases versus controls. To determine the overlap between the eQTLs in our sample (at p≤0.002) and significant GWAS hits (at p≤1E-4) from previously reported alcohol and smoking GWAS [29, 30], we employed the Simes enrichment test [31]. We further tested the overlap between eQTLs obtained from our analyses against eQTLs obtained from GTEx consortium [32]. The significance of this overlap was assessed via a Fisher’s exact test at p ≤0.05 threshold. See S1 File for more details.

MiRNA/mRNA target prediction

The relationship between significant miRNA and mRNA modules from each brain region was examined by performing a Pearson’s correlation on the miRNA and mRNA module MEs using the Stats package (ver. 3.5.1) in RStudio. Significant miRNA/mRNA ME correlations (at FDR of 5%) were followed up with a more detailed series of analyses, in which individual mRNA hub and miRNA expression was correlated via Pearson’s correlations using the miRLAB package in R (ver. 1.14.3).

Results

AD case/control differentially expressed genes

When we assess each covariate’s contribution to the overall gene expression variance, we see that the impact of a given covariate on gene expression is highly variable across individual genes. Thus, our approach to use a stepwise regression to select for the most influential set of covariates that contribute to the highest proportion of the variance and hence minimize model overfitting. The covariates showing the highest mean contribution to the variance were also the same factors that were most frequently incorporated into our regression models (S2 Table).

The bidirectional stepwise regression revealed 3,536 and 6,401 differentially expressed genes (DEG) in PFC and NAc, respectively, at the nominal p ≤0.05, of which 1,279 DEG were shared between the two regions. Among these, 603 and 494 genes were downregulated and upregulated, respectively, and 182 genes were expressed in opposite directions between the two regions. Within the DEG in NAc, nine genes (ADH1B, ADH1C, H2AFZ, EIF4E, FTO, DRD2, SLC39A8, and VRK2) were implicated in the largest and most recent AD GWAS [33]. At FDR of 5%, we identified 1,841 DEG from the NAc and 70 from the PFC. The miRNA regression analysis identified 430 and 170 nominally significant miRNAs in the NAc and PFC, respectively, with 168 miRNAs differentially expressed in NAc at FDR of 5% with no miRNA reaching FDR significance in PFC. To maintain an identical analytical pipeline for both brain regions and optimize the selection for the most influential confounding factors, we co-jointly analyzed the PFC expression data generated in this study with our previously published NAc expression data [18]. We observed a highly significant overlap between the differentially expressed genes identified in NAc from both studies (Fisher’s exact test, p = 1E-10). For detailed information about regression coefficients, the regression models used for each transcript, partitioning of variance, and frequency of covariates incorporated into the analysis, see S2 Table.

Gene network module clustering

In NAc, at a Bonferroni adjusted p ≤0.05, we identified 6 modules significantly correlated with AD case status (Fig 1A). Among these, NAc was the only negatively correlated module, whereas NAc, NAc, NAc, NAc, and NAc were all positively correlated with AD cases relative to controls (Fig 1B). In PFC, we identified 3 modules significantly correlated to AD at Bonferroni adjusted p≤0.05 (Fig 1C). Of these, the PFC module was negatively correlated, while PFC and PFC were positively correlated with AD cases (Fig 1D). To assess the validity of these network modules, we performed a bootstrap resampling that showed consistent module clustering when compared to the original gene networks (S2 File).

Fig 1

WGCNA clustering and module-trait relationships.

NAc cluster dendrogram and module assignment with dissimilarity based on topological overlap. The 6,401 selected transcripts were clustered into 23 distinct modules. NAc module-trait relationship heatmap correlating (Pearson’s) module MEs with AD diagnosis and covariates. Uncorrected p-values are given in parenthesis below each correlation coefficient. 6 AD associated significant modules (NAc, NAc, NAc NAc, NAc, and NAc,) were identified after Bonferroni correcting p-values (* = p≤0.05). PFC cluster dendrogram and module assignment. The 3,536 selected transcripts were clustered into 17 different co-expressed modules. PFC module-trait relationship heatmap created as previously described. We identified 3 AD associated modules (PFC, PFC, and PFC) after Bonferroni correcting p-values (* = p≤0.05).

NAc and PFC network preservation

We performed a network preservation analysis to determine how well co-expressed networks from the PFC are conserved in NAc and vice versa. We focused primarily on the Z and Median Rank network preservation statistics because Z estimates network overlap by also taking into consideration network connectivity. Median Rank being invariant to module size, provides a more accurate estimate of network preservation since larger networks tend to be more conserved due to their size alone. We observed that NAc and NAc showed little to no network preservation (Zsummary <2), NAc, NAc, PFC, and PFC showed moderate levels of network preservation (2< Zsummary <10), and NAc, NAc and PFC showed high levels of network preservation (Zsummary >10) (Fig 2A and 2B). For detailed information about the individual density and connectivity statistics that were used to create the composite network preservation statistics, see the S3 Table.

Fig 2

Network preservation and gene-set enrichment.

NAc Z-summary statistic calculated as an aggregate of network preservation statistics (Preservation level: high = Z>10; moderate = 2

NAc Z-summary statistic calculated as an aggregate of network preservation statistics (Preservation level: high = Z>10; moderate = 2AD associated modules. PFC Z-summary statistic and corresponding GO biological processes term (-log10(FDR) transformed). Venn-diagram of the shared transcripts from highly preserved NAc modules (NAc and NAc) and their corresponding significant PFC modules (PFC and PFC). Brain cell type gene-set enrichment from the NAc and PFC (-log10(FDR) transformed). Colors correspond with their respective modules (NAc NAc, NAc, PFC, PFC, and PFC) with single gene sets enriched in modules.

Identifying biological processes and cell-type enrichment of the AD significant modules

To gain perspective on the biological underpinnings of the significant gene networks from NAc and PFC, we performed a gene-set enrichment analysis, GO biological processes annotation (ShineyGO ver.61) and neuronal cell type enrichment for the two regions. As one of our aims was to identify unique and shared gene networks associated with AD in NAc and PFC, we focused our analyzes on NAc modules that were highly (i.e., NAc and NAc) and poorly (i.e., NAc, and NAc) preserved in PFC. NAc and NAc are primarily associated with the immune response process (FDR ≤0.05) believed to be a consequence of neurotoxicity caused by chronic alcohol abuse (Fig 2A). These modules are enriched among microglia and astrocyte cell types (FDR ≤0.05), which is expected based on the functional properties of the glial cells (Fig 2D). The poorly preserved NAc modules showed enrichment within gene-sets associated with cilia-based cell projection and cell morphogenesis (FDR≤0.05) (Fig 2A).

Corollary, we performed gene-set enrichment analysis on PFC modules, which were highly and poorly preserved in NAc. (Fig 2B). Similar to the NAc and NAc modules, the highly preserved PFC module was associated with immune response processes (FDR ≤0.05) and significant microglial cell type enrichment (FDR ≤0.05) (Fig 2D). PFC and NAc, were moderately preserved with each other (Fig 2C) with PFC showing astrocyte cell type enrichment (Fig 2D). Interestingly, a class of genes in one family of immune response proteins, metallothioneins (MTs), contained in both the PFC and NAc modules, were differentially expressed in both brain regions between cases and controls (Fig 3). For the complete gene set enrichment analyses for all modules associated with AD in NAc and PFC see S4 Table. Since hubs are considered the most important genes for preserving the network’s integrity, when these analyses were further limited only to the hub genes, not surprisingly, we captured the same GO terms and biological process that we observed from the entire module gene lists (S4 Table).

Fig 3

Metallothionein gene expression.

Relative expression of 8 metallothionein cluster genes (MT1E, MT1F, MT1G, MT1H, MT1HL1, MT1X, MT2A, and MT3) comparing AD case to controls for both the NAc and PFC. P-values presented for each transcript are based on our bidirectional stepwise regression.

Hub genes of potential biological significance

To identify candidate hub genes of potential biological significance, we focused on the relationship between intramodular connectivity (i.e., module membership (MM)) and gene significance (GS) to AD case status (S1 File). Of the 459 genes from the 3 significant PFC modules and the 6 significant modules in NAc, we identified 99 and 433 unique hub genes with MM ≥0.80, respectively. We focus on the hub genes due to their biological relevance to AD and predicted role as drivers of expression for the entire module [34]. For full PFC and NAc transcript information regarding individual MM and GS see S5 Table.

Detection of miRNA gene network modules in NAc and PFC

In NAc and PFC, we identified miRNA modules with varying levels of significant correlation to AD case status. The NAc miRNA data revealed 430 nominally significant loci, which clustered in 5 modules ranging from 18 (NAcmi) to 259 (NAcmi) loci in size, of which, at Bonferroni adjusted p ≤0.05, three miRNA modules remained significantly correlated to AD (NAcmi, NAcmibrown, and NAcmi). Of these, NAcmi and NAcmi were negatively correlated, whereas NAcmi was positively correlated within AD (Fig 4A). The 170 miRNA transcripts from the PFC clustered into 6 modules ranging in size from 9 (PFCmi) to 55 miRNA transcripts (PFCmi), of which PFCmi and PFCmi, remain significant at Bonferroni adjusted p ≤0.05; both miRNA modules were negatively correlated with AD (Fig 4D). For further details on the miRNA WGCNA from both brain regions, including individual transcript MM and GS values, please refer to S6 Table.

Fig 4

MiRNA WGCNA and mRNA:miRNA interaction.

NAc miRNA cluster dendrogram and module assignment with module-trait relationship heatmap, both as previously described in Fig 1. Bonferroni adjusted significant (p≤0.05) NAc mRNA/miRNA module ME correlations (Pearson’s). Alcohol and control groups are separated by color to emphasize sample clustering. Significant (FDR≤0.05) correlation (Pearson’s) between mir-449a and selected mRNA transcripts from the low network preserved NAc module. PFC miRNA cluster dendrogram module assignment along with module-trait relationship heatmap. Bonferroni adjusted significant (p≤0.05) NAc mRNA:miRNA module ME correlations (Pearson’s).

MiRNA networks show unique patterns of regulation

In an attempt to identify a higher order system, network levels of interactions, existing between the AD significant mRNA and miRNA modules we correlated their respective module MEs. From the NAc, we identified 2 significant positive mRNA/miRNA ME correlations and 4 negative ME correlations at Bonferroni adjusted p ≤0.05 (Fig 4B). To better understand the biological function of miRNA/mRNA interacting networks at specific loci, we honed on the interaction between individual miRNA/gene pairs. After correlating individual mRNA hubs and miRNA, we identified 1,801 significant mRNA/miRNA interactions (FDR ≤0.10) spanning 318 genes and 68 miRNA loci (S7 Table). Interestingly, we observed 97% (35/36) of the purple mRNA module hub genes to be negatively correlated with either mir-449a or mir-449b from NAcmi (Fig 4C). In PFC, we identified one positive mRNA/miRNA ME correlation and 3 negative correlations at Bonferroni adjusted p≤0.05 (Fig 4E). Individual mRNA/miRNA interaction analysis from the PFC revealed 6 mRNA/miRNA interactions (FDR of ≤0.10) spanning 6 genes and one miRNA transcript, mir-485-5p. For a full list of individual mRNA/miRNA interactions, see S7 Table.

Brain region specific eQTL regulation of differential gene expression

In NAc, we detected a total of 36 mRNA eQTLs spanning 17 unique genes and 9 miRNA eQTLs covering 4 different miRNA (FDR ≤0.10). Of the 17 hubs with significant eQTLs, 7 are from NAc (VRK1, INPP4A, HMP19, DKK3, PCDH8, RNF34, and RASGRP1), 4 from NAc (FCGR3A, CTSS, AASS, and RNASE4), 3 from NAc (DNALI1, CCDC81, and SPAG6), 2 from NAc (HIVEP1 and GNAS), and one from NAc (VAMP5). Within the PFC we identified 34 eQTLs spanning 16 unique genes and 18 miRNA covering 7 different miRNA transcripts (FDR ≤0.10). Of these, 11 genes are from PFC (SERPINH1, CDKN1A, PNP, EMP1, FKBP5, IL4R, TNFRSF10B, RTEL1/TNFRSF6B, SERPINA1, MAFF, and SERPINA2) and 5 from PFC (GAD2, ACTL6B, KCNF1, SEZ6L, and EFNB3). Among our significant eQTLs, we highlight two examples: FCGR3A:rs12087446 (NAc p = 3.24E-07; PFC p = 0.002) from the highly conserved NAc module and DNALI1:rs12119598 (NAc p = 1.94E-09; PFC p = 0.150) from the poorly conserved NAc module. The brain region specific eQTL impact on the expression of these two genes suggests that different genetic mechanisms are likely at play in NAc and PFC that may further shed light on the different behavioral measures encoded by the two brain regions (Fig 5). For the full list of cis-eQTL, please refer to S8 Table. To highlight the potential clinical importance of our findings and provide functional support for previous genetic studies, we also tested for enrichment of our clinically relevant eQTLs (i.e., testing only SNPs that showed a significant (SNP x AD) interaction term) and previously published GWAS of addiction phenotypes. While the overlap did not reach formal significance, likely due to the smaller GWAS sample size, we nevertheless observed suggestive enrichment, i.e., GSCAN drinks per week p = 0.195; GSCAN smoking initiation p = 0.251; GSCAN smoking cessation p = 0.147; and COGA+Irish p = 0.299. Finally, we attempted to replicate all eQTLs in our study, irrespective of their potential disease relevance, in the GTEx database using the Fisher’s exact test. Interestingly, we observed a significant overlap between our eQTLs detected in the PFC (n = 2,368, 6.6% of eQTLs tested, p = 0.003), but not in the NAc (n = 5,436, 3.4% of eQTL tested, p = 1).

Fig 5

Cis-eQTL analysis.

Cis-eQTL boxplot directly comparing AD case/control designation with the FCGR3A:rs12087446 eQTL from the high network preservation NAc/PFC module and the DNALI1:rs12119598 eQTL from the low preservation NAC module, the relative expression is presented on the y-axis and SNP/genotype on the x-axis. Alternative boxplot visualization of the same cis-eQTL directly comparing differences between brain regions.

Discussion

AUD continues to be a growing public health concern with a complex and poorly understood etiology as recreational alcohol use becomes habitual and problematic. The broad goal of this study is to identify the neurobiological processes associated with chronic alcohol use via analyzing brain region-specific gene networks from the NAc and PFC. To understand the human behavior leading to addiction, it is important to investigate how chronic alcohol use impacts expression changes in the evolutionarily newer cortical areas, in contrast to the older, more evolutionarily conserved subcortical brain regions [35]. Here, we attempt to understand the neurobiological underpinnings of alcohol specific reward conditioning in the NAc and disruption of executive function within PFC [6] through identifying gene networks and biological processes associated with AD that are conserved or unique to each brain region. Additionally, we assessed the relationship between the miRNA and mRNA networks significantly correlated to AD based on the miRNA functions to induce mRNA degradation and/or translational inhibition. Finally, we tested the impact of genetic variants on gene expression in a disease dependent manner via AD-mediated eQTL analysis.

Our network analyses are consistent with previously published reports by others and us, showing the upregulation of immune response mechanisms among AD cases as a byproduct of alcohol’s neurotoxic effects [36]. The immune-related modules show significant enrichment for both astrocyte and microglial cell types, which has been validated by previous alcohol studies and the known immune functions of astrocyte and microglia in the brain [37, 38]. More importantly, we observed generalized up-regulation of immune response mechanisms within both the PFC and NAc, suggesting that the neurotoxic response to chronic alcohol use is ubiquitous across cortical and subcortical brain regions. Interestingly, in both brain regions, we further identified differentially expressed genes in the metallothionein cluster (MT1HL1, MT1H, MT1X, MT1E, MT1G, MT1F, MT2A, and MT3). The metallothionein cluster is primarily responsible for maintaining the cellular homeostasis of zinc and copper while also regulating oxidative stress [39]. Zinc is an essential catalytic cofactor for alcohol metabolism via alcohol dehydrogenase [40]. Free or “chelated” zinc ions (Zn2+) are seen in abundance in the brain, specifically at ionotropic glutamate receptors such as the NMDA receptor family. The interaction between Zn2+ and NMDAR activity has shown to be an important contributor to synaptic plasticity through regulating postsynaptic density assembly [41]. It is well understood that chronic alcohol abuse leads to varying degrees of organ-wide zinc deficiency [42]; however, the neurobiological consequences of how zinc deficiency in the brain contributes to AD neuropathology is poorly understood. We believe this interaction between chronic alcohol abuse, metallothionein expression, zinc deficiency, and synaptic plasticity is an important avenue for future research that should be explored.

In addition to identifying dysregulated immune response mechanisms, we validate recent studies showing differential expression among signaling and neurodevelopmental processes within AD cases [13, 15, 16, 18]. However, these processes are less conserved between cortical and subcortical regions, likely due to the different neuronal composition and functional properties of the PFC and NAc [43]. Interestingly, two NAc modules that primarily associate with cilium assembly (NAc) and cellular localization/morphogenesis (NAc) show limited network preservation within the PFC. There has been increasing evidence suggesting primary cilia aid in facilitating extrasynaptic signaling during adult neurogenesis [44, 45], an important aspect of addiction related extracellular membrane plasticity [46]. For example, GRP88, a g-protein coupled receptor and primary cilia enriched gene [47], was linked to increased alcohol seeking behaviors in knock out (KO) mice models [48], further reinforcing the importance of primary cilia in AUD etiopathology. The cilium assembly genes enriched in NAc, were shown to be associated with axonemal dynein assembly (DNAAF1, DNAI2, and DNALI1). A recent gene expression study in adolescent rat hippocampus identified increased expression of two dynein associated genes (dnai1 and dnah5) [49]. One explanation for increased expression of primary cilia associated genes in the NAc relative to the PFC is related to potential discrepancies in adult neurogenesis between subcortical vs cortical brain regions. It is well understood that most adult neuronal stem cells originate in the ventricular–subventricular zones (V-SVZ) and migrate to adjacent cortical and subcortical brain regions as neuroblasts to promote neurogenesis [50]. A recent study showed increased adult neurogenesis of medium spiny neurons within the NAc and that the migration and incorporation of new neurons was experience-based [51]. We believe that the increased expression of genes that encode for the cilia assembly complex may reflective of experience mediated neurogenesis of medium spiny neurons in NAc, except being driven by chronic alcohol consumption instead of pain. These new neurons formed in response to alcohol use may play an important role in the reward response deficits we often associate with addiction and AUD [4].

Other interesting findings arise from our mRNA/miRNA interactions, e.g., when correlating the MEs from mRNA and miRNA modules, we see distinct patterns between cases and controls within both brain regions. Based on the known function of miRNAs in regulating the expression of target mRNAs [52] we can infer these significant miRNA networks may serve as a driving contributor for differential network expression between AD cases and controls. Specifically, 97% (35/36) of the hub genes from the NAc module were significantly negatively correlated with either mir-449a or mir-449b. Mir-449a/b have primarily been studied in the context of spermatogenesis and cellular proliferation in cancer [53-55]. Based on the mRNA-miRNA correlations, our study suggests that mir-449a/b cluster has additional functions related to cellular proliferation in the brain. Among the genes correlated with mir-449a in the NAc, ELAVL4, DPYSL3, and KCNJ6 have shown significant associations with AD in other expression, and genetic association studies [19, 56, 57], as well as being implicated in other substance use disorders [58-61].

In an attempt to understand the causal nature of the gene networks associated with AD, we integrated genetic information via eQTL analysis. We were able to detect a significant number of mRNA and miRNA cis-eQTLs from both brain regions. We selected highly significant eQTLs (FCGR3A (Fc fragment of IgG receptor IIIa):rs12087446 and DNALI1 (dynein axonemal light intermediate chain 1):rs12119598) based on FCGR3A and DNALI1’s role as network hubs to highlight the interaction between AD case status and eQTL while also demonstrating brain region-specific eQTL variation. FCGR3A is one of the low-affinity Fc receptor genes important for NK cell-mediated antibody-dependent cytotoxicity [62] and a hub gene from our highly conserved NAc and PFC modules. The consistent effect of rs12087446 on FCGR3A expression between both brain regions suggests the genetic impact on immune response processes might also be ubiquitous across the brain of chronic alcohol users. Differential FCGR3A expression was recently shown to be associated with both alcohol preference and binge-like behaviors in the ventral tegmental area of rats [63]. In contrast, DNALI1, a hub gene in the cilium assembly enriched NAc module, is under the genetic control of specific eQTL only in NAc but not in PFC, suggesting that changes to cilia organization due to alcohol abuse might be under different genetic control between the two brain regions. We observed suggestive evidence for enrichment between our eQTLs and previously published GWAS of alcohol or other addiction phenotypes, such as smoking. We believe this is primarily due to three factors: 1) low statistical power within our sample to detect genetic signals that would otherwise appear in large-scale GWAS studies, 2) our selective study design focusing only on potentially clinically relevant eQTLs, and 3) the presence of variants with a lower MAF in the GWAS potentially not detectable in our dataset. We further successfully replicated our eQTLs in the GTEx database for PFC, but not NAc. One possible explanation is that the increased number of DEG in the NAc relative to PFC with the fact GTEx does not include AD diagnosed brains in their analyses [64] effectively limits our ability to replicate GTEx eQTLs based on significant and potential subtle non-significant expression changes among AD cases.

Conclusion

The strength of this study lies in our ability to compare and contrast expression changes between subjects with AD and controls within two different brain regions. We successfully identified gene networks and biological processes from both brain regions that were validated by previous AD studies as well implicated a novel biological process (cilia assembly) and gene family (metallothionein cluster) as potentially important for the development of AD. Our mRNA/miRNA interaction analysis pinpointed mir-449a/b cluster as an important regulator of differentially expressed genes between AD cases and controls. Finally, via our eQTL analysis, we provided evidence that mRNA and miRNA expression differences between AD cases and controls might be under brain region specific genetic control. While our sample size could be perceived as a limitation, we mitigated this by utilizing WGCNA to aggregate differentially expressed genes into biologically relevant modules with single expression values, effectively increasing our power to detect significant AD associations within a multivariate framework. Additionally, to increase the power of our study, considering the more prevalent and heavier drinking patterns in men, we assessed the molecular processes of alcohol drinking in male subjects only. While we recognize the importance of comparing the molecular pathology of drinking between the two sexes, we would like to highlight observations from genetic epidemiological studies showing male and female subjects to have a similar genetic predisposition to alcohol abuse [65]. We further recognize that a number of our significant AD associated modules in PFC were also nominally correlated to neuropathology (p≤0.05). This is not entirely unexpected, given the known neuropathological impact of chronic alcohol abuse [66]. Finally, while we understand that the lower RINs from the PFC can be seen as confounding factor, studies have suggested that reliable data can still be obtained from postmortem brain tissue even with suboptimal RNA quality [67-69]. However, our careful analytical design to adjust for the impact of RIN on gene expression maintains the robustness of our results even in the presence of lower RINs.

Overall, the broader impact of our findings is the understanding that chronic alcohol consumption can reinforce addiction behaviors through dysregulating different biological process across various brain regions. This information could potentially lead to more focused therapies for AUD by targeting important brain regions specific neurobiological pathways involved in the development of alcohol addiction. While our results point to certain biological processes that differentiate between the PFC and NAc, these findings require replication in an independent postmortem brain samples spanning other cortical and subcortical brain regions. Additional support for the postmortem brain findings presented here can also be obtained by studying ethanol activity in animal models or neuronal cell cultures. Increased research within the methodological framework of our study can help validate our findings and identify biological processes and genes that play the most significant role in the development of AUD.

Expanded methodology descriptions (bi-directional stepwise regression, WGCNA, addiction GWAS enrichment, and GTEx eQTL replication) and supplemental figure displaying MM by gene significance for AD associated mRNA and miRNA modules.

(DOCX)

Click here for additional data file.

100 permutation robust WGCNA dendrogram clustering.

(PDF)

Click here for additional data file.

Sample demographics.

(XLSX)

Click here for additional data file.

Stepwise regression coefficients, models, covariate frequency counts and variance partitioning (VP).

3.1) mRNA, 3.2) miRNA.

(ZIP)

Click here for additional data file.

Network preservation Z-summary table and supplemental network preservation statistics.

(XLSX)

Click here for additional data file.

Full GO biological processes annotation for each AD associated module from the NAc and PFC.

6.1) NAc, 6.2) PFC.

(ZIP)

Click here for additional data file.

mRNA WGCNA module membership (MM) and gene significance(GS) values with AD associated modules isolated in separate tabs.

7.1) NAc, 7.2) PFC.

(ZIP)

Click here for additional data file.

miRNA WGCNA module membership (MM) and gene significance(GS) values with AD associated modules isolated in separate tabs.

8.1) NAc,8.2) PFC.

(ZIP)

Click here for additional data file.

Top mRNA/miRNA correlations (NAc = top 2000; PFC = Top 500).

(XLSX)

Click here for additional data file.

cis-eQTL analysis results.

(XLSX)

Click here for additional data file.

28 Oct 2020

PONE-D-20-30169

Network Preservation Reveals Shared and Unique Biological Processes Associated with Chronic Alcohol Abuse in NAc and PFC.

PLOS ONE

Dear Dr. Vornholt,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands.

You will see the detailed comments of the two Reviewers below. However, I am particularly concerned with the potential RNA degradation in the PFC samples, as mentioned by Reviewer 2 . Differential degradation could impact brain regional analysis. I hope you can adequately address this concern. If you can, I also suggest to add a more generalized take-home message to Discussion section.

If you feel that are able to adequately address these main concerns, I invite you to submit a revised version of the manuscript that addresses all points raised during the review process.

If you decide to do so, please submit your revised manuscript by Dec 11 2020 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

We look forward to receiving your revised manuscript.

Kind regards,

Andrey E Ryabinin, Ph.D.

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: Rigorous and interesting analysis of mRNA and miRNA expression associated with alcohol dependence in post-mortem male human brain.

Methods.

Very detailed and well described methods.

Might be interesting to the reader to comment on which covariates had a larger effect on gene expression.

S1 Table. Probably need better legend. What do the numbers mean for Hemisphere, Neuropathology, Hepathology (#should be hepatology?), Toxicology, Smoking.

Results.

Figure 3. Significance levels for DE not noted in figure, legend, or text.

Hub genes of potential biological significance. This section is less interesting. Are they enriched for relevant biological processes or GWAS hits? Are they sig. DE between cases and controls? Other than their high degree of connectivity within modules enriched for AD, what is there biological significance? Or maybe integrate into another results section?

Figure 5. Units of expression not clearly described in the legend.

Discussion.

Further discussion of the interesting finding of cilium assembly in the NAc but not in the PFC may be warranted. Adult neurogenesis in the NAc? Neurogenesis of medium spiny neurons in the nucleus accumbens continues into adulthood and is enhanced by pathological pain. García-González et al., Mol. Psych. 2020. PMID: 32612250

Reviewer #2: The report by Vornholt et al. is a descriptive study that utilizing genome-wide mRNA and miRNA expression from postmortem brain of AD cases compared to matched controls. Gene co-expression networks, mRNA/miRNA interactions, and eQTL analyses were used to identify expression patterns that are uniquely and differently dysregulated in the PFC and NAc. The analytical pipelines are well described and mostly robust.

The Methods indicates 42 cases and 42 controls. Then after matching, 72 total cases. Table S1 shows demos for 70 cases (1 case was a statistical outlier). It is not clear why the initial 12 cases were excluded from the original group (84).

It is not clear from the beginning of the Methods section that the NAc data were reported previously by the same group. Utilization of the previously published data was clarified in the first section of the Results but needs to be outlined clearly in the Methods. Also, for clarity, any annotation updates and data handling differences for the dataset should be listed since the data were published in 2015.

There is significant concern about the within subject RNA quality differences for each of the brain regions. A major goal of the reported work was to identify overlapping and unique expression differences in PFC and NAc. There are clear differences in RNA quality between the 2 brain regions. This is likely degradation due to the sample processing at different times. The authors argue reliable data can be obtained from RINs<=4; however, the impact of RNA degradation on downstream expression analysis can be significant, especially when trying to compare multiple brain regions from the same individual. The cross-region analysis is likely compromised.

The major components of the study (descriptive statistics, WGCNA analysis, network preservation, mRNA/miRNA interactions) are potentially interesting; however, there is not a well-integrated overview or thematic model presented from consolidated findings. It is especially hard to understand exactly what the reader will gather from complex mRNA/miRNA interactions.

Functional validation of key findings would significantly strengthen the study. At a minimum, testable hypotheses should be discussed based on the integrated analyses.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

2 Nov 2020

Editors Comment: You will see the detailed comments of the two Reviewers below. However, I am particularly concerned with the potential RNA degradation in the PFC samples, as mentioned by Reviewer 2 . Differential degradation could impact brain regional analysis. I hope you can adequately address this concern. If you can, I also suggest to add a more generalized take-home message to Discussion section.

We appreciate Editor’s feedback and hope to have addressed the comment in respect to PFC RINs below in our response to Reviewer #2. Additionally, we have added a more generalized take-home message at the end of the discussion (pg. 21; line 534-538)

Reviewer #1: Rigorous and interesting analysis of mRNA and miRNA expression associated with alcohol dependence in post-mortem male human brain.

Methods.

Very detailed and well described methods.

Might be interesting to the reader to comment on which covariates had a larger effect on gene expression.

We appreciate the reviewer’s comments and have elaborated in greater details the importance of each included covariate on gene expression across the two brain regions (pg. 6; line 169-172). Briefly, in the stepwise regression model, the importance of a given covariate to affect gene expression was based on the frequency by which each covariate is included in the final model. This additional information is also available in Table S3.

S1 Table. Probably need better legend. What do the numbers mean for Hemisphere, Neuropathology, Hepathology (#should be hepatology?), Toxicology, Smoking.

We apologize for the confusion and have clarified the data labeling in Table S1.

Results.

Figure 3. Significance levels for DE not noted in figure, legend, or text.

We have added the DE significance levels for each MT gene within the figure (pg. 12; line 324-325).

Hub genes of potential biological significance. This section is less interesting. Are they enriched for relevant biological processes or GWAS hits? Are they sig. DE between cases and controls? Other than their high degree of connectivity within modules enriched for AD, what is there biological significance? Or maybe integrate into another results section?

We appreciate the reviewer’s comments and have included the GO biological processes enrichment for just the hub genes in Table S6. Generally, in network analyses, due to their high connectivity, hub genes play an important role in maintaining the integrity of the biological network. Furthermore, the hub genes are the strongest contributors to summarized ME expression and are also considered the most important genes for future downstream, more functional, studies (pg. 13; line 334-335). Therefore, not surprisingly, we see consistent enrichment of the same GO terms between the two tests, i.e., when only the hub genes are tested versus the entire module (pg.12; 319-322). For these reasons, we also focused our eQTL analysis on the hub genes from AD significant modules only.

Figure 5. Units of expression not clearly described in the legend.

We have changed the legend to clarify the x and y-axis for each plot in Figure 5 (pg. 16; line 406).

Discussion.

Further discussion of the interesting finding of cilium assembly in the NAc but not in the PFC may be warranted. Adult neurogenesis in the NAc? Neurogenesis of medium spiny neurons in the nucleus accumbens continues into adulthood and is enhanced by pathological pain. García-González et al., Mol. Psych. 2020. PMID: 32612250.

We welcome the suggestion and have included additional discussion on the implications of cilium assembly within the context of neurogenesis in the NAc (pg. 18; line 458-469)

Reviewer #2: The report by Vornholt et al. is a descriptive study that utilizing genome-wide mRNA and miRNA expression from postmortem brain of AD cases compared to matched controls. Gene co-expression networks, mRNA/miRNA interactions, and eQTL analyses were used to identify expression patterns that are uniquely and differently dysregulated in the PFC and NAc. The analytical pipelines are well described and mostly robust.

The Methods indicates 42 cases and 42 controls. Then after matching, 72 total cases. Table S1 shows demos for 70 cases (1 case was a statistical outlier). It is not clear why the initial 12 cases were excluded from the original group (84).

We appreciate the reviewer’s comments and have clarified this confusion. Essentially, we began with 82 samples (41 cases and 41 controls). These samples were then matched for available covariates presented on pg. 5 lines 140-141, leaving a total of 36 matched cases-controls for downstream analysis. Once the expression data were normalized, one final subject (a case) was removed due to presenting itself as a statistical outlier.

It is not clear from the beginning of the Methods section that the NAc data were reported previously by the same group. Utilization of the previously published data was clarified in the first section of the Results but needs to be outlined clearly in the Methods. Also, for clarity, any annotation updates and data handling differences for the dataset should be listed since the data were published in 2015.

We have clarified in the methods that the raw data generated in NAc were collected and normalized as part of the previous study. However, considering that both brain regions were assayed on the same microarray platform, and to maintain the consistency of our normalization and analytical pipelines, we opted to re-analyze the expression data from the NAc together with PFC using the same analytical approach (see pg. 6, 166-168).

There is significant concern about the within subject RNA quality differences for each of the brain regions. A major goal of the reported work was to identify overlapping and unique expression differences in PFC and NAc. There are clear differences in RNA quality between the 2 brain regions. This is likely degradation due to the sample processing at different times. The authors argue reliable data can be obtained from RINs<=4; however, the impact of RNA degradation on downstream expression analysis can be significant, especially when trying to compare multiple brain regions from the same individual. The cross-region analysis is likely compromised.

We agree with the reviewer’s comment. RNA quality and how it may affect subsequent gene expression analyses has been among the most tackled problems in the postmortem brain research area. While generally there has been a positive relationship between RIN and detection of gene expression differences, studies have also shown that reliable gene expression data can also be obtained from RNA with lower RINs [1, 2]. (pg. 20-21; line 529-533). It is also important to note that differences in RINs between cases and controls has a much more detrimental role than low RINs itself [3]. While the RINs in the PFC regions were lower, we did not detect significant differences in RIN scores between cases and controls. Finally, while we acknowledge that the lower PFC RINs may affect our detection power due to our careful analytical design to adjust for the impact of RIN on gene expression, the robustness of our identified differentially expressed genes will remain.

The major components of the study (descriptive statistics, WGCNA analysis, network preservation, mRNA/miRNA interactions) are potentially interesting; however, there is not a well-integrated overview or thematic model presented from consolidated findings. It is especially hard to understand exactly what the reader will gather from complex mRNA/miRNA interactions.

We have included additional details (pg. 4; line 117-126 and pg.16; line 418-423) to address this comment in order to provide a more meaningful interpretation of the data with respect to alcohol dependence.

Functional validation of key findings would significantly strengthen the study. At a minimum, testable hypotheses should be discussed based on the integrated analyses.

We appreciate the reviewer’s comment. The main focus of our study is to compare the gene expression profiles between two critical for alcohol addiction cortical and subcortical brain regions. We have, however, discussed in greater detail future in vivo experiments using animal models and/or neuronal cell cultures, as well as expanding on the methodological framework of this study by increasing postmortem brain sample sizes as well as include additional comparisons from other brain regions. (pg.21, line 538-544)

References:

[1] Weis S, Llenos IC, Dulay JR, Elashoff M, Martínez-Murillo F, Miller CL. Quality control for microarray analysis of human brain samples: The impact of postmortem factors, RNA characteristics, and histopathology. Journal of Neuroscience Methods 2007; 165: 198–209.

[2] Stan AD, Ghose S, Gao X-M, Roberts RC, Lewis-Amezcua K, Hatanpaa KJ, et al. Human postmortem tissue: What quality markers matter? Brain Research 2006; 1123: 1–11.

[3] Gallego Romero I, Pai AA, Tung J, Gilad Y. RNA-seq: impact of RNA degradation on transcript quantification. BMC Biology 2014; 12: 42.

Submitted filename: Response To Reviewers.docx

Click here for additional data file.

19 Nov 2020

PONE-D-20-30169R1

Network preservation reveals shared and unique biological processes associated with chronic alcohol abuse in NAc and PFC.

PLOS ONE

Dear Dr. Vornholt,

Thank you for resubmitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Specifically, there is still a concern regarding the low RINs. As one of the reviewers suggests, there could be  more quantitative arguments on this issue. In addition, I could not find the database used for sample access.

Please submit your revised manuscript by Jan 03 2021 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

We look forward to receiving your revised manuscript.

Kind regards,

Andrey E Ryabinin, Ph.D.

Academic Editor

PLOS ONE

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: (No Response)

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Partly

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: I Don't Know

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

Reviewer #2: There is lingering concern about the impact of the low RIN samples on the analysis. The authors argue that RIN values are consistently low within PFC (thus impacting cases and controls similarly); however, the concern from the initial review concerned a major goal of the work--to identify overlapping and unique expression differences in PFC and NAc. At a minimum, the proportion of variance explained by the known covariates such as age, gender, RIN and PMI should be reported using violin plots or other graphics. Linear mixed models can then be used to quantify the contribution of multiple sources of variation and identify any covariates that should be corrected in the final analysis.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

22 Nov 2020

Response to Reviewers (Revision 2)

Editor Comments:Specifically, there is still a concern regarding the low RINs. As one of the reviewers suggests, there could be more quantitative arguments on this issue. In addition, I could not find the database used for sample access.

We appreciate the editor’s comments and have uploaded our data to GEO on NCBI. The mRNA and miRNA expression data for the NAc was uploaded as part of a previous study (GSE62699). The PFC mRNA and miRNA data is uploaded and currently being processed, but for the sake of expedience we are resubmitting before obtaining an accession number. We will provide the required accession number once it is released to us.

Reviewer #2: There is lingering concern about the impact of the low RIN samples on the analysis. The authors argue that RIN values are consistently low within PFC (thus impacting cases and controls similarly); however, the concern from the initial review concerned a major goal of the work--to identify overlapping and unique expression differences in PFC and NAc. At a minimum, the proportion of variance explained by the known covariates such as age, gender, RIN and PMI should be reported using violin plots or other graphics. Linear mixed models can then be used to quantify the contribution of multiple sources of variation and identify any covariates that should be corrected in the final analysis.

We appreciate the reviewer’s comment. As suggested by the reviewer, we have provided a graphical representation (Tables S3) of the contribution of each covariate impact on gene expression. We have also provided further justification on using the bidirectional stepwise regression to correct for the covariates’ effect (Pg. 6 Line 172 & 173) on gene expression (Pg. 9 Line 236-241). Furthermore, while the use of linear mixed models is well developed to adjust for hidden covariates effect in gene associations studies, their applicability in gene expression studies might be less evident (for example see Yao at al. (2019)). We also note that by using stepwise regression, the covariates’ effect is assessed individually on each gene, allowing for the inclusion of only the most relevant covariates, thus minimizing model overfitting. Finally, the validity of our model to correctly capture and adjust for the RIN’s effect on gene expression is reflected by the increased number of times RIN is incorporated as a covariate in the regression models assessing gene expression in PFC vs NAc (See Table S3). Therefore, we are confident in our ability to control for the RIN or any other covariates’ effect in our analyses.

Yao C, Joehanes R, Johnson AD, Huan T, Liu C, Freedman JE, et al. Dynamic Role of trans Regulation of Gene Expression in Relation to Complex Traits. Am J Hum Genet 2017; 100: 571–580.

Submitted filename: Response to Reviewers_Rev2.docx

Click here for additional data file.

30 Nov 2020

Network preservation reveals shared and unique biological processes associated with chronic alcohol abuse in NAc and PFC.

PONE-D-20-30169R2

Dear Dr. Vornholt,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Andrey E Ryabinin, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

7 Dec 2020

PONE-D-20-30169R2

Network preservation reveals shared and unique biological processes associated with chronic alcohol abuse in NAc and PFC.

Dear Dr. Vornholt:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Andrey E Ryabinin

Academic Editor

PLOS ONE