Genetic dissection of glutathione S-transferase omega-1: identification of novel downstream targets and Alzheimer's disease pathways.

Alzheimer's disease (AD) is affected by genetic factors. Polymorphisms in the glutathione S-transferase omega-1 (Gsto1) gene have been shown by genetic correlation analyses performed in different ethnic populations to be genetic risk factors for AD. Gene expression profile data from BXD recombinant inbred mice were used in combination with genetic and bioinformatic analyses to characterize the mechanisms underlying regulation of Gsto1 variation regulation and to identify network members that may contribute to AD risk or progression. Allele-specific assays confirmed that variation in Gsto1 expression is controlled by cis-expression quantitative trait loci. We found that Gsto1 mRNA levels were related to several central nervous system traits, such as glial acidic fibrillary protein levels in the caudate putamen, cortical gray matter volume, and hippocampus mossy fiber pathway volume. We identified 2168 genes whose expression was highly correlated with that of Gsto1. Some genes were enriched for the most common neurodegenerative diseases. Some Gsto1-related genes identified in this study had previously been identified as susceptibility genes for AD, such as APP, Grin2b, Ide, and Psenen. To evaluate the relationships between Gsto1 and candidate network members, we transfected astrocytes with Gsto1 siRNA and assessed the effect on putative downstream effectors. We confirmed that knockdown of Gsto1 had a significant influence on Pa2g4 expression, suggesting that Pa2g4 may be a downstream effector of Gsto1, and that both genes interact with other genes in a network during AD pathogenesis.

Introduction

Alzheimer’s disease (AD) is a common neurodegenerative disease that is caused by both environmental and genetic factors. Only a handful of genetic variants are consistently associated with increased susceptibility to either early-onset autosomal-dominant AD or late-onset AD (LOAD). However, LOAD is estimated to be highly heritable (~80%), and much of the phenotypic variance has yet to be explained (Bertram et al., 2008). Untangling the complicated genetics of LOAD is necessary to elucidate disease pathogenesis and identify potential targets for effective treatment and prevention.

A linkage site on chromosome 10q harboring two members of the omega class of glutathione S-transferase genes (GSTO1 and GSTO2) has been reported to be an important potential risk factor for AD. Studies (Li et al., 2003; Wongtrakul et al., 2018) have found significantly decreased GSTO1 expression levels in AD patients and identified a significant association between GSTO1 polymorphisms and the age at onset of AD. Since the initial report of the involvement of GSTO1 in AD pathogenesis (Li et al., 2003), other studies have been published that draw different conclusions (Kölsch et al., 2004; Nishimura et al., 2004; Wang, 2015; Umlauf et al., 2016; Wongtrakul et al., 2018). Various research groups have presented different opinions on the role of GSTO1 polymorphisms in the clinical risk of developing AD. For instance, Ozturk and colleagues (2005) suggested that there is no special association between GSTO1 polymorphisms and AD risk. However, a growing number of studies have provided evidence verifying that GSTO1 variants are risk factors for AD. The rare variant GSTO*E155del, which may decrease enzymatic activity, was linked to a risk of developing LOAD (Piacentini et al., 2012), as were the minor alleles GSTO1*A140D and GSTO2*N142D (Allen et al., 2012). On the basis of the biological functions of the Gsto1 enzyme and its putative association with LOAD, GSTO1 is an especially interesting candidate therapeutic target.

GSTO1 has been demonstrated to be indispensable to the pro-inflammatory response macrophages mount in response to stimulation with lipopolysaccharide that is mediated through Toll-like receptor 4 (Menon et al., 2015). After lipopolysaccharide stimulation, GSTO1-deficient macrophages failed to up-regulate nicotinamide adenine dinucleotide phosphate oxidase 1 expression and produce reactive oxygen species, and also failed to up-regulate several pro-inflammatory factors, including inducible nitric oxide synthase and cyclooxygenase-2 (Menon et al., 2015). GSTO1 has been shown to participate in sodium fluoride-induced cytotoxic effects in the mouse hippocampus cell line HT22 (Wang et al., 2019). Furthermore, inhibition of Gsto1 reduced the expression of factors associated with learning and memory, such as brain-derived neurotrophic factor and cyclic adenosine monophosphate response element binding protein, indicating that GSTO1 inhibition may be an effective treatment for sodium fluoride-induced damage to learning and memory (Wang et al., 2019). In contrast, Gsto1 knockout mice exhibited more severe inflammation than wild-type mice in an inflammatory bowel disease model (Menon et al., 2017). However, the role of the GSOT1-mediated inflammatory response in AD pathogenesis remains unclear. Astrocyte-mediated inflammation and oxidative stress play important roles in neurologic diseases, such as AD, Parkinson’s disease (PD), Huntington’s disease (HD), multiple sclerosis, and amyotrophic lateral sclerosis. In the present study, we attempt to validate candidate Gsto1 downstream genes in primary cultured astrocytes and predict the role of Gsto1-mediated gene regulatory pathways in neurodegenerative diseases such as AD.

