Genetic Dissection of the Gpnmb Network in the Eye

Abstract

Purpose.: To use a systematic genetics approach to investigate the regulation of Gpnmb, a gene that contributes to pigmentary dispersion syndrome (PDS) and pigmentary glaucoma (PG) in the DBA/2J (D2) mouse.

Methods.: Global patterns of gene expression were studied in whole eyes of a large family of BXD mouse strains (n = 67) generated by crossing the PDS- and PG-prone parent (DBA/2J) with a resistant strain (C57BL/6J). Quantitative trait locus (eQTL) mapping methods and gene set analysis were used to evaluate Gpnmb coexpression networks in wild-type and mutant cohorts.

Results.: The level of Gpnmb expression was associated with a highly significant cis-eQTL at the location of the gene itself. This autocontrol of Gpnmb is likely to be a direct consequence of the known premature stop codon in exon 4. Both gene ontology and coexpression network analyses demonstrated that the mutation in Gpnmb radically modified the set of genes with which Gpnmb expression is correlated. The covariates of wild-type Gpnmb are involved in biological processes including melanin synthesis and cell migration, whereas the covariates of mutant Gpnmb are involved in the biological processes of posttranslational modification, stress activation, and sensory processing.

Conclusions.: These results demonstrated that a systematic genetics approach provides a powerful tool for constructing coexpression networks that define the biological process categories within which similarly regulated genes function. The authors showed that the R150X mutation in Gpnmb dramatically modified its list of genetic covariates, which may explain the associated ocular pathology.

Glaucoma is a genetically heterogeneous disease that damages retinal ganglion cells and the optic nerve, resulting in visual field loss. Elevated intraocular pressure (IOP) is a significant risk factor, although it is not a requisite. Millions of people worldwide are affected by this disease, ^{1 –3} which is the second leading cause of blindness in the United States. ⁴ Although adult-onset glaucoma runs in families, no single mode of inheritance can adequately describe its transmission from generation to generation. The expression and severity of glaucoma are now known to depend on the interaction of multiple genes, age, and environmental factors. ⁵  

Pigmentary glaucoma (PG), a form of secondary open angle glaucoma, is the second most common form of glaucoma in young adults. In this disease pigment sloughs off the posterior iris and blocks the drainage of aqueous humor through the trabecular meshwork. In some cases elevated IOP leads to glaucoma, whereas in other cases it does not, ⁶ indicating pigmentary dispersion alone is insufficient to increase IOP or to damage the optic nerve. The DBA/2J (D2) strain has served as the preeminent mouse model for pigment dispersion syndrome (PDS) and PG. ⁷ The eyes of D2 mice have no obvious abnormalities at 3 months of age. However, by 6 months, iris depigmentation develops in a large percentage of D2 mice; by 8 months, an increase in IOP followed by damage to the optic nerve develops. ⁷  

Mutations in two genes—tyrosinase-related protein 1 gene (Tyrp1), which causes iris stromal atrophy, and glycoprotein neuromedin B-associated gene (Gpnmb), which is associated with pigment dispersion—are already known to contribute to pigment dispersion in the D2 mouse. ⁸ Although it is known that Tyrp1 encodes a membrane-bound melanosomal protein with enzymatic activity required for melanin synthesis, ⁹ relatively little is known about Gpnmb function. Gpnmb was first identified in human melanoma cell lines with low metastatic potential ¹⁰ ; however, its expression is not limited to cells of melanocyte lineage. ¹¹ This protein has also been detected in dendritic cells ¹² and a small number of other cell types. ¹³ Previous studies also demonstrated that the presence of the Gpnmb mutation combined with a mutation in Tyrp1, specifically on the D2 genetic background, increases the likelihood of PG. ^{14 –16} Despite these findings, the molecular mechanisms by which Gpnmb contributes to PDS and glaucoma are largely unknown. 

Recombinant inbred (RI) strains of mice are a useful resource with which to identify the genetic sources of variation in phenotype, in our case, the severity of glaucoma. The largest panel of these strains, the BXD family, consists of the inbred progeny of a cross between the wild-type normotensive C57BL/6J strain (B6 or B) that does not develop any form of glaucoma with the glaucoma-prone D2 strain. Collectively, these 80 BXD strains are among the largest RI sets, ^{17 –19} and they have been used extensively in genetic and genomic studies of the eye and central visual system. The initial set of 26 BXD strains was bred using a D2 strain that retained the wild-type Gpnmb allele. In contrast, the set of 50 BXD strains we generated used a D2 strain homozygous for the Gpnmb mutation. We have exploited these two branches of the BXD family to compare and contrast expression networks associated with Gpnmb. ^{17 –20}  

