RNA-Seq Analysis of Genetic and Transcriptome Network Effects of Dual-Trait Selection for Ethanol Preference and Withdrawal Using SOT and NOT Genetic Models.

BACKGROUND: Genetic factors significantly affect alcohol consumption and vulnerability to withdrawal. Furthermore, some genetic models showing predisposition to severe withdrawal are also predisposed to low ethanol (EtOH) consumption and vice versa, even when tested independently in naïve animals.
METHODS: Beginning with a C57BL/6J × DBA/2J F2 intercross founder population, animals were simultaneously selectively bred for both high alcohol consumption and low acute withdrawal (SOT line), or vice versa (NOT line). Using randomly chosen fourth selected generation (S4) mice (N = 18-22/sex/line), RNA-Seq was employed to assess genome-wide gene expression in ventral striatum. The MegaMUGA array was used to detect genome-wide genotypic differences. Differential gene expression and the weighted gene co-expression network analysis were implemented as described elsewhere (Genes Brain Behav 16, 2017, 462).
RESULTS: The new selection of the SOT and NOT lines was similar to that reported previously (Alcohol Clin Exp Res 38, 2014, 2915). One thousand eight hundred and sixteen transcripts were detected as differentially expressed between the lines. For genes more highly expressed in the SOT line, there was enrichment in genes associated with cell adhesion, synapse organization, and postsynaptic membrane. The genes with a cell adhesion annotation included 23 protocadherins, Mpdz and Dlg2. Genes with a postsynaptic membrane annotation included Gabrb3, Gphn, Grid1, Grin2b, Grin2c, and Grm3. The genes more highly expressed in the NOT line were enriched in a network module (red) with annotations associated with mitochondrial function. Several of these genes were module hub nodes, and these included Nedd8, Guk1, Elof1, Ndufa8, and Atp6v1f.
CONCLUSIONS: Marked effects of selection on gene expression were detected. The NOT line was characterized by higher expression of hub nodes associated with mitochondrial function. Genes more highly expressed in the SOT aligned with previous findings, for example, Colville and colleagues (Genes Brain Behav 16, 2017, 462) that both high EtOH preference and consumption are associated with effects on cell adhesion and glutamate synaptic plasticity.

C57BL/6J (B6) and DBA/2J (D2) mouse strains differ markedly in a number of ethanol (EtOH)‐related phenotypes including EtOH preference and acute withdrawal (Chester et al., 2003; Metten et al., 1998). B6 mice have a high preference and low withdrawal sensitivity, while the D2 mice show the opposite profile. Recombinant inbred strains and selected lines derived from B6xD2 founders have repeatedly illustrated a consistent and negative genetic correlation between preference and withdrawal (Metten et al., 1998). Moving to mice derived from more complex cross this negative genetic relationship is diminished. For example, using mice selectively bred from a 4‐way heterogeneous stock (HS) that included the B6 and D2 strains, lines bred for high and low preference showed the expected (negative) acute withdrawal scores. In contrast, animals bred for high and low acute withdrawal did not differ in preference consumption (Hitzemann et al., 2009). Using replicate withdrawal sensitive prone (WSP) and withdrawal sensitive resistant (WSR) mouse lines derived from even more complex HS founders, Crabbe and colleagues (2013) observed that compared to the WSP‐2 line, the WSR‐2 line had a higher EtOH preference; however, when comparing the WSR‐1 and the WSP‐1 lines, preference consumption was higher in the WSP‐1 line. Our overall interpretation of these data is that while the negative correlation between preference and withdrawal may be limited to only some genotypes; however, our understanding of this relationship will inform and direct studies that move to more diverse genotypes. Further, the considerable genetic resources available to study animals derived from B6xD2 crosses (Mulligan et al., 2017) leverage this endeavor.

With these perspectives in mind, Metten and colleagues (2014) initiated a dual selection study in which beginning with a B6xD2 F2 intercross, animals were selectively bred for high preference/low withdrawal and vice versa. The former were denoted as SOTs from the Old English for habitual drunkard; animals from the reverse selection were simply noted as NOTs. After 3 generations of selection, the divergence between the lines was robust and met the expected outcomes. QTL analysis revealed significant loci on chromosomes 2, 4, and 9, in regions previously associated with preference or acute withdrawal. Microarray‐ and network‐based analyses of the ventral striatum suggested that selection had its strongest effect on a network module significantly enriched in kinase‐related annotations.

Iancu and colleagues (2012) were the first to observe that when implementing network‐based gene expression analyses, RNA‐Seq had significant advantages over microarray‐based approaches. The advantages were mainly associated with the more robust variance structure of the RNA‐Seq data; that is, there are practically no floor or ceiling effects in the assessment of gene expression. The variance expansion produces gene coexpression network modules with significantly higher intramodular connectivity. With these data in mind and given a continuing interest in the SOT/NOT phenotypes, a new selection of the SOT/NOT lines was initiated with the primary goal of using RNA‐Seq to detect the effects of selection on the ventral striatum transcriptome. Importantly, the data presented here extend the observations of Metten and colleagues (2014) in several important aspects and perhaps most importantly in showing that selection has marked effects on genes associated with mitochondrial function.