In the current study, BXD recombinant inbred (RI) mice were used to dissect the genetic regulation of Gsto1 and to identify other pathway members. BXD RI mice were generated from two independent advanced intercrosses between C57BL/6J (B6) and DBA/2J (D2) progenitor strains. Briefly, the progenitors strains (B6 and D2) were mated, followed by more than 20 generations of sib-matings. This resulted in a panel of inbred strains with fixed genotypes at each locus, with the parental B6 and D2 alleles segregated among the strains. BXD RI mice have been well characterized, and comprehensive gene expression profile and behavioral trait data have been collected. In addition, substantial genotype data are available for BXD mice. Thus, genotype, gene expression profile, and phenotype data for the BXD RI strain, which are essential tools for genetic research, are already available and make this mouse strain suitable for systems genetics analysis, such as analysis of sequence polymorphisms and expression variations and correlation analysis with behavioral traits (Dai et al., 2009; Wang et al., 2010; Chen et al., 2018). BXD parental strains have been used to investigate genetic regulators of AD, due to the vast genetic differences between B6 and D2, such as differences in amyloid β (Aβ) levels and processing—B6 is less responsive to Aβ immunization compared with other mouse strains (Spooner et al., 2002)—and a significant difference in amyloid precursor protein (APP) expression was detected in the hippocampus (Wang et al., 2010), which indicated that genetic differences between the B6 and D2 strains may be risk factors for AD. In addition, BXD RI mice have been successfully used to identify a number of quantitative trait loci (QTLs) related to AD (Chen et al., 2018), PD (Dai et al., 2009), and anxiety and stress responses (Wang et al., 2012; Cai et al., 2020). In the current study, we integrated gene expression data from BXD RI mice with linkage analysis to dissect the mechanisms modulating expression of Gsto1 and novel downstream targets.

Materials and Methods

GeneNetwork datasets

The GeneNetwork website (www.genenetwork.org) is a public web source that stores phenotype and transcriptome data from various BXD RI mouse tissues, including the hippocampus. The gene expression data used in the current study can be accessed using the dataset “Hippocampus Consortium M430v2.0 (Jun06) RMA” that we generated through a collaborative effort for the purpose of studying neurodegenerative diseases (Dai et al., 2009; Wang et al., 2010; Chen et al., 2018). This dataset includes steady-state transcript abundance measurements collected from the hippocampus of the parent strains (C57BL/6J and DBA/2J), F1 hybrids (B6D2F1 and D2B6F1), and 67 BXD RI strains using the Affymetrix Mouse Genome array (Dai et al., 2009; Wang et al., 2010; Chen et al., 2018).

Quantitative trait loci mapping

To dissect the mechanisms underlying variation in Gsto1 expression, genome-wide QTL mapping was carried out using QTL Reaper, as previously described (Chesler et al., 2005; Dai et al., 2009; Lu et al., 2016; Cai et al., 2020). There are four options for QTL mapping on the GeneNetwork website: interval mapping, marker regression analysis, composite interval mapping, and pair-scan analysis. In this case, interval mapping was used to compute linkage maps for the entire genome. The log of odds (LOD) score was used to assert that a causal relation exists between a chromosomal location and a phenotypic variant, such as Gsto1 expression variation.

Analysis of allele specific expression by SNaPshot

Variable expression levels between Gsto1 alleles in the hybrid F1 offspring were quantified using a combination of quantitative real-time polymerase chain reaction (qPCR) and single-base primer extension (SNaPshot, Applied Biosystems, Foster City, CA, USA), as described previously (Ciobanu et al., 2010). In brief, we isolated RNA from the murine hippocampi and genomic DNA from reciprocal F1 hybrids (B6D2 F1 and D2B6 F1 hybrids). These DNA samples were used as controls and tested at the same time as the RNA pools. We selected the single nucleotide polymorphism rs4232251, located in the 3′UTR of Gsto1, to design SNaPshot extension primers (). Fold change in B and D allele expression in both genomic DNA and cDNA pools were calculated, and Student’s t-test was carried out to validate the polarity of parental alleles.

Correlation analysis

We calculated coefficients between Gsto1 and the above-mentioned genome-wide dataset (Hippocampus Consortium M430v2.0 (Jun06) RMA) by Pearson correlation analysis to generate a list of genetically correlated genes, referred to here as co-expressed genes. Genes with an expression level greater than eight and a statistically significant (P < 0.05) correlation with Gsto1 were selected for further analysis.