The purpose of this investigation was to use a systematic genetics approach to construct coexpression networks for both wild-type and mutant Gpnmb. Using expression QTL (eQTL) mapping, we were able to identify a cis-eQTL that modulates Gpnmb expression levels. Additionally, by segregating strains based on the presence or absence of the R150X mutation in Gpnmb, we were able to construct distinct and nonoverlapping coexpression networks of genes associated with the mutant and wild-type transcripts. 

Materials and Methods

BXD Strains

Throughout this study, mice were handled in a manner consistent with the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research and the Guide for the Care and Use of Laboratory Animals (Institute of Laboratory Animal Resources, Public Health Service Policy on Humane Care and Use of Laboratory Animals), and all studies were approved by the Animal Care and Use review board of the University of Tennessee Health Science Center. A total of 284 adult mice of both sexes were used for this study: 268 BXD mice (two male and two female from each of 67 BXD strains), eight mice from parental strains (two male and two female from B6 and D2 parents), and eight F1 hybrid mice (two male and two female from each F1 hybrid mating). Animals were weaned at 25 days of age and housed in same-sex cages with two to five mice per cage until the time of testing. Animals had free access to standard laboratory chow and water. Animals were maintained on a 12-hour light/12-hour dark cycle. Room temperature ranged from 20°C to 24°C. 

RNA Isolation and Microarray Hybridization

Naive mice of 67 BXD strains, their parental strains (B6 and D2), and F1 hybrids (B6D2F1 and D2B6F1) were killed by cervical dislocation between 2 and 4 months of age. ¹⁹ It is critical to note that this is well before any overt sign of glaucoma develops in D2. The expression data so critical in our work were, therefore, presymptomatic, and differences were presumably associated with liability differences (causal factors) rather than downstream reactive changes after tissue damage. Eyes were quickly removed and processed using our published protocols. ¹⁹ Total RNA samples were hybridized to microarrays (M430V2; Affymetrix, Santa Clara, CA). Posthybridization staining and washing were performed according to manufacturer's protocols. The arrays were scanned and images were quantified. Two arrays were prepared and analyzed for each strain and consisted of four eyes from both male and female mice. 

Microarray Data Normalization

Normalization was performed using the robust multichip average method using our published protocols. ²¹ This method, typically used when analyzing data from high-throughput microarray chips, involves adjusting for background noise, quantile normalization, and summarization. ²² Expression levels below 7 are typically close to background noise levels. The units for these measurements are arbitrary and are used in our study for measuring proportional changes in expression. The advantage of this modified Z score is that a twofold difference in expression level corresponds approximately to a 1-unit difference. 

Next-Generation RNA Sequencing

RNA sequencing using a genetic analysis platform (SOLiD Analyzer; ABI, Carlsbad, CA) was performed according to the manufacturer's instructions. Briefly, the following steps were performed: total RNA from eyes of B6 and D2 strains were extracted (RNeasy kit; Qiagen, Valencia, CA); total RNA was depleted of rRNA (RiboMinus eukaryote kit; Invitrogen, Carlsbad, CA); the libraries were prepared and bar coded (SOLiD Whole Transcriptome Analysis and SOLiD Transcriptome Multiplexing Kits; ABI; the multiplexed transcriptome libraries were amplified on magnetic beads by emulsion PCR, and the DNA sequences on the beads were then sequenced using a ligation-based assay (SOLiD 3 system; ABI); approximately 50 million oriented 35 nt reads from both B6 and D2 were mapped against known mouse RefSeq transcripts (www.ncbi.nlm.nih.gov/refseq) and the NCBI36/mm8 and NCBI37/mm9 whole genome assemblies (Parktek Genomics Suite; Parktek Inc., St. Louis, MO); and the expression level across the Gpnmb transcript was determined for both parental strains. All RNA-seq data for both strains across the entire genome can be viewed directly on a copy of the UCSC Genome Browser we maintain at http://ucscbrowser.genenetwork.org. 

Western Blot Analysis

Standard Western blot analysis protocols were used to quantify GPNMB protein levels in the eyes of representative BXD strains with low (i.e., BXD51, BXD55, and BXD98) and high (i.e., BXD22, BXD60, and BXD80) expression levels of Gpnmb. Reduced protein samples (25 μg/lane) were separated on high-performance gels (NuPAGE 4–12%; Invitrogen). Anti-GPNMB (R&D Systems, Minneapolis, MN) was used at 0.2 μg/mL, and immunopositive bands were visualized using an horseradish peroxidase–conjugated secondary antibody and substrate (SuperSignal West Femto Maximum Sensitivity Substrate; Pierce, Rockford, IL). Blots were scanned using a Kodak image station (2000MM; Eastman Kodak, Rochester, NY). 