Materials and Methods

Animal Husbandry

All mice used for the short‐term selection (STS) were obtained from the colony at OHSU, an AAALAC Animal Care‐approved facility. A description of STS including the total number of mice behaviorally tested is included in the Supplemental Material. Mice were maintained in polycarbonate/polysulfone cages and housed on Eco‐Fresh bedding (Absorption Corp., Ferndale, WA, USA) with free‐access to tap water and Purina 5001 chow (PMI Nutrition International, Brentwood, MO, USA). Offspring were weaned and group‐housed (2 to 5/cage) with same‐sex littermates. Procedure and colony rooms were kept at a temperature of 21 ± 1°C with lights on from 0600 to 1800 hour. Adult animals (8+ weeks of age) from STS 1‐3 were used for behavioral testing. Naive animals from STS generation 4 (S4) were euthanized by cervical dislocation, and the brain was quickly removed and stored at −80°C. All procedures were approved by the OHSU Institutional Animal Care and Use Committees in accordance with the United States Department of Agriculture and the United States Public Health Service guidelines.

Dual Behavioral Trait Testing

EtOH Preference

Male and female mice were individually housed and habituated for 1 week. EtOH consumption and preference were measured using a 24‐hour access, 2‐bottle free‐choice procedure. Graduated glass cylinders (25‐ml) containing EtOH (Decon Laboratories, Inc., King of Prussia, PA) or tap water were fitted with rubber stoppers and stainless‐steel sipper tubes (Phillips, Crabbe, Metten, Belknap, 1994). In order to acclimate mice to the taste and effects of EtOH, the animals received access to progressively increasing concentrations of EtOH (0, 5, and 10% v/v) every 4 days. The position of the graduated cylinders was reversed every other day to control for side preferences. EtOH intake, averaged over the second and fourth day of 10% EtOH access, was used to calculate EtOH consumption in g/kg. EtOH preference, measured as the EtOH consumed (in ml) divided by total fluid (water + EtOH) intake, was also determined and averaged over the second and fourth day of 10% EtOH consumption.

Acute EtOH Withdrawal

One week after completing the EtOH preference test mice was weighed and scored twice to measure baseline handling‐induced convulsion (HIC) severity (Kozell et al., 2008). Immediately after the second baseline HIC measurement, mice were injected with 4 g/kg EtOH intraperitoneally (20% v/v in physiological saline). Mice were HIC‐scored beginning 2 hours after injection and hourly through 12 hours. After the EtOH injection, HIC scores were suppressed for several hours, then usually exceeded baseline at 4 hours and increasing for several hours, peaking ~7 to 8 postethanol. Withdrawal scores from 2 to 12 hours were calculated. Withdrawal scores for individual animals were indexed as the total area under the curve (AUC) regressed on the average baseline (RESIDUAL WD).

In order to select breeding pairs for subsequent selection of SOT and NOT generations, a rank ordering of preference (greater % of total liquid consumed as EtOH solution) was used followed by residual withdrawal severity scores. The males and females (20 pairs) with the highest preference/lowest withdrawal values were paired, avoiding brother–sister mating, to create the SOT line; similarly, those mice (20 pairs) with the lowest preference/highest withdrawal scores were paired to create NOT line. For each generation (S1 to S3) at least 200 offspring from SOT and NOT lines combined were tested as described above. Active selection ended at S3 with S4 alcohol‐naive pups was used for RNA sequencing. The overall analysis strategy is illustrated in Fig. 1.

Fig. 1

Schematic summarizing methodology. Details of the methods used for experimental animals and data analysis are summarized and shown in order of process.

The MegaMUGA genotyping array (GeneSeek [Neogen], Lincoln NE) was used to determine the effects of selection on the SOT and NOT genotypes (N = 84 animals). Animals genotyped were also those used for RNA sequencing (see below). For this analysis, genetic distances were calculated for each pair of samples determined by the number of identical alleles at each marker. This was summed over all markers. Between‐group and within‐group distances were calculated. Finally, an analysis of molecular variance approach was used (Excoffier et al., 1992).

RNA Sample Preparation

Naïve adult S4 mice were euthanized by cervical dislocation. Whole brains were removed and frozen immediately using dry ice that was followed up by storage at −80°C. Under RNAse‐free conditions, the brains were thawed on ice and hand dissected. The midhypothalamus was used as a caudal marker for the ventral striatum, and tissue was isolated from a 1‐mm coronal slice of brain tissue. The striatum is somewhat teardrop‐shaped, so beginning at the medial–ventral aspect of the striatum, the dissection was 1 mm dorsal, followed by a cut on the lateral boundary of the striatum with a final cut along the lateral–ventral striatal boundary. The isolated tissue was immediately placed into Trizol (1 ml, Invitrogen; Carlsbad, CA). RNA was extracted by following the manufacturer’s directions for Trizol extraction and then stored in RNAse‐free water at −80°C until transfer to the OHSU Massively Parallel Sequencing Shared Resource (MPSSR) Core.