The literature was also searched to identify genes most likely to be linked to Gsto1 function. Literature coefficients that evaluate co-citations can be calculated online (www.genenetwork.org) and were used to examine the correlation coefficient (r) value of genes already described using similar terminology in published papers (Homayouni et al., 2005; Cai et al., 2020). Genes with an r value > 0.5 were considered to have a high literature correlation, and are referred to here as co-cited genes. The Chilibot system (http://chilibot.net) and the ALZGENE database (www.alzgene.org) were used to look for genes associated with AD that had been previously reported. Chilibot searches the PubMed literature database (abstracts) for specific relationships between proteins, genes, or keywords, such as Gsot1 and its co-expressed genes. ALZGENE is a collection of published genetic association studies related to Alzheimer’s disease, with random-effects meta-analyses for polymorphisms with genotype data in at least three case-control samples. For example, we can determine whether a Gsto1 co-expressed gene is associated with AD based on the information provided by ALZGENE.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was carried out using the DAVID tool (https://david.ncifcrf.gov/) to classify biological processes and the functional categories from a list of genes. Co-expressed (P < 0.05) and co-cited (P > 0.3) genes were selected as the input gene list.

Identification of Gsto1 downstream genes and associated phenotypes

Putative Gsto1 targets were identified as follows: (1) all genes whose expression level was greater than eight units and that had a significant trans-expression quantitative trait loci (eQTL) on chromosome 19 located in a 10-Mb region (between 42 and 52 Mb) close to Gsto1 were ranked; (2) the correlation coefficient between Gsto1 expression and the expression of the above trans-regulated genes was calculated, and only genes whose expression was highly correlated with that of Gsto1 were retained; and (3) a partial correlation analysis was performed to exclude the effect of linkage disequilibrium and to estimate the degree of correlation that remained (cis-regulated genes located near Gsto1 were used as control genes in this analysis).

Data from an RNA-sequencing study of freshly isolated astrocytes and microglia from APP/PS1 transgenic and wild-type (WT) mice were used to identify differentially expressed genes between the AD and WT groups. These data are accessible through GEO Series accession number GSE137028.

The GeneNetwork website contains extensive phenotypic datasets ranging from behavioral to morphological to pharmacological. To identify phenotypes associated with Gsto1 variation, we queried the BXD phenotype database in GeneNetwork, which contains nearly 3000 phenotypes, to look for the phenotypes that are most closely related to hippocampal expression of Gsto1 (probe set 1416531_at).

Gsto1 knockdown in primary astrocytes and qPCR

Primary cultures of cortical astrocytes were obtained from 40 newborn C57BL/6J mice (specific pathogen-free) from the Laboratory Animal Center of Nantong University in China (license No. SYXK (Su) 2017-0046). This study was approved by the Laboratory Animal Ethics Committee of Jiangsu Province, China (approval No. S20190728-301) on July 28, 2019. As previously described (Li et al., 2019; Yang et al., 2020), the newborn mice were euthanized by exposure to carbon dioxide (CO2) for 35 minutes. The CO2 displacement rate is 30–70%, which ensures that the rats have lost consciousness before pain can occur.

Each brain was isolated, the meninges were removed, and the cortices were dissected and chopped into fine pieces and then enzymatically digested using trypsin. Cells were cultured by inoculation into poly-L-lysine-coated flasks. Dulbecco’s modified Eagle’s medium (Gibco, Grand Island, NY, USA) containing 1% (v/v) antibiotics (penicillin/streptomycin; Invitrogen, Carlsbad, CA, USA) and 10% fetal bovine serum (Sigma-Adlrich, St. Louis, MO, USA) was used as the culture medium. The medium was changed every 3 days. Mixed glial cultures were maintained in culture for 2–3 weeks at 37°C in a humidified atmosphere composed of 5% CO2 and 95% air. After culturing for 11–14 days, the cells were shaken overnight at 260 r/min at 37°C, and then resuspended and re-seeded to remove microglia.

To observe the effect of changes in Gsto1 expression on its candidate downstream genes, Gsto1 was knocked out. A small interfering RNA (siRNA) duplex specific for mouse Gsto1 and a non-silencing negative control RNA were synthesized (Biomics Biotechnologies Co., Nantong, China) (). siRNA transfection was carried out in accordance with the manufacturer’s instructions. Briefly, when the cells reached 50% confluence, 30 nM siRNA and SuperFectinTM II reagent (Pufei Co., Ltd., Shanghai, China) were added to the appropriate groups. After 48 hours, mRNA was extracted from the cells in the different groups, and qPCR was used to estimate the siRNA knockdown efficiency. Three biological replicates were tested in this study. The 2–ΔΔCt method (Chen et al., 2018; Cai et al., 2020) was used to calculate relative expression levels.

Statistical analysis

Gene expression data were analyzed using GraphPad Prism 5.0 (GraphPad Software, Inc., La Jolla, CA, USA). A P-value of < 0.05 was considered statistically significant. Normally distributed data are represented as the mean ± standard error of mean. The data were analyzed using Student’s t-test. To map the eQTL for Gsto1 expression variation, interval mapping was used to generate linkage maps for the entire genome. The permutation test available on the GeneNetwork website was used to determine the approximate LOD score for the entire genome (Chesler and Williams, 2004). A P-value of < 0.05 was defined as the significance threshold, and P < 0.63 was defined as the suggestive threshold. Pearson’s correlation analysis was used to analyze sample correlation and literature correlation values on the GeneNetwork website. A P-value of < 0.05 was defined as a statistically significant correlation.

Results

Gsto1 levels vary widely in the hippocampus of BXD RI mice

The expression of Gsto1 on the Affymetrix MOE430V2 array available from GeneNetwork was interrogated with two probe sets. Both probe sets are well correlated in the hippocampus (r = 0.75, N(strain) = 71) and target different parts of the Gsto1 gene. Probe set 1416531_at targets the last two coding exons and the proximal 3′ untranslated region (UTR), and probe set 1456036_x_at targets a more distal region of the 3′ UTR. Probe set 1416531_at showed 1.6-fold expression variation and an average expression level of 10.6. Probe set 1456036_x_at exhibited 2.9-fold variation across BXD RI mice, with an average expression level of 10.5 ().

.

F1 hybrids and parental strains had intermediate mRNA expression levels. The y-axis denotes expression levels (mean ± SEM), which are normalized on a log2 scale. Strains are indicated along the x-axis. Gsto1: Glutathione S-transferase omega-1.

Variation in Gsto1 expression is regulated by local sequence variants

Next, we evaluated whether the variations in Gsto1 expression were due to gene variants segregating among the BXD population, using QTL mapping. A significant cis-eQTL (probe set 1416531_at, LOD score = 19, P < 0.01; probe set 1456036_x_at, LOD score = 15, P < 0.01) was detected in the hippocampus close to where Gsto1 is located (Chr19, Mb = 47.938787) (). We identified a region with 1.5 LOD score confidence for the hippocampal eQTL between 47.1 and 48 Mb, close to the location of Gsto1 itself. As shown in , Gsto1 was consistently detected in the cis-eQTL by at least one probe set for all expression datasets surveyed. In most cases, Gsto1 expression variation was highly correlated with B allele inheritance. The exceptions were the hypothalamus and amygdala, where higher expression was contributed by the D allele. This significant and robust cis-eQTL suggested that sequence variants within or near Gsto1 may induce expression variation.

.

Genetic regulation of Gsto1 expression probe sets 1416531_at (top panel) and 1456036_x_at (bottom panel) is shown for BXD hippocampi. The upper x-axis shows the chromosome (1–19, and X), the lower x-axis shows the location in megabases, and the location of Gsto1 is indicated by the triangle below the x-axis. The strength of association between gene expression and genotype is expressed as a log of odds (LOD) score and is marked by a solid blue line. The red and gray lines represent the significance and suggestive thresholds, respectively. The left y-axis shows the LOD score in blue, and the right y-axis shows the additive effect. The red and green lines show the effect of the D or B allele, respectively, on trait values. Gsto1: Glutathione S-transferase omega-1.

Summary of the genetic regulation of .

Intensity indicates strength of association, from a log of odds (LOD) score of 0 (black) to a suggestive (light colors) or significant (bold colors) LOD score. Cool blue indicates increased expression contributed by the B allele, and warm red indicates increased expression contributed by the D allele. Arrowheads show the location of Gsto1 on chromosome 19 (Chr 19). Amy: Amygdala; Cb: cerebellum; Gsto1: glutathione S-transferase omega-1; Hip: hippocampus; Hyp: hypothalamus; NAc: nucleus accumbens; PFC: prefrontal cortex; PG: pituitary gland; Str: striatum.

An allele-specific expression assay was carried out to validate the Gsto1 cis-eQTL. We detected a significant difference between B and D alleles in reciprocal F1 strains (P = 0.00016), with the presence of the B allele correlating with increased Gsto1 expression. These results definitively prove that sequence variants within Gsto1 or cis-regulatory regions cause differences in Gsto1 mRNA expression levels among BXD strains ().

Gsto1 co-expression and network analysis

We queried the above-mentioned hippocampus consortium dataset (Hippocampus Consortium M430v2.0 (Jun06) RMA) to identify transcripts that are co-expressed with Gsto1 and may function in similar biological networks. Pearson correlation analysis revealed that 5043 probe sets (corresponding to 4127 unique genes) correlated with Gsto1 (P-value of sample correlation < 0.005). Among the 5043 candidates, 2741 probe sets (representing 2168 unique genes) exhibited significant co-citation with Gsto1 (literature correlation > 0.3). We identified significant enrichment of the annotated genes in 37 KEGG pathways (), including AD, PD, and HD pathways, suggesting that these genes may participate together in the development of these neurodegenerative diseases. In addition, several cell signaling pathways (e.g., mitogen-activated protein kinase (MAPK) and adipocytokine) and metabolism-associated pathways (e.g., glutathione metabolism, alanine, aspartate and glutamate metabolism, glutathione metabolism and histidine metabolism) were also enriched for genes that were co-expressed and co-cited with Gsto1.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of genes that are co-expressed with .

Co-expressed and co-cited genes (P-value of sample correlation < 0.05, literature correlation > 0.3) were selected as the input gene list for KEGG pathway analysis. Thirty-seven KEGG pathways with significant enrichments were identified and sorted by enrichment score (–Log (P value)). Gsto1: Glutathione S-transferase omega-1.

We identified 17 transcripts that are indicated by at least one database (Chilibot or ALZGENE) to be involved in AD. These genes include APP (Guyant-Maréchal et al., 2007), Grin2b (Jiang and Jia, 2009), and Ide, which the ALZGENE database reports as participating in the pathogenesis of AD. In addition, several other genes, including Atp5g2, Casp8, Hsd17b10, Ndufb3, Ndufs3, Ndufs7, Ndufv1, Ndufv2, Plcb4, Ppp3ca, Ppp3r1, and Psenen, are identified as AD-associated genes in PubMed or the Chilibot database ().

Variations in Gsto1 expression lead to changes in downstream gene expression and phenotype

Gsto1 expression variation has been shown to be regulated by a cis-eQTL, suggesting that it is caused by sequence variants within or near Gsto1, indicating that Gsto1 could be a genetic modifier. We reasoned that transcripts regulated by a trans-eQTL located on Chr19 (42–52 Mb)—a region that precisely overlaps the Gsto1 cis-eQTL—could be downstream targets of Gsto1.

The following criteria were used to identify potential downstream targets of Gsto1: 1) genes with expression levels greater than 8 in the hippocampi of BXD mice; 2) genes whose expression variation correlated significantly (P < 0.05) with that of Gsto1; and 3) genes that are trans-regulated by the locus nearby Gsto1. In total, eight transcripts met the above requirements ().