Expression QTL (eQTL) Mapping and Heritability Calculation

On the Affymetrix array (MOE 430 version 2; Affymetrix), Gpnmb was represented by one probe set, 1448303_at, located on chromosome 6 at approximately 49 Mb. Our previous experience has shown us that the SNP position and the number of SNPs overlapping probe sets can have a significant impact on the detection of cis-acting expression differences (see Ref. 23 for an example). Therefore, before we performed the eQTL analysis, we looked for SNPs in the region of Gpnmb to which the probe set 1448303_at corresponded and determined that no SNPs were present in that region (data not shown). 

QTL mapping was carried out using GeneNetwork (www.genenetwork.org) ²⁴ and the Haley-Knott regression method. ²⁵ BXD24 was excluded from the study because the substrain we received from The Jackson Laboratory (Bar Harbor, ME; stock ID 000031), which is now known as BXD24/TyJ-Cep290rd16/J, has undergone retinal degeneration because of a spontaneous mutation in Cep290. ²⁶ Simple interval mapping was carried out at regular chromosomal intervals to identify eQTLs that regulate Gpnmb expression. Significance levels and confidence intervals were estimated by permutation and bootstrapping, respectively. We also used composite interval mapping to factor out genetic variance associated with major eQTLs and to detect potential secondary eQTLs that might otherwise be masked. Each of these analyses produced a likelihood ratio statistic (LRS) score, providing us with a measure of the strength of linkage between our observed phenotype—in this case, Gpnmb expression level variation—and genetic markers. Candidate genes were identified based on associated LRS scores, mean expression levels, and the presence of polymorphisms. 

Heritability of the expression level of Gpnmb in each data set was calculated using the formula of Hegmann and Possidente ²⁷ : h ² = 0.5Vg/(0.5Vg + Ve), where h ² is the heritability, Vg is the genetic variance, and Ve is the environmental variance. The factor of 0.5 in this ratio is applied to adjust for the twofold increase of genetic variance among inbred strains relative to outbred populations. ¹⁹  

Correlation and Heat Map Analyses, Gene Ontology Tree Machine, and GeneNetwork Construction

A key aim of this work was to identify genes that were coregulated with Gpnmb and shared a functional relationship. A secondary goal was to determine whether Gpnmb was involved in modulation of the immune system in the whole eye. As a first step in this process, we segregated our expression data into two subsets based on the presence or absence of the R150X mutation in Gpnmb. Analyses were performed on each subset individually. The transcript level of our gene, Gpnmb, was compared with all 43,000 probe sets to produce a set of genetically correlated transcripts. Genetic correlative analysis was calculated using Spearman's rank correlations, which was computed using the tools at GeneNetwork. 