RNA Sequencing

The MPSSR verified intact RNA quality using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA), and the 84 highest quality samples were selected for library construction. Ovation RNA‐Seq v2 kit (NuGen Technologies, San Carlos, CA) was used to create amplified cDNA (total of ~4 ng of total input RNA). Libraries were constructed using the TruSeq DNA library prep kit (Illumina, Inc., San Diego, CA). Samples were multiplexed 6 per lane on an Illumina HiSeq 2500 sequencer yielding 30 to 50 million single‐end reads (100 bp) per sample. Base call data were assembled into Fastq files using the Illumina CASAVA suite. FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used for quality assessment on the raw sequence reads. Reads were then aligned to the mouse genome (mm10) using STAR 2.5.3a (Spliced Transcripts Alignment to a Reference [Dobin et al., 2013]) allowing for a maximum of 3 mismatches per 100 bp read. Nonuniquely mapped reads were excluded from the alignment process. Before alignment, the total reads for each sample ranged from 37 to 56 million. From those, ~90% were uniquely mapped (i.e., from 33 to 49 million reads) to the genome. Read counts per gene were then obtained using feature‐Counts 1.6.0 (SubRead Package). Differential expression (DE) analysis between SOT and NOT lines was performed with the R (3.5.1) edgeR 3.24.3 package utilizing a generalized linear model and controlling for gender‐specific expression differences. The read density threshold for inclusion in the data analyses for genes was 1 count‐per‐million (CPM) average reads across samples. Three samples were removed from subsequent data analysis. Two SOT samples were misidentified in the genotyping array, and 1 SOT sample was an outlier (>2 SD from mean) in an Interarray correlation plot. Expression data are available at NCBI Gene Expression Omnibus (GEO accession number http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE140567).

Consensus Network Analysis

Coexpression Network Construction and Effects of Selection

Weighted gene coexpression network analysis (WGCNA; Iancu et al., 2012; Langfelder and Horvath, 2008) was used to construct a coexpression network (Figs 1 and 2). A consensus module approach was used followed by an assessment of selection effects on network structure (Ando et al., 2015; Gill et al., 2010; Iancu et al., 2013; Ideker and Krogan, 2012). Details are found in Hitzemann and colleagues (2019). A flowchart showing the overall data analysis strategy is shown in Fig. 1. Genes within the SOT and NOT networks were classified as Hub genes if the relative intramodular was greater than 0.8 (see Table S4).

Fig. 2

Schematic summarizing module and network details following WGCNA. Only results for DE, DV, and DW are shown where NOT > SOT. For SOT > NOT network, there were only 2 genes showing intramodular connectivity (Ralgapa2 and Zhx1), neither were hub nodes. For DW, there were 9 genes with higher relative connectivity in the SOT line; however, no common annotation was detected. Genes that are bolded and blue are those included in the GeneMANIA analysis of DE‐ and DW‐enriched modules. GeneMANIA coexpression network annotation indicated annotation found in results for both DE and DW.

DE, DV, and DW Analyses