Evaluation of .

(A) Eight transcripts are trans-regulated by a locus near the Gsto1 gene on Chr19 (arrowhead). The different colors indicate the intensity of the genetic association and the additive effect of the allele polarity. Cool and warm colors indicate that B and D alleles increase expression, respectively. (B–G) Gsto1 and its downstream candidate gene expression in the astrocytes and microglia of Alzheimer’s disease (AD) and wild-type (WT) mice. The data were derived from an RNA-sequencing study of freshly isolated astrocytes and microglia from APP/PS1 transgenic and wild-type mice. They can be accessed through GEO series accession number GSE137028. Compared with WT mice, Gsto1 was up-regulated in the astrocytes (B) and microglia (C) of AD mice, but the difference was not statistically significant. Expression of the Gsto1 downstream candidate genes Pa2g4 and Ncoa6ip did not change significantly in astrocytes (D, F) or microglia (E, G) in the AD group. Gsto1, Pa2g4, and Ncoa6ip exhibited similar trends in expression variation in AD mouse astrocytes. (H) Gsto1 expression was significantly down-regulated after siRNA treatment in primary cultured astrocytes (*P < 0.05 and **P < 0.01). (I) Gsto1, Pa2g4, and Ncoa6ip expression in primary astrocytes after Gsto1 siRNA transfection compared with the scrambled siRNA control group. As expected, Gsto1 levels were greatly reduced. Knockdown of Gsto1 resulted in an ~50% decrease in Pa2g4 levels (**P < 0.01) and had little effect on Ncoa6ip levels. Data are expressed as the mean ± SD from three separate experiments, and were analyzed by Student’s t-test. AD: Alzheimer’s disease; Gsto1: glutathione S-transferase omega-1; Pa2g4: proliferation-associated 2G4; WT: wild-type.