After removing Riken clones, putative intergenic sequences, and poorly annotated and nontranslated genes, the remaining list of transcripts with mean expression levels greater than 7 were analyzed by Gene Ontology (GO) enrichment analysis using WebGestalt (http://bioinfo.vanderbilt.edu/webgestalt). ²⁸ WebGestalt allows users to input lists of highly correlated genes through the Web interface, identifies GO terms that are significantly associated with the input gene lists, and visualizes the enriched GO terms in a directed acyclic graph. For the enrichment analysis, all probe sets on the array (M430v2; Affymetrix) were included in a reference set for the hypergeometric test in GO enrichment analysis. P values generated from the hypergeometric test were further adjusted to account for multiple comparisons using the Benjamini and Hochberg correction, ²⁹ as implemented in Gene Ontology Tree Machine. 

The top candidates from the statistically significant GO categories were used to construct network graphs of genes whose expression levels were highly correlated with Gpnmb. Each node in the graph represented an individual transcript, and the thickness of the lines joining nodes encoded Pearson correlation coefficients. The top correlates were also used to construct heat maps through the heat map link on GeneNetwork. 

As a positive control to demonstrate the ability of GO enrichment and correlative network analyses to detect the interactions with the immune system, identical analyses were performed using leukocyte surface antigen CD53. CD53 contributes to the transduction of CD2-generated signals in T cells and natural killer cells and, therefore, has direct functional relevance to immune-mediated events. ³⁰  

Results

Variation of Gpnmb Transcript and Protein Levels in Eyes of BXD Mice

The transcript expression level of Gpnmb varied significantly among BXD strains. The average expression in the BXD strains was 10.9 ± 0.2 (mean ± SEM) and ranged from a low of 7.9 ± 0.3 in BXD50 to a high of 12.2 ± 0.1 in BXD32 (Fig. 1A). In the B6 parent, Gpnmb had an average expression level of 12.0 ± 0.1. This was significantly greater than that of the mutant D2 parent (8.2 ± 0.1; P < 0.01). The original BXD lines were generated in two phases at The Jackson Laboratory. All but two of those BXD lines have high expression levels of Gpnmb, indicating that they have a wild-type transcript that is able to hybridize to the 1448303_at probe. The two lines with low Gpnmb levels (BXD40 and BXD42) were generated after the inadvertent fixation of a spontaneous mutation in Gpnmb arose in the D2 stock colony. ³¹ All the new lines generated at UTHSC since the late 1990s segregate for the wild-type B alleles and the mutant D alleles of Gpnmb. Therefore, some lines have high expression levels of Gpnmb and others have low expression levels. These data were used to segregate our expression data into two subsets: one with wild-type Gpnmb and one with the R150X mutation. 

Figure 1.

 

Expression levels of Gpnmb across strains. (A) Rank-ordered mean Gpnmb levels across 67 BXD RI strains, their parental strains, and F1 crosses. N, two male and two female mice per strain. Values denote normalized relative expression levels on a log₂ scale (mean ± SEM). All but two RI strains (BXD40 and BXD42) generated by the Jackson Laboratory (JAX) have high levels of Gpnmb. The RI lines we generated at UTHSC express both high and low levels of Gpnmb. (B) RNA-seq data showing the level of Gpnmb across all 11 exons from D2 (bottom, white bar) and B6 (middle, black bar) parental strains. The level of Gpnmb message was greater in the B6 mouse across the entire transcript. (C) Whole eye lysates from BXD lines with high and low mRNA expression levels underwent Western blot analysis with an anti-GPNMB antibody. This representative blot shows the absence GPNMB protein at the expected molecular weight in lines with low mRNA expression levels.

Genomewide RNA sequencing data from eyes have been generated in our laboratory using RNA sequencing with a genetic analysis platform (SOLiD Analyzer; ABI). Counts were mapped to the B6 reference genome to reveal the difference in the expression of genes of importance in the eye. The expression level of Gpnmb was assessed across the entire transcript. Data illustrated that the expression of Gpnmb was abundant in B6 mice. In contrast, the expression of Gpnmb in D2 was low across the entire transcript from the 5′ to the 3′ UTR regions. Specifically, the cumulative number of counts for Gpnmb from B6 mice is 3450, whereas that from D2 mice is 424 (Fig. 1B). 

Western blot analysis showed that in BXD strains with the R150X mutation, GPNMB levels were nondetectable, whereas in those with wild-type Gpnmb, protein levels were high. This finding was consistent across each of the three representative low- and high-expressing lines. A representative Western blot is shown in Figure 1C. 

eQTL Mapping for Gpnmb and Heritability Calculation

Variation in Gpnmb expression in the eye among the UTHSC BXD strains is associated with a highly significant eQTL with an LRS of 110.7 [logarithm (base 10) of odds] score of 24) on chromosome 6 within 3 Mb of the location of its cognate gene (Fig. 2). Composite interval mapping controlling for the Gpnmb interval revealed no secondary loci modulating Gpnmb expression (data not shown but easily recomputed using the BXD eye expression database in GeneNetwork; see Ref. 19). Heritability was calculated for both BXD subfamily (BXD1–42 and BXD43–103) data sets. Among BXD strains with wild-type and mutant Gpnmb, heritability of mRNA variation was 20.3% and 30.0%, respectively. 

Figure 2.

 

Graphic illustration of LRS scores for Gpnmb expression in eyes of BXD mice. A significant eQTL for Gpnmb is present on chromosome 6 at the location of the gene itself. Dark gray trace: LRS scores across the genome; horizontal lines: transcript-specific significance thresholds for significant (LRS∼17) and suggestive (LRS∼10) based on results of 1000 permutations of the original trait data.

Biological Process Enrichment and Coexpression Networks

After filtering the top 500 transcripts with expression levels that correlated with wild-type Gpnmb, a list of 406 transcripts was submitted for GO enrichment analysis, and a directed acyclic graph that grouped the list of transcripts according to biological processes was produced (Fig. 3A). Ten of the 37 (27%) categories in the tree were statistically significant for enriched transcripts. The categories that were statistically significant formed two distinct branches of the GO enrichment graph. The first branch included the categories of “localization” (103 genes), “localization of cell” (29 genes), “cell motion” (29 genes), “regulation of cell motion” (15 genes), “cell motility” (24 genes), “cell migration” (24 genes), and “regulation of cell migration” (14 genes). The second branch included “melanin metabolic process” (5 genes), “melanin biosynthetic process” (5 genes), and “melanin biosynthetic process from tyrosine” (4 genes). Supplementary Appendix S1 lists the genes in each enriched GO category. 

