Identifying a novel biological mechanism for alcohol addiction associated with circRNA networks acting as potential miRNA sponges.

Our lab and others have shown that chronic alcohol use leads to gene and miRNA expression changes across the mesocorticolimbic (MCL) system. Circular RNAs (circRNAs) are noncoding RNAs that form closed-loop structures and are reported to alter gene expression through miRNA sequestration, thus providing a potentially novel neurobiological mechanism for the development of alcohol dependence (AD). Genome-wide expression of circRNA was assessed in the nucleus accumbens (NAc) from 32 AD-matched cases/controls. Significant circRNAs (unadj. p ≤ 0.05) were identified via regression and clustered in circRNA networks via weighted gene co-expression network analysis (WGCNA). CircRNA interactions with previously generated mRNA and miRNA were detected via correlation and bioinformatic analyses. Significant circRNAs (N = 542) clustered in nine significant AD modules (FWER p ≤ 0.05), within which we identified 137 circRNA hubs. We detected 23 significant circRNA-miRNA-mRNA interactions (FDR ≤ 0.10). Among these, circRNA-406742 and miR-1200 significantly interact with the highest number of mRNA, including genes associated with neuronal functioning and alcohol addiction (HRAS, PRKCB, HOMER1, and PCLO). Finally, we integrate genotypic information that revealed 96 significant circRNA expression quantitative trait loci (eQTLs) (unadj. p ≤ 0.002) that showed significant enrichment within recent alcohol use disorder (AUD) and smoking genome-wide association study (GWAS). To our knowledge, this is the first study to examine the role of circRNA in the neuropathology of AD. We show that circRNAs impact mRNA expression by interacting with miRNA in the NAc of AD subjects. More importantly, we provide indirect evidence for the clinical importance of circRNA in the development of AUD by detecting a significant enrichment of our circRNA eQTLs among GWAS of substance abuse.

INTRODUCTION

Alcohol is among the most readily available and commonly abused recreational drugs worldwide with substantial socioeconomic and public health implications. The shift from recreational alcohol use to problematic drinking resulting in alcohol use disorder (AUD) is dependent upon genetic and environmental factors. AUD is moderately heritable (~49%) ; however, the genetic mechanisms underlying this heritability are poorly understood. Although the alcohol dehydrogenase cluster on Chromosome 4 has been among the most consistently replicated genetic loci associated with AUD, molecular studies from the mesocorticolimbic (MCL) system of human postmortem brains and animal models have implicated additional AUD risk genes involved in neurosignaling, synaptogenesis, and immune response. , The limited overlap between molecular and genetic studies has hindered our understanding of the link between AUD‐associated genetic loci and gene expression changes in the brain. Broadly, the human transcriptome can be divided into coding and noncoding, with the noncoding transcriptome (represented by a large set of noncoding RNA [ncRNA] species characterized by their minimal or complete lack of protein‐coding abilities and gene regulatory functions , ) being a largely unexplored domain of the human genome with a potentially substantial impact on the neuropathology of AUD. Among these, a particular class of ncRNA, termed circular RNA (circRNA), have been implicated in the development of alcoholic hepatitis in mouse models. ,

CircRNAs are abundantly and dynamically expressed throughout the mammalian central nervous system (CNS). , They primarily arise from pre‐mRNA splicing events in which the 5′ and 3′ ends of introns or alternatively spliced exons are covalently linked to form closed‐loop structures. Although several hypotheses have been proposed to explain the mechanisms by which circRNAs regulate gene expression, a commonly accepted one, based on experimental observations, is the miRNA sponge hypothesis. MiRNAs regulate gene expression mainly through binding to the 3′ untranslated regions (UTRs) of their target genes, leading to translational repression and mRNA degradation. CircRNAs serve as competing endogenous RNAs (ceRNA) for miRNA by competing with miRNA response elements (MREs) in the 3′ UTRs of mRNA. This leads to miRNA sequestration by circRNA and decreased miRNA–target interactions, effectively increasing gene expression as a result. More importantly, circRNA's ability to function as miRNA sponges has shown to be emerging biomarkers for both neurodegenerative and neuropsychiatric disorders. , ,

With their varied spatiotemporal expression in the brain, circRNAs were implicated in the etiology of neurodegenerative and neuropsychiatric disorders. , , To test whether these recent observations also extend to alcohol dependence (AD), we assessed the genome‐wide expression of circRNA, miRNA, and mRNA in the nucleus accumbens (NAc) from subjects with AD followed by weighted gene co‐expression network (WGCNA) and bioinformatic and statistical analyses (Figure 1). We assess the functional impact of circRNA expression in NAc based on (1) its role as part of the reward pathway central to the neuropathology of addiction, (2) availability of miRNA and mRNA data generated in NAc on the same subjects, and (3) the known characteristic of the NAc as a hub for adult neurogenesis, an important neurobiological process believed to be partially regulated circRNA interactions. Finally, we applied an expression quantitative trait loci (eQTL) analysis to identify genetic elements affecting circRNA expression and the ability to interact with miRNA and mRNA. With this study, our main goals were to identify the potential regulatory mechanisms by which circRNAs affect the expression of risk AUD genes and provide a methodological framework for exploring circRNA, miRNA, and mRNA interactions in future postmortem brain studies. This is the first study to examine the effect of circRNA on mRNA expression via miRNA sponge interactions in NAc from chronic alcohol abusers to the best of our knowledge.