To address any confounding factors caused by linkage disequilibrium—the close physical proximity of several cis-modulated genes—at the Gsto1 locus, we carried out a partial correlation analysis to eliminate interference from 18 neighboring cis-regulated genes. After the partial correlation analysis, only Pa2g4 and Ncoa6ip retained their high correlation with Gsto1. As shown in Figure and , Gsto1 was up-regulated in the astrocytes and microglia of AD mice. The data were derived from a recent RNA-sequencing study of freshly isolated astrocytes and microglia from APP/PS1 transgenic and wild-type mice (Pan et al., 2020). The expression of the Gsto1 downstream candidate genes Pa2g4 and Ncoa6ip was not significantly altered in astrocytes (Figure and ) or microglia (Figure and ) in the AD group. In contrast, Gsto1, Pa2g4, and Ncoa6ip expression levels exhibited similar trends in variation in AD mouse astrocytes. To validate the potential Gsto1 downstream targets, we examined the expression levels of Pa2g4 and Ncoa6ip in primary astrocytes after treatment with Gsto1 siRNA. Gsto1 knockdown resulted in a significantly decrease in Pa2g4 expression, indicating a significant regulatory effect of Gsto1 on Pa2g4 expression. No significant difference was observed for Ncoa6ip expression (). These results indicate that Pa2g4 is a downstream member of the Gsto1 network.