Figure 3.

 

Genetic associations with wild-type Gpnmb. (A) In BXD RI mice with wild-type Gpnmb, GO enrichment analysis illustrates that most biological processes to which transcripts correlated with Gpnmb expression belong include cell motion, cell migration, and melanin synthesis. Enriched categories reaching statistical significance are indicated with a bold outline. (B) Genetic coexpression network generated from genes correlated with wild-type Gpnmb shows that the expression level of Gpnmb was directly linked to five nodes. Other correlations with Gpnmb within the network occurred through intermediary transcripts. Each transcript is shown as a node, and the Pearson correlation coefficient between nodes is indicated by a line. Bold lines: coefficients of 1.0 to 0.7; normal lines: coefficients of 0.7 to 0.5; dashed lines: coefficients of 0.5 to 0.0.

In generating the network graph composed of genes contained within the statistically significant GO enrichment, we limited ourselves to the significant categories with the exception of “localization.” The network graph function of GeneNetwork is limited to 100 traits, and that category alone contained 103 transcripts; most of the genes in that category were replicated in the other significant categories along that same branch of the biological process tree. The network graph indicated that the expression level of Gpnmb was directly linked to five nodes. Other correlations with Gpnmb within the network occurred through intermediary transcripts. The expression levels of the five genes involved in melanin biosynthesis or metabolism (Dct, Si, Tyr, Slc45a2, Myo5) were correlated with each other and with other genes involved in cell migration and motion (Fig. 3B), indicating the coregulation of both branches of the GO enrichment analysis. 

In BXD RI mice with a mutant allele of Gpnmb, 286 of the top 500 transcripts genetically correlated with Gpnmb remained after the data set was filtered. The GO enrichment analysis generated from this data set contained 35 categories, nine (26%) of which were statistically significant (Fig. 4A). The graph contained three branches with distinct biological functions. The first branch contained the following statistically significant categories: “biopolymer modification” (42 genes), “cellular protein metabolic process” (54 genes), “protein modification process” (41 genes), “posttranslational protein modification” (38 genes), “phosphorus metabolic process” (32 genes), “phosphate metabolic process” (32 genes), and “protein amino acid phosphorylation” (26 genes). Each of these contained various kinases or nucleotide exchangers involved intracellular signaling. The second branch contained “stress-activated protein kinase signaling pathway” (6 genes), and the third branch contained “sensory perception of light stimulus” (8 genes) and “visual perception” (8 genes). Five of the six genes in the “stress-activated protein kinase signaling pathway” were shared with statistically significant pathways of the first branch of the tree, suggesting significant functional overlap of these two branches. Interestingly, none of these statistically significant biological process categories were shared by the GO enrichment graph generated by genetic correlates of wild-type Gpnmb (compare Figs. 3A and 4A). Supplementary Appendix S2 lists the genes in each enriched GO category. The network graph generated from the list of transcripts within the enriched GO categories illustrated that Gpnmb was directly linked to 10 nodes in the network. All other interactions occurred through intermediaries (Fig. 4B). Because of the known contribution of both Gpnmb and Tyrp1 to PDS and PG, we specifically looked for the coregulation of both genes. Surprisingly, Tyrp1 was not included in the list of top 2000 correlates. 

Figure 4.

 

Genetic associations with mutant Gpnmb. (A) In BXD RI mice with mutant Gpnmb, GO enrichment analysis illustrates that the majority of biological processes to which transcripts correlated with Gpnmb expression belong include posttranslational protein modifications, stress-activated kinase signaling, and visual perception. Enriched categories reaching statistical significance are indicated with a bold outline. (B) A genetic coexpression network generated from genes correlated with mutant Gpnmb shows that Gpnmb was directly linked to 10 nodes in the network. All other interactions occurred through intermediaries. Each transcript is shown as a node, and the Pearson correlation coefficient between nodes is indicated by a line. Bold lines: coefficients of 1.0 to 0.7; normal lines: coefficients of 0.7 to 0.5; dashed lines: coefficients of 0.5 to 0.0.

The heat maps generated from the lists of top correlates show very different patterns. In the heat map of correlates of wild-type Gpnmb, tight bands of correlates are present on chromosomes 5, 9, and 19, and more dispersed bands are present on chromosomes 8 and 10 (Fig. 5A). In contrast, the data generated from the mutant strains show tight bands on chromosomes 4, 5, and 8. 

Figure 5.

 

Heat maps of Gpnmb covariates. (A) In BXD RI mice with wild-type Gpnmb, tight bands of correlates are present on chromosomes 5, 9, and 19, and more dispersed bands are present on chromosomes 8 and 10. (B) In RI strains with mutant Gpnmb, the heat map shows tight bands of correlates on chromosomes 4, 5, and 8.