FIGURE 1

Framework for circRNAs as miRNA sponges and study design flowchart. (A) CircRNAs are primarily formed through back splicing of unspliced transcripts in which introns or a combination of exons and introns have their 3′ and 5′ ends covalently bonded to form closed‐end loops. (B) Under normal circumstance, miRNA will bind to 3′ UTR of mature mRNAs, leading to mRNA degradation or translational repression; however, in the presence of circRNA with complementary sequences, miRNAs are sequestered away from their target mRNAs, leading to increased gene expression. (C) Flowchart depicting the steps and analyses used to determine significant circRNA–miRNA–mRNA interactions in this study

METHODS AND MATERIALS

Tissue processing and RNA extraction

Postmortem NAc from 42 AD cases and 42 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 on Alcohol Abuse and Alcoholism, and the New South Wales Department of Health. As part of a previous study, several criteria were used to exclude samples with (1) agonal state, (2) substantial brain damage, (3) history of infectious disease, and (4) postmortem interval (PMI) > 48 h (Table S1). Samples were further matched for RNA integrity number (RIN) (mean = 6.9, ±0.84), sex (all male), ethnicity (100% Caucasian), brain pH, and postmortem interval to minimize covariates' effect on expression, resulting in 18 matched case–control pairs (n = 36). Total RNA from flash‐frozen NAc was extracted and purified via mirVANA‐PARIS kit (Life Technologies, Carlsbad, CA) following the manufacturer's protocol. RNA integrity (RIN) and concentrations were assessed via Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA) and Quant‐iT Broad Range RNA Assay kit (Life Technologies), respectively.

Microarrays and expression normalization

Genome‐wide circRNA, miRNA, and mRNA expression was assessed on three different platforms: (1) Arraystar Human Circular RNA Array spanning 13,617 circRNA probes, (2) Affymetrix GeneChip miRNA 3.0 Array spanning 1733 mature miRNAs, and (3) Affymetrix GeneChip Human Genome U133A 2.0 array containing 22,214 probe sets spanning ~18,400 unique mRNAs. Raw expression data from each assay were background corrected, log2 transformed, and quantile normalized via Partek Genomics Suite v6.23 (PGS; Partek Inc., St. Louis, MO) and the limma package (Version 4.0) in R. To exclude outliers that could impact downstream analyses, three samples were removed from the circRNA normalized dataset, leaving 17 cases and 16 controls (n = 33), and one sample was removed from both the miRNA and mRNA normalized datasets, resulting in 17 cases and 18 controls (n = 35). Here, we validate only the circRNA array by assessing the expression of three randomly selected circRNA from the entire microarray at the Arraystar facilities via quantitative PCR (qPCR) because the mRNA and miRNA expression arrays were validated previously. The assessed genes showed a high mean correlation (Kendall tau r = 0.87 [SD ± 0.021]) between the two platforms (Figure S1).

Identifying differential transcript expression

We assessed the relationship between transcript expression and AD status in RStudio (ver. 1.2.1335). Differential circRNA expression was assessed via robust linear regression in the MASS package (v.7.351.5) with smoking and RIN included as covariates in the model as these were shown to have a greater impact on circRNA expression, , compared with other demographic and postmortem covariates. Differentially expressed miRNA and mRNA were previously identified via a bidirectional stepwise regression adjusting for demographic and postmortem covariates in the Stats package (v.3.6.1).

WGCNA

To better capture the system‐wide network interactions between the circRNA and increase power to detect circRNAs associated with AD, we partition the nominally significant (p ≤ 0.05) circRNAs into co‐expression networks using WGCNA in RStudio (v.1.69). Our criteria to include nominally significant genes were based on retaining genes with (1) smaller effect sizes, albeit true positive signals; (2) exclude genes not likely associated with AD (i.e., reduce the noise); and (3) provide a sufficient number of genes for the network analysis. A detailed description of WGCNA is outlined in literature , , and in Methods S1. The module eigengenes (MEs), a single aggregate expression value for each gene module, were correlated to AD case status and available demographic/biological covariates. To validate the gene networks associated with AD in WGCNA, we performed a bootstrap resampling of 100 iterations with replacement (Figure S2).

CircRNA hub gene prioritization

CircRNA hubs representing potential drivers of expression for entire modules were identified from the absolute value of Pearson's correlation coefficient between MEs and individual gene expression. This value of intramodular connectedness, denoted as module membership (MM), was used to define circRNA hubs as transcripts significantly correlated with AD (p ≤ 0.05) and a MM ≥ 0.70 within the significant AD modules.

Correlations analysis between circRNA, miRNA, and mRNA