Next, we investigated the effect of the marked variation in Gsto1 expression on the ~3000 published BXD phenotypes curated in GeneNetwork. Hippocampal expression of Gsto1 (1416531_at) is significantly correlated with several central nervous system traits, such as glial acidic fibrillary protein levels in the caudate putamen (Record 15158; n = 10; r = 0.87, P = 0.0004), cortical gray matter volume (Record 10754; n = 29; r = 0.60, P < 0.0004), and hippocampus mossy fiber pathway volume (Record 12588; n = 26; r = 0.60, P = 0.0008). However, only cortical gray matter volume maps back to the Gsto1 locus, suggesting that this phenotype is downstream of Gsto1 expression variation.

Discussion

There is evidence that Gsto1 defends against the negative effects of oxidative stress, inflammation, and exposure to toxins. Alterations in Gsto1 function may be a risk factor for AD (Kölsch et al., 2004). In the present study, we systematically analyzed the causes and consequences of variations in expression at the Gsto1 locus and identified genes that are co-expressed with and possibly downstream targets of Gsto1, which may themselves be involved in AD and other disease pathways. Our findings provide new insight into the possible functions of Gsto1 and its associated genes in AD.

Differential Gsto1 expression has been identified in multiple brain regions of BXD strains, which could be due in part to sequence variation between their parent strains (C57BL/6J and DBA/2J) and which segregates among the BXD RI strains, especially within cis-regulatory regions of the Gsto1 gene and even within Gsto1 itself. Once a cis-modulated Gsto1 transcript has been validated, meaning that differences in its expression are likely due to polymorphisms within the gene itself, it becomes a valuable molecular resource, as variations in its expression level can be used to identify other genes in the same network and associated phenotypes. In the current study, we identified 4127 unique genes that were highly co-expressed with Gsot1 and 2168 genes that were highly co-cited with Gsto1 in the hippocampus. Interestingly, this network is significantly enriched for the most common neurodegenerative diseases, including AD, PD, and HD (KEGG terms: mmu05010, mmu05012 and mmu05016), which suggests that common variants in these genes may contribute to all three of these conditions. In addition, the mmu00250 (alanine, aspartate and glutamate metabolism) and mmu00480 (glutathione metabolism) pathways were also enriched in this group, indicating that these genes may be involved in the alterations in glutathione-dependent enzyme activities and participate in the regulation of glutathione homeostasis, which has been increasingly implicated in the induction and progression AD, PD, and HD (Johnson et al., 2012). Among the 34 genes (corresponding to 76 probe sets) enriched in the AD pathway, 17 were indicated by at least one database to be involved in AD. These genes include: 1) APP, which encodes the protein APP; most mutations in APP increase the production of Aβ42 and contribute to AD pathology; 2) Grin2b, which encodes the NMDA receptor subunit NR2B; an important missense mutation in GRIN2B that has been identified in patients with AD, representing a potential key susceptibility locus for genetic predisposition to AD; 3) Ide, which encodes an insulin-degrading enzyme that may be potential therapeutic target, as the expression of this gene adversely affected by dichlorodiphenyltrichloroethane, which could inhibit the clearance and extracellular degradation of Aβ peptides (Li et al., 2015); and 4) Psenen, which encodes presenilin enhancer gamma secretase subunits that are necessary for the maturation and additional stabilization of the γ-secretase complex. In addition, Psenen depletion reduces the total level of γ-secretase complex, and Psenen overexpression increases Aβ levels (Schafer et al., 2015).

Linkage disequilibrium is an important factor affecting the accuracy of QTL mapping. In this study, partial correlation analysis was carried out to exclude the effects of linkage disequilibrium and to estimate the correlation that remains after controlling for selected control genes. Several downstream genes and phenotypes also highlighted a small number of intriguing Gsto1 genetic targets. Some of these genes and phenotypes, such as hippocampal and cortical morphometric phenotypes, have been previously demonstrated to be involved in AD pathology, while others, such as Pa2g4, have not. AD progression involves profound degeneration of cortical and hippocampal structures. Several central nervous system-related traits are correlated with Gsto1 and support a role for this gene in the brain structures that are most sensitive to AD, including cortical gray matter volume and hippocampus mossy fiber pathway volume. Both traits are positively correlated with Gsto1 expression, and cortical volume maps back to the Gsto1 locus. Pa2g4 was confirmed as a downstream target of Gsto1. Silencing of Gsto1 in astrocytes resulted in decreased expression of Pa2g4 mRNA, supporting a significant influence of Gsto1 on Pa2g4 expression. Pa2g4, also known as ErbB3 binding protein (Ebp1), is an RNA-binding protein that modulates translation initiation through inhibition of eIF2α phosphorylation (Squatrito et al., 2006). Pa2g4 is abundantly expressed in neocortex neurons that form early in development and regulates their morphology, thereby affecting cell adhesion molecule synthesis and thus indicating an important role for this protein in neurogenesis and neurodevelopment (Kraushar et al., 2021). The role of Pa2g4 in AD and neurodegenerative disorders is unknown and warrants further investigation.