The Gpnmb networks generated using the BXD strains that were stratified based on the presence or absence of the R150X mutation are specific to the eye. When similar analyses were performed on other tissues, completely different networks were generated. For example, in the lung, in which Gpnmb is expressed at moderate levels, enriched GO categories coregulated with wild-type Gpnmb included genes involved in sphingoid biosynthetic processes, whereas those coregulated with mutant Gpnmb included ATP biosynthetic process and RNA splicing. 

As a positive control, we performed identical analyses with Cd53, if it was present, to demonstrate the robustness of GeneNetwork and its ability to detect networks of coregulated genes involving the immune system. Cd53 (leukocyte surface antigen) is a tetraspanin involved in the immune system. Through its interactions with integrins, Cd53 mediates signal transduction in T cells and natural killer cells.³⁰ In BXD RI mice with wild-type alleles of Gpnmb, 368 genes remained after filtering. The GO enrichment graph produced from these transcripts that correlated with Cd53 expression contained 10 biological process categories with significant P values (10 of 62 [16%] total categories; Supplementary Appendix S3A). These categories were highly enriched for “immune system process” (43 genes), “regulation of immune system process” (25 genes), “positive regulation of immune system process” (21 genes), “immune response” (32 genes), “regulation of immune response” (18 genes), “positive regulation of immune response” (21 genes), “activation of immune response” (13 genes), “phagocytosis engulfment” (6 genes), “response to external stimulus” (40 genes), and “response to wounding” (29 genes). The GO enrichment graph generated from the genes within the statistically significant biological process categories indicated that Cd53 was centrally located in the network and was directly linked to many of the genes at high levels of correlation (Fig. 5B). Supplementary Appendix S4 lists the genes in each enriched GO category. 

The strong relationship between Cd53 and the immune system remained prominent in BXD RI mice with Gpnmb^R150X (Supplementary Appendix S5). Ten of the 53 categories (19%) in the tree reached statistical significance: “antigen processing and presentation” (9 genes), “antigen processing and presentation of exogenous antigens” (7 genes), “antigen processing and presentation of peptide antigen” (6 genes), “antigen processing and presentation of exogenous peptide antigen” (6 genes), “phagocytosis engulfment” 5 genes), “dsRNA fragmentation” (21 genes), “production of small RNA involved in gene silencing by RNA” (21 genes), “regulation of production of small RNA involved in gene silencing by RNA (21 genes),” and “regulation of response to stimulus” (21 genes). As expected, nearly all categories in both GO enrichment graphs generated from the two Gpnmb data sets were shared, and one (“phagocytosis engulfment”) was statistically significant. These data indicate that Cd53 is highly correlated with many biological processes of the immune system in BXD RI mice carrying both wild-type and mutant alleles of Gpnmb. These data also demonstrate that GeneNetwork can readily detect interactions of a gene with the immune response if it is present. Supplementary Appendix S6 lists the genes in each enriched GO category. 

Discussion