Figure 2 details the network and enrichment annotation for DE, DV, and DW. The DE was determined as described in Colville and colleagues (2017). In brief, we used edgeR, with “tagwise” as optional dispersion; a false discovery rate (FDR) of <0.05 was used to set a threshold for significance. The DE genes are described as more highly expressed in 1 line; importantly, higher expression is not synonymous with upregulation. Apparent higher expression (SOT > NOT) could be the result of lower gene expression in the NOT line or vice versa for NOT > SOT. A “var.test” procedure in the statistical package R was used to determine gene DV (differential variation). Differential wiring (DW) was determined using the difference in correlation strength as evaluated using the Levene test function from the R (3.0.2) package “car” with the search for Pearson correlations, which was restricted between individual genes differing by ≥0.5 before a power transformation. This procedure was used to identify the number of changed edges for each gene. To determine whether some genes had a disproportionately higher number of changing edges, the binomial test was used with the following parameters: The average incidence of changing edges (the rate of the binomial test) was computed by dividing the number of changing edges (p < 0.01) by the total number of network edges. The number of trials (for each gene) was equal to the number of edges. The number of “successes” was equal to the number of significantly changing edges (p < 0.01). Genes identified as DW were enriched in changing edges. For the analysis of network connectivity and Hub genes, we assessed DE, DV, and DW as shown in Tables S6, S8, and S9, respectively. GO annotation of affected modules is shown in Table S7. Further analysis of network modules enriched in Hub genes was performed using CNS‐related expression datasets within GeneMANIA (Ward‐Farley et al., 2010; https://genemania.org; datasets included in the analysis are listed in GeneMANIA as: Horvatta‐Barlow 2005, Zapata‐Barlow 2005, Akahoshi‐Ishii 2011; Jacob‐Benoist 2010; Guan‐Troyanskaya 2008; Zhang‐Hughes 2004 and Siddiqui‐Merra 2005).

Module Characterization

To assess the functional significance of all modules, gene ontology (GO) enrichment analysis using the GO‐stats R package (Ashburner et al., 2000; Falcon and Gentleman, 2007) was utilized. A graph decorrelation program (Alexa et al., 2006) was used due to the nested structure of GO terms. A ranking procedure was implemented using differential network results at both module and gene summary levels to develop a comprehensive gene screening procedure. In order to implement a ranking procedure, we incorporated differential network results at the module and gene summarization levels into an inclusive gene screening protocol. Modules enriched in gene or edge changes were the primary focus of further annotations. At the individual gene level, the emphasis was on module hubs with normalized intramodular connectivity above 0.8. The GOrilla algorithm (Eden et al., 2009) was utilized for illustrating a representation of GO annotation enrichment and for probing for annotation enrichment of selected groups of genes against a background set of all network genes.

Results

Dual‐Trait Selection

Figure 3 shows the response to selective breeding of the SOT and NOT lines at the third generation (S3). EtOH consumption and withdrawal severity were significantly different between SOT and NOTs. Consumption was significantly higher in the SOTs (6.8 ± 0.5 vs. 2.4 ± 0.3 g/kg, respectively; p ≤ 2.5 × 10−11). SOT preference and consumption scores were significantly different from NOT scores at each generation, F(1, 131) ≥ 12.7, p ≤ 4.6 × 10−4. Figure S1 shows EtOH consumption and preference for each selection generation. In a comparison of selection testing for both traits between S0 (founding population of F2 mice) and S3 mice, S3 SOT mice drank 3.45 g/kg more than S0, and NOT mice consumed 0.94 g/kg less than S0. Withdrawal severity at S3 was significantly higher in the NOTs compared to SOT (corrected AUC scores 4.6 ± 0.4 vs. 2.0 ± 0.3, respectively; p ≤ 2.5 × 10−11). Residual withdrawal scores were also significantly different between SOT and NOT lines (1.9 ± 0.4 vs. −1.7 ± 0.3, respectively; p ≤ 1.6 × 10−11). Withdrawal severity scores across the generations (S0 to S3) are shown in Fig. S2. By S3, withdrawal severity had decreased in SOT mice and increased in NOT mice by approximately the same extent (~50% from mean S0 values).

Fig. 3

Differences in phenotype due to selection. EtOH drinking and residual withdrawal scores in the S3 generation of SOT‐ and NOT‐selective lines relative to the SO means (solid line). (A) Consumption is expressed as average g/kg of 10% EtOH consumed per day. The gray bar indicates increased drinking in the SOT line, and black bar shows decreased drinking in the NOT line. (B) Withdrawal is expressed as the residual from a linear regression of postethanol injection HIC scores on baseline HICs. The black bar indicates increased withdrawal in NOT line, and the gray bar shows decreased withdrawal in SOT line. All values were significantly different from the S0 generation mean values (p’s < 0.05). A residual of 0 is indicated by the dotted line.

Effect of Selection on SOT and NOT Genotypes

Genotyping data are summarized in the MDS plot (Fig. 4). The data illustrate that SOTs and NOTs were genetically distinct at S4. The SOT mice were genetically closer to the B6 founder, while the NOT mice were genetically closer to the D2 founder. Two SOT mice plotted to the NOT domain (data not shown) and were removed from all subsequent analyses. Despite the segregation of the selected lines, significant genetic variation was retained.

Fig. 4

Multidimensional scaling (MDS) plot of genetic differences due to selection. Total genetic distance between samples was computed by summation, over markers, of the number of different alleles. The collection of distances is projected on 2 dimensions for visualization. The figure illustrates the genetic segregation of the SOT (green) and NOT (red) animals. As expected, SOT animals are genetically close to B6, while NOTs are closer to D2; the distance between B6 and D2 is the largest possible distance given the design of our F2 intercross and selective breeding.

Focusing on the SNP markers common in this study and Metten and colleagues (2014), there were 33 markers in this study that exceeded the critical LOD threshold of 10 (see Hitzemann et al., 2009 for details). Significant QTLs were detected on Chr 1, 2, 4, 5, 13, and 19 (see Table S1). The QTL on Chr 4 overlaps with the Chr 4 QTL found in Metten and colleagues (2014). The Chr 1 and 19 QTLs are located in regions that have been previously reported to contain withdrawal QTLs (Buck et al., 1997, 2002). The Chr 5 and 13 QTLs are new with this dataset. The Chr 2 QTL does not overlap with that described by Metten and colleagues (2014), even though among the SNP markers in common, they reported 85 markers that exceeded the LOD threshold with 58 of those on Chr 2. The QTLs on Chr 9 and 12 were not confirmed in this study.

Overview of Gene Expression Data

The ventral striatum transcriptome from a total of 41 SOT (19 females, 22 males) and 40 NOT (18 females and 22 males) mice was sequenced essentially as described elsewhere (e.g., Colville et al., 2017), yielding ~30 million mapped reads per sample. The data were culled to those genes with an average of 1 CPM (N = 14,507; Table S2). The means ± SD, fold change between SOT and NOT animals, and DE (p‐value and adjusted p‐value) are also found in Table S2. A total of 1,816 genes were detected as differentially expressed (FDR < 0.05); 1,079 genes had higher expression in the NOT mice and 737 genes had higher expression in the SOT mice. GO annotation of the differentially expressed genes (DEGs) revealed that for the NOT‐overexpressed genes there was a trend to enrichment in genes associated with nucleoside triphosphate biosynthetic process (unadjusted p < 4 × 10−5) and related processes (Table S3); none of the annotations met the minimal FDR threshold of 0.05. For the SOT‐overexpressed genes, there was a significant enrichment in genes associated with cell adhesion (FDR < 2 × 10−3) and a variety of membrane annotations including postsynaptic membrane (FDR < 3 × 10−3). The genes in this category are listed in Table S3 and include Comt, Dlg2, Dlgap2, Gabrb3, Glra2, Grid1, Grin2b, Grin2c, Grm3, Mpdz, Nbea, Ntrk2, and Pcdhb16. Of the 58 genes with a cell adhesion annotation, 23 were protocadherins.

The genes in Table S2 were entered into a consensus module network analysis, which parses the data from each selected line into network modules (Colville et al., 2017, 2018). The list of genes was culled to remove those genes that contribute the least to network connectivity (the bottom 10%). The culling procedure reduced the number of genes to 9,065, and it was this list of genes that was used for all subsequent analyses. This maneuver largely removes genes that are assigned to the “gray” module in the standard WGCNA (Langfelder and Horvath, 2008). The genes in Table S4 are assigned to color‐coded modules; color has no meaning. The consensus network was annotated for enrichment in genes associated with neurons, astrocytes, and oligodendrocytes (Table S5). Ten of the modules were enriched in neuronal genes; 2 of these modules (green and greenyellow) overlapped with enriched annotations for astrocytes and/or oligodendrocytes. Seven modules each were enriched in genes with an astrocyte or oligodendrocyte annotation; one of the modules overlapped (yellow).

Network—Differential Expression

Within the culled gene network, 764 genes showed a higher expression in the NOT versus SOT and 366 genes were overrepresented in the SOT versus NOT (Table S6). The number of DEGs compared to those in Table S2 was reduced by 29% for NOT > SOT and 50% for SOT > NOT. The culled DEGs were distant leaf network nodes which are the genes that are the least well connected to the network. Distant leaf nodes have, in general, reduced variance in gene expression, which facilitates the detection of DE.

The integration of the DE and network data focused on detecting DE nodes that showed a change in relative intramodular connectivity of ≥0.5 (see Colville et al., 2017). Relative connectivity was scaled from 0 to 1, and genes with a relative connectivity of ≥0.8 were considered hub nodes. For genes more highly expressed in the NOT line, 21 showed an increase in connectivity of >0.5 (Table S6). This gene list was then parsed into (i) the hub nodes (N = 8) and (ii) the nodes in the red (N = 8) plus turquoise (N = 4) modules (6 of these nodes were hub nodes). A summary of the hub nodes showing the strong effect of selection are listed in both Table 1 and Fig. 2. The annotation of the red plus turquoise grouping was somewhat stronger and that is reported here. The 12 genes were entered into GeneMANIA which detected that they were nested in a larger coexpression network structure of 31 genes (Fig. 5). Of these 31 genes, 12 were in the red module and 18 were in the turquoise module. As shown in Table S7, both of these modules were enriched in mitochondrial annotation. The 31 gene set (listed in the legend to Fig. 5 and Table S6) was also annotated using the GOrilla algorithm (Eden et al, 2009). The results were similar to those for the red and turquoise modules (see Table S7); significant enrichments were detected for NADH dehydrogenase complex assembly (FDR < 2 × 10−10), mitochondrial respiratory chain complex I assembly (FDR < 8 × 10−11), proton transmembrane transporter activity (FDR < 5 × 10−3), and mitochondrial protein complex (FDR < 8 × 10−22). The derived 31 gene lists contained 9 NADH–ubiquinone oxidoreductase flavoproteins and 3 mitochondrial ribosomal proteins.

Table 1

Differentially Expressed (DE), Differentially Variable (DV), and Differentially Wired (DW) Hub Nodes Strongly Affected by Selection

DE—NOT Connectivity > SOT Connectivity

Atp6v1f, Eif5a, Elof1, Guk1, Minos1, Nedd8, Ndufa8, Wdr45

DE—SOT Connectivity > NOT Connectivity

Ptp4a1

DV—NOT Connectivity > SOT Connectivity

Wdr45

DW—NOT Connectivity > SOT Connectivity

Eif5a, Elof1, Guk1, Nedd8, Ndufa8, Nfkbib, Minos1, Msrb2, Wdr45

DW—SOT Connectivity > NOT Connectivity

Gtf2a1, Ptpn4

Strong change = a change in relative intramodular connectivity of ≥0.5.

Fig. 5

GeneMANIA representation of mitochondrial gene cluster associated with selection. Twelve nodes in the red and turquoise modules (6 of which were hub genes) were entered into GeneMANIA, which detected that these 12 genes were nested within a larger co‐expression network containing 31 genes. Purple lines indicate gene co‐expression, blue lines indicated co‐localization, and orange lines indicated protein–protein interactions. The hatched circles identify the genes from our analysis, whereas the filled circles indicate the additional genes within the network identified by GeneMANIA. The genes in the network include the following: Anapc11, Atp5k, Atp5o, Atp6v1f, Cox5b, Cox7a2, Elof1, Gcsh, Guk1, Minos1, Mrpl48, Mrpl51, Mrps21, Naa38, Ndufa2, Ndufa3, Ndufa5, Ndufa8, Ndufb10, Ndufb5, Ndufb11, Ndufb6, Ndufc1, Ndufv2, Nedd8, Pam16, Sdhb, Timm8b, Uqcc2, Uqcr11, and Zswim7.

The genes more highly expressed in the NOT line showed enrichment (corrected p < 0.05) in 4 modules: black, gray60, purple, and red. The red module is noted above. Annotations for the black, gray60, and purple modules are found in Table S7. The purple module was also enriched in mitochondrial annotations, for example, mitochondrial part (FDR < 2 × 10−3).

For the genes more highly expressed in the SOT line (Table S6), only 2 (Ralgapa2 and Zhx1) showed a change in intramodular connectivity of ≥0.5 and neither was a hub node. The SOT highly expressed genes were enriched in 2 modules: greenyellow and yellow (Table S7). The greenyellow module was significantly enriched in annotations associated with the synapse including postsynaptic (FDR < 3 × 10−4). Genes with this annotation included Camk4, Gabbr1, Gabra4, Gabrg3, Gphn, Grid2, Grm3, Grm5 Itpr1, Nlg1, and Pde4b. None of the genes with a postsynaptic annotation were hub nodes.

Network—Differential Variation (DV)

Eighty‐four genes showed significant DV (Table S8); however, only one of these genes (Wdr45; see Table 1) showed a change in connectivity of ≥0.5. The DV genes were enriched in 2 modules (black and purple). The purple module is described above. Annotation for the black module was weak, but there was a trend to annotations associated with oxidoreductase activity.

Network—DW

Consistent with previous observations (Colville et al., 2017, 2018), selection had strong effects on DW; 1,717 genes were significantly affected (Table S9). Forty‐three genes showed a difference in relative intramodular connectivity of ≥0.5, and for 34 of these genes, connectivity was higher in the NOT line. Of these 34 genes, 10 were hub nodes (see Table 1). The 10 hub nodes were entered into GeneMANIA, which revealed the genes nested in a larger cluster of 31 genes (not shown). Twenty‐three of these genes were also found in the DE‐derived cluster. The annotation of the DW cluster was similar to that for the DE cluster, with strong annotations for mitochondrial‐related processes and components (Table S7). The DW genes were enriched in 2 network modules, red and turquoise, which have been described above. For the 9 genes with higher relative connectivity in the SOT line, no common annotation was detected.

Discussion

We recognize that there are numerous approaches to analyzing genome‐wide RNA‐Seq surveys of the brain transcriptome. The approach we have taken here is one that we have previously found useful (Colville et al., 2017, 2018; HItzemann et al., 2019). The focus is on 3 parameters (DE, DV, and DW), which can be generated completely independent of the network analyses. In the current study, DV was the least informative parameter. This differs markedly from the results in Colville and colleagues (2017, 2018) where a large number of genes exhibited sizable changes in variance associated with selection for EtOH preference. However, a key difference is that these studies used HS‐collaborative cross (HS‐CC) mice as the founders rather than B6xD2 F2s as founders. Using SNPs as a surrogate for genetic diversity, the HS‐CC is approximately an order of magnitude more diverse than the F2 animals (Roberts et al., 2007), which may account for the difference in the power of the DV parameter. We used the WGCNA to limit the field of investigation to those genes (nodes) with a detectable intramodular connectivity by discarding the distant leaf nodes. These “culled” nodes may be useful biomarkers, but the prediction would be that their manipulation would only have very limited effects on the SOT/NOT phenotypes. Rather we focused on those nodes that showed a large change in intramodular connectivity and were hub nodes. This approach was especially useful in the characterization of the nodes where connectivity was markedly greater in the NOT line for both DE and DW genes. The analysis revealed that from both perspectives, selection was associated with marked effects on mitochondrial genes and largely the same genes, which include clusters of NADH–ubiquinone oxidoreductase flavoproteins and mitochondrial ribosomal proteins.

Recent publications have highlighted gene expression and network co‐expression differences between D2/B6‐derived mouse strains exhibiting severe and mild EtOH withdrawal symptoms (Metten et al., 2014; Walter et al., 2017). Metten and colleagues (2014) found that 1 module in particular was significantly disrupted between SOT and NOT strains. They reported that this module (light pink) was enriched in genes with GO annotations that included MAP kinase signaling. However, further analysis of the data presented in the Supplemental Tables from the Metten and colleagues (2014) manuscript also revealed annotations including cellular components mitochondrion, respiratory chain, mitochondrial membrane part, mitochondrial respiratory chain, mitochondrial inner membrane, mitochondrial respiratory chain complex I, and biological processes electron transport chain and cell redox homeostasis. Walter and colleagues (2017) assessed gene expression differences in whole brain from reciprocal D2 and B6 background chromosome 1 congenic strains, which differed in EtOH withdrawal severity. This study found that there was enrichment in KEGG pathways which included oxidative phosphorylation and ribosome.

Network enrichment analysis of frontal cortex from alcohol‐exposed mice showed enrichment in 1 module for electron transport chain pathways (Nunez et al., 2013), suggesting that perturbations that exist as a predisposition to alcohol withdrawal severity might persist following EtOH exposure. Consistent with this idea, Berres and colleagues (2017) found that exposure to alcohol during chick embryo development resulted in DE in KEGG pathways that include ribosome and oxidative phosphorylation.

Several studies have assessed differential gene expression or networks using postmortem tissue from human alcoholics and controls (for review, see Warden and Mayfield, 2017). Liu and colleagues (2006) found DE in mitochondrial genes in the superior frontal cortex of alcoholics compared to controls. In a study using prefrontal cortical tissue, Zhang and colleagues (2014) found DEGs in GO categories that included oxidative reduction, mitochondrion, mitochondrial part, mitochondrial envelope, mitochondrial membrane, and mitochondrial inner membrane mitochondrial genes. Flatscher‐Bader and colleagues (2006) reported DEGs involved in oxidative phosphorylation and energy production in the prefrontal cortex but not in nucleus accumbens or ventral tegmental area. However, Mamdani and colleagues (2015) reported DEGs in GO annotations that included oxidative phosphorylation and mitochondrial function, in the nucleus accumbens. The postmortem human studies assessed self‐reported EtOH exposure or consumption, but not specifically withdrawal symptoms.

As noted above, the inverse reciprocal relationship between preference and withdrawal appears to diminish as 1 moves from relatively simple to more complex founder populations for selection. However, here we wish to emphasize that some of the genes associated with high preference in the current study (the SOT line) are very similar to the results for animals selected for high preference from genetically complex HS‐CC founders (Colville et al., 2017, 2018), even though how the similar genes were detected differs. Colville and colleagues (2017, 2018) observed that when comparing the high versus the low preference lines, there were relatively few DEGs in the 3 brain regions studied: the prelimbic cortex, the nucleus accumbens shell, and the central nucleus of the amygdala. The authors concluded was that the variance structure in the genetically diverse HS‐CC founders was not complimentary to detecting DEGs; that is, the changes in gene expression were relatively small when compared to the variance in gene expression. In contrast to the current study, DV was a key parameter in detecting the transcriptional effects of selection. Across the 3 brain regions, a similar network module was highly enriched in DV genes; the module was enriched in genes associated with cell adhesion and synaptic signaling. The cell adhesion genes included clustered protocadherins, and the synaptic annotation included genes associated with GABA and glutamate signaling and plasticity. The roles of the protocadherins in the adult brain are not entirely clear, although Suo and colleagues (2012) have shown that both the α and γ protocadherin clusters are involved in the inhibition of Pyk2 (protein tyrosine kinase 2), which results in the disinhibition of Rac1 (Ras‐related C3 botulinum toxin substrate 1) that can facilitate the proper assembly of dendritic spines (see figure 8 in Suo et al., 2012). Among the synaptic genes, Colville and colleagues (2018) emphasized the effect on Dlg2 and NMDA receptor subunits, for example, Grin2b. Dlg2, which encodes for PSD‐93, was a hub node in the commonly affected module. The current study engaged a similar set of genes, but the detection parameter was DE, not DV. It is of interest to note that none of these DE genes were hub nodes. In the absence of our previous data, we would have likely concluded that the cell adhesion and synaptic genes more highly expressed in the SOT line (compared to the NOT line) are of interest but probably not important selection drivers. However, such a hub‐centric conclusion seems premature and may suggest that it would be wise to employ multiple strategies for network construction and hub node detection.

Two of the genes more highly expressed in the SOT line, Mpdz and Nbea, are of special interest. Mpdz has been previously identified as a quantitative trait gene for EtOH withdrawal (Buck et al., 1997; Fehr et al., 2002; Milner et al., 2015; Shirley et al., 2004). It encodes the multiple PDZ domain–containing protein (MPDZ/MUPP1), which is a cytoplasmic scaffolding protein containing 13 PDZ domains. It is involved in cytoskeletal and cell signaling pathway organization and has an important role in synaptic plasticity (Jones et al., 2009; Krapivinsky et al., 2004). Nbea, which encodes for neurobeachin, had decreased expression and was hypermethylated in nucleus accumbens core of rhesus macaque chronically self‐administering EtOH (Cervera‐Juanes et al., 2017). These authors proposed that decreased Nbea expression was associated with a decrease in GABAergic tone and increased NMDA‐mediated hyperexcitability. Subsequently, they found that Nbea knock‐down mice exhibit increased EtOH consumption (Cervera‐Juanes and Cuzon‐Carlson, unpublished observations). Although in the current study we found an increase in Nbea expression as a risk factor, the data point to Nbea mediated receptor trafficking as being a key element in excessive EtOH consumption.

Metten and colleagues (2014) and the current study differ in the analytic approach (microarray vs. RNA‐Seq) for detecting the transcriptional effects of SOT/NOT selection. With some exceptions, the data presented here differ from the previous reported results. Bottomly et al. (2011) compared RNA‐Seq to both Illumina and Affymetrix microarrays for analyzing the striatal transcriptome. There were substantial differences in the detection of DE genes between the arrays and RNA‐Seq although the patterns of DE genes were relatively similar. These differences are likely associated with the marked differences in the variance structure of the arrays and RNA‐Seq. However, we believe that an alternative explanation should also be considered, namely that the 2 selections while methodologically similar and producing similar selected lines actually have engaged different sets of genes. We have noted this previously in the selection of the high‐drinking‐in the‐dark (HDID) mouse lines; the HDID‐1 and HDID‐2 lines are clearly genetically distinct (Iancu et al., 2013). A similar phenomenon may be at work when comparing Metten and colleagues (2014) to the current study. Evidence supporting this position is seen in the QTL analyses. While both studies detected a QTL on Chr 4, the other QTLs detected were different. Of interest, the current study detected previously reported QTLs on Chr 1 and 19 (Buck et al., 1997, 2002) that were not detected in Metten and colleagues (2014).

In conclusion, the current study was undertaken to complement the observations in Metten and colleagues (2014) principally by leveraging the numerous advantages of RNA‐Seq over microarray‐based data (e.g., Hitzemann et al., 2013, 2014) to detect the effects of the dual selection on the brain transcriptome. The strong mitochondrial signal that emerged extends earlier observations (Walter et al., 2017) that mitochondrial oxidative homeostasis is associated with the genetic vulnerability to EtOH withdrawal; the data are also consistent with the detection of the QTL on distal Chr 1. The DE genes between the SOT and NOT lines included Dlg2, Nbea, several glutamate receptor subunits, and several cell adhesion genes, in particular numerous cadherins and protocadherins. Our data cannot discriminate as whether these affects are associated with preference and/or withdrawal. However, similar results were obtained when selecting for EtOH preference (Colville et al., 2017, 2018). Importantly, the founder populations for our current and previous study differ markedly in terms of genetic diversity (see Roberts et al., 2007), suggesting the results are not limited to B6 x D2 genotypes.

Conflict of interest

The authors declare that they have no competing financial interests.

Author Contributions

RH and KJB performed the concept of manuscript and designed the study. DL performed the sample preparation. SE, DL, and KS involved in the behavioral testing. PD contributed to RNA‐Seq and WGCNA analysis and involved in the manuscript preparation. RH, LBK, and DL involved in the expression and additional molecular analyses.

Fig. S1. Ethanol consumption and preference for each selection generation.

Fig. S2. Ethanol withdrawal severity for each selection generation.

Click here for additional data file.

Table S1. Significant Quantitative Trait Loci detected.

Click here for additional data file.

Table S2. Gene expression above the threshold of >1 count per million.

Table S3. GO Annotation of differentially expressed (DE) genes.

Table S4. Network modules after the least connected genes were removed. Intramolecular connectivity was calculated for the Consensus, SOT and NOT Networks. The difference between SOT and NOT intramolecular connectivity is indicated as Δ Conn.

Table S5. Networks with cell‐type specific genes expression.

Table S6. Identification of Hub Nodes in DE Networks. Modules enriched for Hubs, enrichment of cell‐type specific genes expression and GeneMANIA derived gene clusters.

Table S7. GO Annotation for affected Network modules. GeneMANIA derived gene clusters for DE and DW Networks.

Table S8. Differential Variability (DV) in Networks, and enrichment of cell‐type specific expression.

Table S9. Identification of Hub Nodes in DW Networks. Modules enriched for Hubs, enrichment of cell‐type specific genes expression and the GeneMANIA derived gene clusters.

Click here for additional data file.