This study had several limitations. The genetic and bioinformatics analyses that we carried out were preliminary. Although we identified genes that may be synergistically involved in AD pathogenesis with Gsto1, including genes that are co-expressed with Gsto1 and its downstream regulatory targets, such as Pa2g4, further investigation is needed to elucidate the relationships among Gsto1, Pa2g4, and the other genes in the proposed genetic network. Furthermore, it will also be important to elaborate the specific roles of Gsto1 and its co-expressed and downstream genes that are involved in AD pathogenesis, which can be further examined in AD transgenic mice, thus providing more information for AD prevention and therapy.

Overall, the findings reported here support that Gsto1 expression variation is regulated by a cis-eQTL in BXD mice. Genes that are co-expressed with Gsto1, some of which were already known to be involved in AD and some of which were identified here for the first time, may participate in AD via novel interactions with Gsto1. Additionally, the validated Gsto1 downstream target Pa2g4 may represent a novel target for AD treatment and prevention.

Table 1

Primers used in analysis of allelic-specific expression and qPCR, and chemically synthesized RNA duplexes

	Forward primer	Reverse primer	Extension primer
Gsto1 #	GCA AAT GCT GTT TCC CCT CA	CGC TTA GGC CAG AAG GTC AAA T	GAA AGG AGC CAG TGA GAA TAC TGG A
Gsto1 siRNA_1*	GGC ACG AAG UCA UCA AUA U dTdT	AUA UUG AUG ACU UCG UGC C dTdT	NA
Gsto1 siRNA_2**	CCU ACC GCG AGU ACU UAA A dTdT	UUU AAG UAC UCG CGG UAG G dTdT	NA
Gsto1 siRNA_3***	CUC CUU CUA UGG UUG AUU A dTdT	UAA UCA ACC AUA GAA GGA G dTdT	NA
Control siRNA****	UUC UCC GAA CGU GUC ACG UTT	ACG UGA CAC GUU CGG AGA ATT	NA
Gsto1 †	AAT GCT GTT TCC CCT CAC TG	GGG ATG GCA GTG AAG ACT GT	NA
Pa2g4 †	ATG TGG ATG GCT TCA TTG	GGT CTG TAG GGT TCT GGA TA	NA
Ncoa6ip †	TCT GAG CCT TGG AAC CTC C	TCA TTT GAG CCA GGG ATT T	NA
Actb †	CCA CCA CAG CTG AGA GGG AA	AGC CAC CGA TCC ACA CAG AG	NA

All primer sequences are shown from 5’ to 3’. # These primer pairs (forward and reverse) were used to amplify the 3’UTR of Gsto1, and the PCR product was analyzed using an extension primer. Subsequent extension with DNA polymerase added a single fluorescent triphosphate complementary to the nucleotide at the polymorphic site, and allele-specific expression was calculated. *, **, and *** indicate three pairs of chemically synthesized RNA duplexes that were specific for the mouse Gsto1 sequence. **** indicates a non-silencing negative control siRNA. † indicates primers used for qPCR. Actb: β-Actin; Gsto1: glutathione S-transferase omega-1; NA: not applicable; Pa2g4: proliferation-associated 2G4; qPCR: quantitative real-time polymerase chain reaction; siRNA: small interfering RNA; UTR: untranslated region.

Table 2

Validation of cis-eQTL by measurement of difference in allele-specific expression

Sample	Signal intensity (peak height)		Ratio	Fold difference	P-value of Student’s t-test
	B	D
B6D2 male RNA pool	5101	4908	1.03932		0.01675
B6D2 female RNA pool	3824	3578	1.06875		0.02888
D2B6 male RNA pool	3013	2614	1.15264		0.06169
D2B6 female RNA pool	1700	1571	1.08211	1.08571	0.03427
B6D2_genomic DNA sample 1	3132	4324	0.72433		–0.14006
B6D2_genomic DNA sample 2	3376	4612	0.732		–0.13549
B6D2_genomic DNA sample 3	3494	4803	0.72746		–0.13819
B6D2_genomic DNA sample 4	2820	3856	0.73133	0.72878	–0.13589
All RNA sample				1.49	0.00016

Using C57BL/6J mice as the maternal parent and DBA/2J as the paternal parent, the inbred offspring were D2B6; similarly, using DBA/2J as the maternal parent and C57BL/6J mice as the paternal parent, the inbred offspring were B6D2. Ratio: Ratio of the signal intensity of the B allele to the signal of the D allele; Fold difference: DNA samples were used as controls, and then calculate the B, D allele signal ratios in RNA samples to DNA samples. eQTL: Expression quantitative trait loci.

Table 3

Co-expressed genes enriched in the AD pathway are also highly co-cited with Gsto1 and AD