When considering the causes of blindness, it is impossible to ignore the importance of glaucoma. Although several studies have elucidated the genetic mutations resulting in some forms of glaucoma (summarized at http://neibank.nei.nih.gov/cgi-bin/CDL/showDiseaseLoci.cgi?dt=Glaucoma), most of these studies have considered individual genes rather than associated networks of coregulated genes. The systematic genetics approach has proven to be a powerful tool for identifying candidate genes and for the construction of genetic networks that regulate nervous system function. ²⁴ We combined this successful methodology with our powerful genetic resource of RI strains of mice to reveal the network of genes that covary with Gpnmb, a gene that contributes to PDS and PG in the D2 mouse. 

In our study, we excluded the BXD24 strain from all analyses given that the Cep290 mutation in its genome has a very strong effect on the relative expression levels of all transcripts in the eye. This is because the Cep290 mutation causes a total loss of photoreceptors before 2 months of age, thus distorting the normalization of all transcript levels in the eye. In contrast, strains with mutations in Gpnmb maintain their entire complement of ocular cells, and, therefore, the normalization of transcript levels is not grossly distorted. In fact, the slight variations in relative transcript levels between BXD strains with wild-type versus mutant Gpnmb allowed us to determine with great sensitivity other transcripts whose levels are modulated by the presence of the R150X mutation in Gpnmb. 

The Gpnmb gene is predicted to encode a transmembrane protein with homology to the melanosomal protein silver (Si). The GPNMB protein encoded by this gene is an intracellular, endosomal/melanosomal protein that is predicted to be important for the biosynthesis of melanin and the development of the retinal pigment epithelium and the iris, the pigmented tissues within the eye. ¹³ It is also expressed by low-metastasizing melanoma cells. ¹⁰ Although this knowledge is certainly useful, the focus is narrow and may overlook coregulated genes that are altered by SNPs or mutations in Gpnmb. Our study sought to expand on this functional knowledge by identifying loci that regulate Gpnmb expression and other genes whose expression levels covary with Gpnmb. Because previous investigations have suggested that Gpnmb plays a significant role in the immune response in the eye, specifically anterior chamber-associated immune deviation, ¹⁵ one of our goals was to determine whether Gpnmb participated in immunomodulatory pathways in the eye. In this investigation, we were unable to identify any link between Gpnmb expression levels and the regulation of immune responses. However, this finding does not rule out a biological interaction at the protein level. 

In our expression analysis, the abundance of a transcript was treated as a quantitative trait, allowing us to perform conventional linkage analysis to find the genetic loci that affect transcript abundance. However, transcript levels for a given gene differ from traditional phenotypes in that each transcript has a corresponding encoding gene with a known position in the genome. As mapping studies reveal the location of QTLs, their expression can be classified as cis (acting within 3 Mb of the location of the gene itself) or trans (acting elsewhere in the genome). Our investigation revealed a very large cis-eQTL (LRS score of 110.7) that regulates the expression level of Gpnmb. Comparison of this score with those of all 45,101 probes on the Affymetrix chip (M430V2) revealed that only 427 probes have cis-eQTLs of 110.7 or greater. This places the cis-eQTL in the upper 0.95% of cis-eQTLs in the mouse genome. 

Binding of a transcription factor to the promoter region of a gene is traditionally regarded as the main mechanism by which gene expression is regulated. Recent studies ^32,33 have shown, however, that polymorphisms within both intronic and exonic regions also play an important role in regulating gene expression. In our study, we found multiple SNPs in the 3′-UTR, intronic, and exonic regions of the Gpnmb gene (Supplementary Appendix S7). One of the SNPs encoded for the stop codon in exon 4, previously reported by Anderson et al. ³⁴ The SNP at 48.995384 Mb has been previously documented to introduce a stop codon that encodes for a truncated protein of 150 amino acids, which is less than half the full-length protein. It was predicted that the mutant mRNA undergoes nonsense-mediated decay ³⁴ and the transcript is not translated into protein, thus making it a loss-of-function mutation. Our RNASeq data demonstrate that the D2 mouse has very low levels of expression across the entire transcript, which supports the hypothesis that the Gpnmb^R150X mutant transcript is degraded and not expressed. Our Western blot data further corroborate this prediction and demonstrate that in BXD lines with low transcript levels of Gpnmb, no GPNMB protein is expressed, similar to what was documented in the D2 parental strain. ¹⁵  

To reveal the possible mechanisms by which Gpnmb functions in the eye, we stratified the RI lines based on the presence or absence of the Gpnmb^R150X mutation and performed analyses on each data set separately. Importantly, control of the expression levels of genetic correlates was heritable for strains with both wild-type and mutant alleles. Interestingly, there was no overlap of biological process categories in the GO enrichment graphs generated from transcripts correlated with wild-type or mutant Gpnmb. Furthermore, there was no overlap in any of the chromosomal regions in which Gpnmb correlates are localized, as shown by the heat maps. These data demonstrate that the loss-of-function mutation dissociates Gpnmb from its normal biological functions and molecular networks, which may result in abnormal interactions at both the transcriptional and the protein levels. 

Given its localization in pigmented tissues, ^13,15 it was not unexpected that wild-type Gpnmb correlated with other genes (i.e., Dct and Si) involved in melanin synthesis. Our finding is supported by a recent paper demonstrating, with a melanocyte cell line and traditional microarray analysis, that Gpnmb is coregulated with Dct and Si. ³⁵ Other categories with which Gpnmb also correlated in our study included genes involved in cell motility and migration. This finding is supported by a recent publication in which an association between Gpnmb and cell adhesion, a component of migration, was documented. ³⁶ Collectively, these data indicate that our systematic genetics approach is powerful in predicting physiological functions of a gene based on other transcripts similarly regulated at the mRNA level. 

Novel biological process categories represented in the GO enrichment analysis from the mutant Gpnmb data set are involved in posttranslational modifications and stress activated signaling, each of which includes kinases such as tyrosine protein kinase, serine/threonine kinase, receptor-associated kinase, and guanine nucleotide exchange. Another enriched category of the mutant Gpnmb correlate graph relates to visual perception. Seven of the eight genes in this category were positively correlated with Gpnmb expression, which meant these seven genes were downregulated in the mutant. Four of the downregulated gene products—Kcnv2, Opn1sw, Rho, and Cnga1—function as part of the phototransduction cascade, whereas the other three—Rom1, Opn1sw, and Rho—function as structural proteins. Several very careful studies in the D2 mouse have demonstrated age-related changes in visual acuity, visual function, and photoreceptor morphology, including inner and outer segment length and number of cells in the outer nuclear layer, ^37,38 though the structure of the outer segment is not compromised. ³⁹ Although ROM1 is a structural protein, elegant studies have demonstrated that this protein is involved in determining the fine structure of the outer segment as opposed to its binding partner, RDS/peripherin, which mediates folding of the outer segment membranous discs. ⁴⁰ Our data lend support for these functional and anatomic studies of the outer retina in the D2 mouse. They also further document the sensitivity of a systematic genetics approach; we are able to predict structural and functional abnormalities in D2 mice before they become manifest. 

Comparison of the outcomes of our present investigation with those of other investigators using microarray analysis of gene expression and mutations in Gpnmb and Tyrp1 yielded differences among all studies. For example, using iris samples taken from B6.Tyrp1^b Gpnmb^R150X B6 mice, Anderson et al. ⁴¹ identified an overexpression of several transcripts, including several crystallins (Gja3 and B3gnt5) compared with B6 mice. In contrast, Vetter et al. ⁴² found an upregulation of transcripts involved in the immune response along with a downregulation of crystallins in retinas from D2 mice at an advanced age. All differences can be readily explained in light of the differences in tissues used in each study and the genetic backgrounds of the mice. Nonetheless, each study provides valuable information. 

Given that mutations in both Gpnmb and Tyrp1 contribute to the glaucoma phenotype of D2 mice, it was unexpected that the expression level of Tyrp1 did not correlate with the expression level of mutant Gpnmb. Our unpublished data demonstrate that in BXD strains, the iris transillumination defect maps to Tyrp1 as early as 1 to 2 months. In contrast, the defect does not map to Gpnmb until mice are older than 13 months old. These data collectively suggest that although mutations in both genes increase the severity of the iris defect in PDS, they are differentially regulated and function within distinct pathways in the eye. 

To verify the ability of GeneNetwork to identify interactions of a gene of interest with the immune system, if such an interaction took place, we performed identical correlation analysis with Cd53, a gene expressed by T cells and natural killer cells. In RI lines with both wild-type and mutant Gpnmb, Cd53 was highly correlated with many members of various immunologic biological functions. The biological function categories enriched for the immune system included regulation of immune system processing, activation of immune response, regulation of acute inflammatory response to antigenic stimulus, and regulation of type I hypersensitivity. Moreover, the strength of the correlations was very high. 

In summary, our study demonstrates that Gpnmb is regulated by a very large cis-eQTL and that SNPs within Gpnmb itself, or another gene located very near Gpnmb on chromosome 6, are responsible for the modulation of its expression and that the truncated transcript attributed to the R150X mutation is not expressed at a significant level and likely undergoes nonsense-mediated decay. Genetic correlates of wild-type Gpnmb are involved in melanin synthesis and cell migration; supportive evidence of both functions is found in the literature. The R150X mutation disrupts the functional networks with which Gpnmb is normally associated. The new networks include stress activation, phosphorylation, and visual perception. Within a complex tissue such as the eye or the retina, neither wild-type nor mutant Gpnmb is highly correlated with genetic networks involved with the immune response. It is possible that an enriched population of ocular dendritic cells or antigen-presenting cells is required to reveal this connection. 

Supplementary Materials

Text s01, DOC - Text s01, DOC 

Text s02, DOC - Text s02, DOC 

Text s03, DOC - Text s03, DOC 

Text s04, DOC - Text s04, DOC 

Text s05, DOC - Text s05, DOC 

Text s06, DOC - Text s06, DOC 

Text s07, DOC - Text s07, DOC 

Footnotes

 Supported by National Institutes of Health Grants EY017841, EY021200, AA014425, AA017590, and DA021131; National Eye Institute Core Grant P30EY013080; Natural Science Foundation of Jiangsu China Grant BK2008186; UT-ORNL Governor's Chair in Computational Genomics; and an unrestricted grant from Research to Prevent Blindness.

Footnotes

 Disclosure: H. Lu, None; X. Wang, None; M. Pullen, None; H. Guan, None; H. Chen, None; S. Sahu, None; B. Zhang, None; H. Chen, None; R.W. Williams, None; E.E. Geisert, None; L. Lu, None; M.M. Jablonski, None

The authors thank Liyuan Li and Jesse Ingles for their assistance in generating the Western blot analysis.