To test the miRNA sponge hypothesis framework, we first explored correlations between the differentially expressed mRNA, miRNA, and circRNA hubs. We used only subjects with complete data across all three expression platforms (i.e., 17 AD cases and 15 controls). The circRNA–miRNA–mRNA correlations were based on Pearson's product moment generated in the miRLAB package (ver. 1.14.3) in RStudio. All significant circRNA–miRNA, negative miRNA–mRNA, and positive circRNA–mRNA correlations, respectively, were extracted at a false discovery rate (FDR) of 10% and retained for follow up analyses.

Computational prediction of circRNA–miRNA interactions

Next, the correlation analyses were supplemented with bioinformatic predictions to further increase the reliability of our findings by identifying circRNA, miRNA, and mRNA pairs with converging evidence for interaction from both the correlation and in silico analyses. Thus, the circRNA–miRNA correlations were complemented with computational predictions using STarMir in the Sfold application suite (http://sfold.wadsworth.org/cgi-bin/index.pl). STarMir calculates probability scores for binding predictions of shared seed sequences between circRNA and miRNA based on logistic regression models developed from crosslinking immunoprecipitation (CLIP) studies. Based on STarMir's recommendations, the logarithm of odds for predicted binding (logit probability score) ≥ 0.50 was considered significant.

Prediction of miRNA–mRNA target interactions

Similarly, the miRNA–mRNA correlations complemented miRNA target predictions from the multiMiR package (v.1.6.0) in Rstudio. MultiMiR is a curated database of miRNA–mRNA target predictions that integrates both computational prediction algorithms (DIANA‐microT‐CDS, ElMMo, MicroCosm, miRanda, miRDB, PicTar, PITA, and TargetScan) and experimentally validated miRNA–target interactions (miRecords, miRTarBase, and TarBase).

Moderation analysis

To test whether miRNA expression moderates the relationship between circRNA and target mRNAs we utilized the Stats package in RStudio. Within the linear regression model that included mRNA expression as the independent variable, moderation was modelled as the interaction (circRNA × miRNA) term effect on gene expression adjusting for smoking and RIN's confounding effects. Significance was based on an FDR ≤ 0.10 threshold. For a more detailed methodological explanation, see Methods S2.

Gene‐set enrichment analyses

To determine the biological function of genes that participate in our identified circRNA–miRNA–mRNA interactions, we performed a GO biological processes gene‐set enrichment via ShinyGO (v.0.61) at each stage in our analyses. ShinyGO utilizes a hypergeometric distribution to determine significant enrichment at an FDR ≤ 0.10.

CircRNA eQTL analysis and enrichment in genome‐wide association study of substance abuse

To explore if known AUD‐associated genetic variants impact the expression of differentially expressed circRNA hubs, an eQTL analysis was performed, followed by a genome‐wide association study (GWAS) enrichment analysis. The postmortem sample was genotyped as part of a larger GWAS study. Monomorphic single nucleotide polymorphisms (SNPs) and those with excessive missingness (>20%) were filtered out. Only local, cis‐eQTLs within 500 kb of each circRNA hub's start/stop position were mapped, and variants in linkage disequilibrium (LD) (R 2 ≥ 0.7) were subsequently pruned via Plink v1.9. We utilized the MatrixEQTL package (ver. 2.3) in RStudio to perform two cis‐eQTL analyses using a linear regression framework after adjusting for relevant covariates. Firstly, we utilized the “modelLINEAR” argument to identify significant eQTLs irrespective of AD case status. Secondly, we incorporated the “modelLINEAR_CROSS,” which introduces a (SNP × AD) interaction term to assess whether any significant eQTL impacts circRNA expression in a disease dependent manner. The overlap (i.e., enrichment) was tested between our eQTLs and recent GWAS of substance abuse, including alcohol and smoking, using two mutually complementing tests (Cauchy combination [CC] and Simes , ) adjusting for multiple testing and LD (R 2 ≥ 0.50). For more details, see Methods S3.

RESULTS

CircRNAs are organized in networks associated with AD

At the nominal p ≤ 0.05, our gene expression analysis revealed 542 differentially expressed circRNAs between AD cases and controls (Figure 2A), with none of them achieving significance at FDR ≤ 0.10. The lack of significant circRNA detected at the FDR of 10% can be attributed to low power (even in the presence of rigorous case–control matching) due to the relatively low expression of individual circRNA transcripts. Among others, low circRNA expression could lead to increased within‐group variance due to purely stochastic events. Thus, to detect potential circRNA drivers of the AD phenotype, we employed WGCNA, which aggregates individual transcripts into co‐expressed modules. In addition to increased statistical power, WGCNA also highlights network interactions at a system level. Except for one module (M3), the nominally significant circRNAs clustered in 10 modules, nine of which significantly correlated with AD (Bonferroni adj. p ≤ 0.05) (Figure 2B). Of the significant modules, M1, M2, M4, and M5 were positively correlated, whereas M6–M10 were negatively correlated with AD status. From these modules, we identified 137 hub genes, which were selected for downstream statistical and bioinformatic analyses against mRNA (n = 3575) and miRNA (n = 264) significantly associated with AD at FDR ≤ 0.10 as a part of a previous study on the same subjects (Figure 2A). Please see (1) Table S2 for regression models and coefficients, (2) Figures S2 and S3 for circRNA regression QQplots and top differentially expressed transcripts, and (3) Table S3 for WGCNA results.

FIGURE 2

Differentially expressed transcripts and circRNA WGCNA results. (A) Volcano plots describing the relationship between regression estimates and −log10(p) for each transcript level in our analysis (circRNA, miRNA, and mRNA). Dashed lines correspond with the significance threshold of p ≤ 0.05 and FDR ≤ 0.10. (B) WGCNA module clustering dendrogram from our nominally AD significant (p ≤ 0.05) circRNA transcripts. (C) Heat plot comparing the correlation (Pearson's) of our identified circRNA module MEs to AD diagnosis and all other available covariates. In respect to AD diagnosis, the top value represents the correlation coefficient, and the bottom value represents uncorrected p‐values. For covariates: *p ≤ 0.05 and **p ≤ 0.005

CircRNA–miRNA–mRNA show complex interaction patterns that are associated with AD

We tested the circRNA ability to interact with miRNA and thus indirectly affect the miRNA target's expression in a disease‐dependent manner. Assuming circRNAs act as miRNA sponges to impact mRNA expression, we posit that the most relevant downstream biological interactions will be represented by negative miRNA–mRNA and positive circRNA–mRNA correlations. We also include negative circRNA–miRNA correlations to capture the predicted inverse relationship between circRNA–miRNA expression as previously described. , Thus, we first performed three independent correlation analyses (circRNA–miRNA, miRNA–mRNA, and circRNA–mRNA) followed by tests to identify the intersection between the significant correlations corrected at FDR of 10%. In the circRNA–miRNA (circRNA n = 137; miRNA n = 264) analysis, we identified 48 significant negative circRNA–miRNA correlations. The miRNA–mRNA (miRNA n = 264; mRNA n = 3,575) analysis revealed 46,501 significant negative correlations. Finally, the circRNA–mRNA (circRNA n = 137; mRNA n = 3575) analysis revealed 2221 significant positive correlations. From the intersection of these analyses, we identified a total of 2480 overlapping correlations, which were then used in all subsequent follow‐up analyses. For a full list of correlation coefficients, see Table S4.

Binding predictions supplement intersecting circRNA–miRNA–mRNA correlations

To reinforce and complement our correlation analyses, the 2480 overlapping circRNA–miRNA–mRNA correlations were further screened computationally to identify predicted circRNA–miRNA and miRNA–mRNA interacting pairs. Based on STarMir's algorithm, no circRNA–miRNA binding predictions with a score greater than our significance threshold (logit probability ≥ 0.50) were detected when circRNA–miRNA correlations were considered in isolation. However, by expanding the circRNA–miRNA binding predictions to include circRNA–miRNA pairs correlated with the same mRNA, we identified 365 circRNA–miRNA–mRNA trios with intersecting negative miRNA–mRNA correlations, positive circRNA–mRNA correlations, and predicted circRNA–miRNA binding. We further narrow down the 365 interactions via selecting the best miRNA–mRNA target predictions to identify the most reliable 47 circRNA, miRNA, and mRNA participating in a three‐way interaction. For a complete list of binding predictions, see Table S5.

Moderation analysis reveals circRNA impact mRNA expression via potential miRNA sequestration

The impact of miRNA sequestration on mRNA expression from these 47 circRNA, miRNA, and mRNA was formally tested in a linear regression model adjusting for AD status and controlling for covariate effects (RIN and smoking history). The miRNA sequestration by circRNA was assessed by introducing a (circRNA × miRNA) interaction term in the model. At FDR ≤ 0.10, we identified 23 interactions that show a significant moderation effect on mRNA expression (Table 1). Interestingly, among these 23 interactions, circRNA‐406702–miR‐1200 stood out by affecting the expression of the largest set of mRNA (n = 17), the four most significant of which (HRAS, PRKCB, HOMER1, and PCLO) that have been linked to chronic alcohol abuse in previous studies are highlighted in Figure 3. For full moderation regression coefficients, see Table S6.

TABLE 1

Top circRNA–miRNA–mRNA interactions

circRNA–miRNA–mRNA interactions			miRNA–mRNA Cor		circRNA–mRNA Cor		Circ × mi interaction		circRNA–miRNA binding		miRNA–mRNA target Predic
circRNA	miRNA	mRNA	Coef	FDR	Coef	FDR	Estimate	FDR	Logit prob	Seed	Predicted	Experimental
ASCRP3011575	hsa‐miR‐1200	ACTR2	−0.4251	0.0909	0.5044	0.0998	0.8127	0.0898	0.6955	offset‐6mer		X
ASCRP3011575	hsa‐miR‐1200	ASTN1	−0.5425	0.0260	0.5057	0.0987	2.0620	0.0070	0.6955	offset‐6mer		X
ASCRP3011575	hsa‐miR‐1200	ATP2B2	−0.5330	0.0293	0.6079	0.0440	2.5879	0.0077	0.6955	offset‐6mer		X
ASCRP3001917	hsa‐miR‐4310	CELF1	−0.4440	0.0773	0.5324	0.0812	0.6517	0.0898	0.7029	7mer‐m8	X	X
ASCRP3011575	hsa‐miR‐1200	E2F3	−0.5165	0.0361	0.5530	0.0708	0.9043	0.0913	0.6955	offset‐6mer	X
ASCRP3010153	hsa‐miR‐3187‐3p	GPD2	−0.5295	0.0307	0.5681	0.0622	−0.7690	0.0380	0.6087	7mer‐A1	X
ASCRP3011575	hsa‐miR‐1200	HOMER1	−0.5717	0.0169	0.6755	0.0237	1.7868	0.0077	0.6955	offset‐6mer	X
ASCRP3011575	hsa‐miR‐1200	HRAS	−0.4523	0.0715	0.6014	0.0470	1.8261	0.0077	0.6955	offset‐6mer	X
ASCRP3011575	hsa‐miR‐1200	IMP4	−0.5091	0.0394	0.5483	0.0739	1.5305	0.0008	0.6955	offset‐6mer		X
ASCRP3011575	hsa‐miR‐1200	IPCEF1	−0.5550	0.0218	0.5356	0.0801	2.0946	0.0086	0.6955	offset‐6mer	X
ASCRP3011575	hsa‐miR‐1200	LDB2	−0.4212	0.0939	0.6082	0.0440	3.5155	0.0086	0.6955	offset‐6mer	X
ASCRP3005132	hsa‐miR‐665	MLEC	−0.5167	0.0361	0.5433	0.0762	−0.6227	0.0134	0.6796	6mer	X
ASCRP3011575	hsa‐miR‐1200	NDST3	−0.5663	0.0185	0.5400	0.0779	2.7353	0.0096	0.6955	offset‐6mer	X
ASCRP3012325	hsa‐miR‐361–5p	NEK7	−0.5016	0.0430	0.6356	0.0342	2.2585	0.0898	0.7444	offset‐6mer	X	X
ASCRP3013378	hsa‐miR‐571	NR3C1	−0.4644	0.0639	0.5547	0.0697	−1.5579	0.0550	0.7044	offset‐6mer		X
ASCRP3011575	hsa‐miR‐1200	OSBPL8	−0.5532	0.0223	0.5434	0.0761	1.9161	0.0077	0.6955	offset‐6mer	X
ASCRP3011575	hsa‐miR‐1200	PCLO	−0.5141	0.0371	0.5592	0.0677	2.5603	0.0070	0.6955	offset‐6mer	X
ASCRP3011575	hsa‐miR‐1200	PRKCB	−0.5650	0.0188	0.5873	0.0535	5.3517	0.0009	0.6955	offset‐6mer	X
ASCRP3011575	hsa‐miR‐1200	RAB11FIP2	−0.5455	0.0250	0.5758	0.0583	0.9061	0.0104	0.6955	offset‐6mer	X
ASCRP3011575	hsa‐miR‐1200	RANBP2	−0.4977	0.0450	0.5858	0.0543	0.8803	0.0368	0.6955	offset‐6mer		X
ASCRP3011575	hsa‐miR‐1200	RFC2	−0.4256	0.0906	0.5162	0.0915	0.8356	0.0380	0.6955	offset‐6mer		X
ASCRP3011575	hsa‐miR‐1200	SSX2IP	−0.7589	0.0001	0.5269	0.0847	2.0238	0.0009	0.6955	offset‐6mer	X

Note: Significant circRNA–miRNA–mRNA interactions that survive all of our bioinformatic and statistical tests (i.e., negative miRNA–mRNA correlation, positive circRNA–mRNA correlation, circRNA–miRNA predicted binding, miRNA–mRNA target prediction, and moderation regression).

FIGURE 3

CircRNA‐406702–miR‐1200 interacting transsynaptic signaling‐associated genes. (A) Boxplot showing relative microarray expression differences between AD cases and controls for miR‐1200. (B) Diagram of predicted binding loci between circRNA‐406702 and miR‐1200. (C) Boxplot showing relative expression differences between AD cases and controls for circRNA‐406702. (D) Correlation plots displaying the significant negative relationship between miR‐1200 and interacting transsynaptic signaling‐associated genes (HRAS r 2 = −0.45; PRKCB r 2 = −0.57; HOMER1 r 2 = −0.57; PCLO r 2 = −0.51). (E) Correlation plot displaying significant positive relationship between circRNA‐406702 and select genes (HRAS r 2 = 0.61; PRKCB r 2 = 0.59; HOMER1 r 2 = 0.68; PCLO r 2 = 0.56). (F) Boxplots for differential mRNA expression between AD cases and controls and diagram of miRNA predicted binding to the 3′ UTR of target genes

CircRNAs interact with genes associated with neuronal function

At each stage of our analyses, we consistently identified significant enrichment (FDR ≤ 0.10) of genes involved in cellular localization, synaptic transmission, neural development, and response to organic stimuli gene sets (Figure 4). The 22 circRNA–miRNA–mRNA interactions also revealed significant enrichment (FDR ≤ 0.10) for GO biological processes associated with regulation of DNA metabolism, anatomical structure homeostasis, regulation of biosynthesis, dendritic spine organization, and anterograde transsynaptic signaling. For full gene‐set enrichment, see Table S7.

FIGURE 4

Identification of significant circRNA–miRNA–mRNA interactions and GO biological processes enrichment. (A) Breakdown of the number of significant circRNA–miRNA–mRNA interactions and unique genes at each step in our analysis (overlapping positive circRNA–mRNA and negative miRNA–mRNA correlations [Step 1], circRNA–miRNA binding predictions [Step 2], miRNA–mRNA binding predictions [Step 3], and moderation regression [Step 4]) ending with circRNA‐406702–miR‐1200 interacting mRNA. (B) GO biological processes enrichment for each set of unique genes at each step in our analysis. The genes or the number of genes from our list is presented within each histogram of the associated gen ‐set

Genetic variants potentially impact circRNA expression

Our linear eQTL analysis revealed three significant circRNA eQTLs at an FDR ≤ 0.10 (Figure 5A). When we repeated the eQTL analysis taking into consideration the interaction (AD × genotype) term, we detect seven additional significant eQTLs (FDR ≤ 0.10) independent from our main analysis, which were associated with one circRNA (circRNA‐080252). After expanding our original linear eQTL analysis to incorporate results at a more relaxed significance threshold (unadj. p ≤ 0.002), we identify an additional 96 eQTLs used in the downstream enrichment analysis. Among these, we identified multiple circRNAs that participated in significant circRNA–miRNA–mRNA interactions at various stages in our multistep analyses (Figure 5B,C). For full eQTL results, see Table S8.

FIGURE 5

Significant circRNA cis‐eQTLs. (A) eQTLs that survive FDR ≤ 0.10 significance threshold. (B) eQTLs from circRNA–miRNA–mRNA trios with negatively correlated miRNA–mRNA, positively correlated miRNA–mRNA, predicted circRNA–miRNA binding, and miRNA–mRNA predicted interactions. (C) eQTLs for circRNA that participate in circRNA–miRNA–mRNA interactions that survive all our bioinformatics and statistical tests

CircRNA‐associated SNPs are enriched within AUD and smoking GWAS

We employed the CC and Simes , tests to detect eQTLs (n = 96) and SNPs in LD with them (r 2 ≥ 0.50; n = 1558) that were enriched among the significant (p ≤ 5E‐4) loci from recent GWAS of AUD and smoking (GWAS & Sequencing Consortium of Alcohol and Nicotine Use (GSCAN) and Psychiatric Genetics Consortium AUD GWAS (PGC‐AUD) ). Adjusting for multiple testing and background enrichment, we observed significant enrichment for our eQTLs in GSCAN cigarettes per day (CC p = 1.41E‐05; Simes p = 2.74E‐05), GSCAN smoking initiation (CC p = 0.025), and PGC‐AUD European ancestry (CC p = 0.034).

DISCUSSION

In recent years, circRNAs have become increasingly relevant for the study of neuropsychiatric and neurodegenerative disorders as potential regulators of the brain's complex and unique transcriptome. Here, we examined the intersection between the pathological characteristics of AD (aberrant neurogenesis and accelerated neurodegeneration , ) and circRNA's predicted role in regulating neurobiological processes within the context of circRNA–miRNA–mRNA interactions.

Our study relied upon a series of experimental, statistical, and bioinformatics tests to narrow down well over a billion possible interactions between circRNA, miRNA, and mRNA to highly specific three‐way interactions within the miRNA sponge hypothesis that survive several layers of correction for multiple testing. Among our most significant circRNA–miRNA interacting pairs (i.e., circRNA‐406702–miR‐1200), we observed a unique set of genes negatively correlated with miR‐1200 and positively correlated with circRNA‐406702. Some of these (such as HRAS, PRKCB, HOMER1, PCLO, ASTN1, and ATP2B2) are enriched within gene sets associated with synaptic transmission/development, highlighting their potential importance to the neuropathology of AD. HRAS, a small GTP‐binding protein, interacts with downstream PI3K, AKT, and mTORC1 as part of a neurosignaling pathway (“Go” pathway) believed to be important for promoting neuroadaptations associated with excessive alcohol consumption and relapse. This is supported by studies showing that HRAS expression is increased among mice strains consuming alcohol in high quantities, as well as in the NAc of rats with an extended history of excessive consumption followed by periods of abstinence. However, in contrast to the animal‐based studies, in our sample, we observed decreased HRAS expression in AD subjects. A possible explanation would be that the ligand‐gated ion channels mediating HRAS activity become desensitized due to chronic receptor activation after years of alcohol abuse, which cannot be easily replicated in animal models. , PRKCB (protein kinase C beta), another gene implicated in our study, is an isoform of the protein kinase C (PKC) family. This set of proteins is shown to be essential for the development of AD through their interaction with CREB‐BDNF neurosignaling pathway, which was reported to be associated with synaptic plasticity. , , More importantly, genetic variants nearby PRKCB have been significantly associated with comorbid bipolar disorder, substance use disorder (SUD), and alcohol cue‐elicited brain activation.

Among the other genes interacting with circRNA‐406702–miR‐1200 are HOMER1 and PCLO, which encode for proteins playing an important role at the synapse. HOMER1 encodes for one of the Homer scaffolding proteins (Homer1/2), which link metabotropic glutamate receptors (mGlu1/5) to the postsynaptic density. Both HOMER1 and one of the mGlu receptor, GRM5, have been consistently implicated as potential therapeutic targets for the treatment of AD due to their role in regulating alcohol‐facilitated neuroplasticity. , Additionally, it has been shown that a polymorphism (rs7713917) in the regulatory region of HOMER1 can help predict increased alcohol consumption in adolescents years later. PCLO codes for the Piccolo protein, a scaffolding protein at the active zone of the presynaptic cytomatrix, an area where neurotransmitters are released. Intronic SNPs within the PCLO gene have been one of the most studied genetic variants associated with major depressive disorder. Functional studies have suggested that these polymorphisms may play an active role in emotional memory processing, , and previous research has indicated that deficits in emotional processing is a hallmark of AD. This deficit then may lead to enhanced emotional reactivity to positive and negative stimuli during periods of drinking and periods of withdrawal, effectively reinforcing continued alcohol abuse. , Importantly, PLCO and HOMER1 have both been implicated as differentially expressed in multiple gene expression studies of AD. , , , , Finally, ASTN1 (astrotactin 1) is a gene that codes for a protein receptor important for glial‐guided neuron migration. In the context of AD, a family‐based linkage study has shown that ASTN1 is significantly associated with AD in multiplex families. Overall, the results from our study provide further support for research suggesting circRNA play an important, yet still underexplored, role in neuronal function.

Some of the miRNAs implicated at various steps in our circRNA analysis, although not all of them directly associated with AD, show significant associations with alcoholic liver disease, brain function, and neuropsychiatric disorders. Among the several miRNA identified from our significant circRNA–miRNA interactions, miR‐665 is significantly upregulated in the prefrontal cortex (PFC) of alcoholics, and miR‐361‐5p shows increased expression in the PFC of early‐stage AD mouse models. The maternal expression of another miRNA from our study (miR‐3119) was shown to increase following alcohol consumption during pregnancy. Two other miRNAs (miR‐1200 and miR‐3187‐3p) have been implicated in various neurobiological processes relevant to AD etiology. Of these, miR‐1200 has been predicted to regulate neuronal connexins 36, 45, and 57 in humans, mice, and rats. Connexins (Cx) are essential for gap junction function at electrical synapses, with Cx36 shown to be associated with various rewarding effects of alcohol intoxication in knockout (KO) mice. Another report has suggested that miR‐3187‐3p expression changes modify the neuronal cell response to oxidative stress. Increased oxidative stress is a well‐known consequence of alcohol's neurotoxic effects in the brain with multiple studies from our group and others identifying increased expression of immune and stress response genes in the postmortem brains of chronic alcohol users. , , Finally, miR‐571 has shown to be an important biomarker for alcohol‐related liver disease. We further show that miR‐571 interacts significantly with NR3C1, a highly pleiotropic glucocorticoid receptor necessary for stress response and reported to be significantly associated with AD. ,

In respect to our cis‐eQTL analysis, we identify genetic variants that impact the expression of our circRNA hubs. Though no specific polymorphisms at the genome‐wide significance level (p ≤ 5E‐8) in GWAS of AUD were replicated among our eQTLs, we observed significant enrichment at a lower significance threshold (p ≤ 5e‐4) using two separate genomic enrichment tests using recent GWAS of AUD and smoking. , Possible explanations for this observation are (1) the limited power of our postmortem brain sample, (2) different methods for phenotypic classification between our study and GWAS, and (3) GWAS of AUD that are still underpowered. However, most likely, with increased postmortem brain sample sizes and deep phenotyping of subjects with chronic alcohol abuse, we may begin to see a meaningful overlap between the results from these two methods. Nevertheless, the importance of identifying eQTL enrichment among GWAS signals from our eQTL analysis is threefold: first, help validate the clinical relevance of these large association studies by providing a functional explanation for AUD associated GWAS signals; second, reinforce such identified eQTLs and SNPs in LD as likely candidates for future, more targeted, follow‐up analyses; and third, provide limited genetic context for the causal relationship between gene expression changes and AUD. Our study also highlights a potentially novel neurobiological mechanism of alcohol addiction by demonstrating that alcohol abuse may impact known AD risk genes by altering circRNA expression and circRNA's ability to act as miRNA sponges. Our circRNA eQTL study further suggests that we must be careful when interpreting GWAS signals given that genetic variants impacting the expression of proximal circRNA can alter the expression of distal genes through epistatic interaction between circRNA and miRNA.

Limitations and future direction

Our study does also have a few limitations. First, it is possible that by focusing solely on the circRNA and miRNA interactions, we may have overlooked other molecular mechanisms (i.e., epigenetic factors) that potentially can also affect the functions of risk AD loci. Second, although the use of male subjects only can be perceived as a limitation, this was a deliberate choice in order to increase our statistical power by removing sex‐based variability. Genetic epidemiological studies have shown that male and female subjects have a similar genetic predisposition to alcohol abuse. Also, we acknowledge our limited ability to differentiate between expression changes that represent the neuropathological consequence of chronic alcohol abuse or predictive factors for the development of AUD. However, this does not diminish the scientific and clinical relevance of our findings because identifying genes and gene networks that promote or maintain problematic drinking behaviors can inform novel therapeutic treatments for AUD. We also integrate genetic information via eQTL in order to tease apart the relationship between predictive genetic factors, AUD, and expression. ,

To the best of our knowledge, ours is the first to specifically investigate the effect of circRNA and miRNA interactions on gene expression in NAc from subjects with AD. We are confident that this pilot study will facilitate and promote future studies to corroborate our findings by experimentally validating these results and further exploring them in the context of increased and more diverse experimental settings. Future studies might include exploring the relationship between circRNA–miRNA–mRNA interactions and other alcohol‐related traits such as alcohol consumption. Including alcohol consumption is beyond the scope of this study given that our samples are obtained from individuals diagnosed with chronic and severe AUD and that the genetic correlation between consumption and severe AUD is relatively weak. We do feel alcohol consumption would be an important factor to include when exploring more intermediate/moderate drinking phenotypes in future studies. It is also important to acknowledge the relevant crossover between the identified genes participating in circRNA–miRNA–mRNA interactions and genes implicated in studies from animal models of AUD. Much of the functional neurobiological work has been performed in animal models, whereas large‐scale GWAS have attempted to identify heritable genetic variants associated with AUD and other addictive behaviors in human populations. , However, the translation between human genetic studies and functional animal studies has been limited. This is most likely due to the complex nonlinear relationships between environmental and genetic factors in humans that are difficult to recapitulate in animal models of AD and the increased transcriptome complexity in the human brain relative to animal brains and other tissues. , Thus, follow‐up functional studies in animal models or proxy tissues are necessary to corroborate these findings. Moreover, we believe that our study will be the steppingstone on which future studies will expand our integrative analytical approach to incorporate other brain regions and psychiatric phenotypes.

CONFLICT OF INTEREST

The authors report no biomedical financial interests or potential conflicts of interest.

Data S1. Supplemental Methods:Additional methodological explanations for WGCNA, Moderation regression and GWAS enrichment.

Click here for additional data file.

Figure S1: Validation of circRNA microarray via qPCR. The validity of the Arraystar Human Circular RNA Array was assessed by performing qPCR on 3 randomly selected circRNA. Overall, we identify a high mean correlation (Kendall tau r = 0.87 [SD ± 0.021]) between the two platforms.

Figure S2: Quantile‐Quantile plot of circRNA regression p‐values. –Log10 transformed p‐values from our analysis were compared to the expected distribution of p‐values via QQplot.

Figure S3: Boxplots of the top differentially expressed circRNA. Relative expression is presented on the y‐axis and AD case status on the x‐axis. A) circRNA‐001328, B) circRNA‐100417 C) circRNA‐103510, D) circRNA‐406203, E) circRNA‐100641, F) circRNA‐015279.