Symbol	Probe ID	Mean expression	Sample r	Sample P(r)	Lit Corr	Chilibot	AlzGene	PubMed
App	1458525_at	9.41	–0.31	0.01	0.31	Yes	Positive	Verified
App	1438373_at	5.19	0.27	0.03	0.31	Yes	Positive	Verified
Atp5f1	1426742_at	12.15	0.24	0.05	0.37	NA	NA	NA
Atp5g2	1415980_at	11.99	0.25	0.04	0.36	NA	NA	Verified
Capn2	1416257_at	10.01	0.27	0.03	0.35	NA	NA	NA
Casp8	1424552_at	7.16	0.35	0	0.33	Yes	Negative	Verified
Cox6a2	1417607_at	8.67	0.48	0	0.39	NA	NA	NA
Cox6a2	1457633_x_at	6.75	0.28	0.02	0.39	NA	NA	NA
Ern1	1438997_at	6.51	–0.28	0.02	0.33	NA	NA	NA
Grin2a	1434638_at	13.13	–0.29	0.02	0.31	NA	NA	NA
Grin2b	1442370_at	9.56	–0.28	0.02	0.33	Yes	Positive	Verified
Hsd17b10	1448286_at	10.06	0.34	0	0.58	Yes	NA	Verified
Ide	1423120_at	10.04	–0.33	0.01	0.59	Yes	Positive	Verified
Ndufa6	1448427_at	12.17	0.25	0.05	0.43	NA	Negative	NA
Ndufa7	1436567_a_at	11.89	0.35	0	0.41	NA	NA	NA
Ndufa7	1428360_x_at	12.39	0.35	0	0.41	NA	NA	NA
Ndufa7	1450818_a_at	12.57	0.33	0.01	0.41	NA	NA	NA
Ndufa7	1422976_x_at	11.69	0.3	0.01	0.41	NA	NA	NA
Ndufa7	1455749_x_at	8.97	0.26	0.03	0.41	NA	NA	NA
Ndufb3	1416547_at	12.03	0.29	0.02	0.35	NA	NA	Verified
Ndufb4	1428075_at	12.47	0.45	0	0.35	NA	NA	NA
Ndufb4	1428076_s_at	12.3	0.35	0	0.35	NA	NA	NA
Ndufb5	1448589_at	12.65	0.36	0	0.31	NA	NA	NA
Ndufb9	1452184_at	13.56	0.29	0.02	0.35	NA	NA	NA
Ndufs3	1423737_at	12.69	0.34	0	0.47	NA	NA	Verified
Ndufs4	1448959_at	11.19	0.31	0.01	0.49	NA	Trend	NA
Ndufs6	1433603_at	11.57	0.32	0.01	0.39	NA	NA	NA
Ndufs7	1451312_at	11.85	0.34	0	0.43	Yes	Negative	Verified
Ndufs8	1423907_a_at	12.74	0.29	0.02	0.44	NA	NA	NA
Ndufs8	1423908_at	10.48	0.28	0.02	0.44	NA	NA	NA
Ndufv1	1415967_at	12.22	0.27	0.03	0.46	NA	NA	Verified
Ndufv2	1428179_at	13.02	0.24	0.05	0.68	NA	NA	Verified
Plcb4	1446947_at	6.82	–0.3	0.01	0.32	NA	NA	Verified
Plcb4	1441531_at	5.26	–0.25	0.04	0.32	NA	NA	Verified
Ppp3ca	1440051_at	9.72	–0.41	0	0.36	NA	NA	Verified
Ppp3cb	1459378_at	5.12	–0.32	0.01	0.37	NA	NA	NA
Ppp3cb	1446149_at	5.08	–0.31	0.01	0.37	NA	NA	NA
Ppp3cc	1420743_a_at	7.61	0.3	0.01	0.55	NA	NA	NA
Ppp3r1	1421786_at	9.41	–0.33	0.01	0.42	NA	NA	Verified
Psenen	1415679_at	11.47	0.33	0.01	0.31	Yes	Positive	Verified
Sdha	1426688_at	12.97	0.27	0.02	0.51	NA	NA	NA
Sdha	1433293_at	7.56	–0.26	0.03	0.51	NA	NA	NA
Sdhc	1448630_a_at	10.79	–0.25	0.04	0.46	NA	NA	NA
Uqcrb	1416337_at	11.79	0.32	0.01	0.42	NA	NA	NA
Uqcrb	1430326_s_at	12.08	0.28	0.02	0.47	NA	NA	NA

Shading denotes genes whose expression in the hippocampus was below background levels (average log2 expression < 7). AD: Alzheimer’s disease; ALZGENE: field synopsis of genetic association studies in AD, “positive” indicates significant (P < 0.05) association in at least one of the performed analyses, and “negative” indicates no evidence for significant association, while “trend” indicates results in between “positive” and “negative”; Chilibot: scan genes for relations with AD; Gsto1: glutathione S-transferase omega-1; Lit Corr: literature correlation between the genes listed in the table and Gsto1 (probe set 1416531_at) based on the NCBI Rif; NA: not applicable; PubMed: verified co-citation with AD; Sample P(r): P-value of sample r; Sample r: genetic correlation between the expression of Gsto1 and the expression of all other transcripts across the mouse genome.