Figure S4: Robust WGCNA dendrogram module clustering. To ensure network robustness and minimize the potential effect of outlier samples on network structure, we used the robust ‘bootstrapped’ version of WGCNA (rWGNCA). We performed 100 iterations in which networks were created after randomly subsetting 2/3 of the total sample. The resulting 100 networks were merged into one large, final consensus network with the individual sub‐networks showing reasonably high consistency with the final networks.

Click here for additional data file.

Table S1: Sample demographics.

Click here for additional data file.

Table S2: CircRNA, miRNA, and mRNA regression coefficients, as well as the covariates used for each model based on bidirectional stepwise regression for miRNA and mRNA.

Click here for additional data file.

Table S3: CircRNA WGCNA results which include the module membership (MM) and gene significance (GS) for each transcript.

Click here for additional data file.

Table S4: Direction specific Pearson's correlations coefficients for each RNA pair (negative circRNA:miRNA correlation, negative mRNA:miRNA correlation and positive circRNA:mRNA 850 correlation)

Click here for additional data file.

Table S5: CircRNA:miRNA binding prediction results obtained via STarMir and our miRNA:mRNA target prediction obtained from MultiMir.

Click here for additional data file.

Table S6: Linear regression coefficients for the circRNA x miRNA interaction term from our moderation analysis.

Click here for additional data file.

Table S7: Full GO biological processes gene‐set enrichment from our unique mRNA at each stage in our analysis.

Click here for additional data file.

Table S8: CircRNA eQTL results from both our linear and AD x genotype interaction (linear cross) analysis. Additionally, we present the full list of SNPs used for our GWAS enrichment analysis.

Click here for additional data file.