Genomic locus modulating corneal thickness in the mouse identifies POU6F2 as a potential risk of developing glaucoma

Central corneal thickness (CCT) is one of the most heritable ocular traits and it is also a phenotypic risk factor for primary open angle glaucoma (POAG). The present study uses the BXD Recombinant Inbred (RI) strains to identify novel quantitative trait loci (QTLs) modulating CCT in the mouse with the potential of identifying a molecular link between CCT and risk of developing POAG. The BXD RI strain set was used to define mammalian genomic loci modulating CCT, with a total of 818 corneas measured from 61 BXD RI strains (between 60–100 days of age). The mice were anesthetized and the eyes were positioned in front of the lens of the Phoenix Micron IV Image-Guided OCT system or the Bioptigen OCT system. CCT data for each strain was averaged and used to QTLs modulating this phenotype using the bioinformatics tools on GeneNetwork (www.genenetwork.org). The candidate genes and genomic loci identified in the mouse were then directly compared with the summary data from a human POAG genome wide association study (NEIGHBORHOOD) to determine if any genomic elements modulating mouse CCT are also risk factors for POAG.This analysis revealed one significant QTL on Chr 13 and a suggestive QTL on Chr 7. The significant locus on Chr 13 (13 to 19 Mb) was examined further to define candidate genes modulating this eye phenotype. For the Chr 13 QTL in the mouse, only one gene in the region (Pou6f2) contained nonsynonymous SNPs. Of these five nonsynonymous SNPs in Pou6f2, two resulted in changes in the amino acid proline which could result in altered secondary structure affecting protein function. The 7 Mb region under the mouse Chr 13 peak distributes over 2 chromosomes in the human: Chr 1 and Chr 7. These genomic loci were examined in the NEIGHBORHOOD database to determine if they are potential risk factors for human glaucoma identified using meta-data from human GWAS. The top 50 hits all resided within one gene (POU6F2), with the highest significance level of p = 10−6 for SNP rs76319873. POU6F2 is found in retinal ganglion cells and in corneal limbal stem cells. To test the effect of POU6F2 on CCT we examined the corneas of a Pou6f2-null mice and the corneas were thinner than those of wild-type littermates. In addition, these POU6F2 RGCs die early in the DBA/2J model of glaucoma than most RGCs. Using a mouse genetic reference panel, we identified a transcription factor, Pou6f2, that modulates CCT in the mouse. POU6F2 is also found in a subset of retinal ganglion cells and these RGCs are sensitive to injury.


Introduction
Since the early Ocular Hypertension Treatment Studies (OHTS) [1] and subsequent independent findings of others [2,3], central corneal thickness (CCT) was identified as a factor related to the risk for developing primary open angle glaucoma (POAG). In these studies, CCT was a powerful predictor for developing POAG, with thinner corneas being associated with an increased risk of developing POAG [1-3] and this risk was independent of the confounding effects of CCT on intraocular pressure measurements [1,3]. The thinner CCT is also associated with an increased severity of visual field loss and a more rapid progression of the disease [4][5][6]. Furthermore, ethnic differences in CCT are correlated with increased risk of developing POAG and an increased severity of the disease [7]. Thus, there is a profound link between CCT and the risk of developing POAG.
In humans, there is a considerable variation in CCT, ranging from under 450 μm to over 650 μm [8][9][10] with a mean CCT of approximately 550 μm [8][9][10][11][12][13]. The variation in CCT is a highly heritable ocular trait, the genomic contribution of CCT is estimated to be near 90% [8,9,11,14]. Two studies on monozygotic and dizygotic twins from two different human populations confirmed the heritability, with the Chinese population having a heritability of 0.88 [8] [9] and the population in Australia and the United Kingdom having a heritability of 0.95 [8]. In addition, there is considerable variation in CCT between different ethnic groups. Aghaian et al. [15] found that African-Americans (mean CCT of 524 μm) and Japanese (mean CCT of 538 μm) had significantly thinner corneas than the general population. Others [10] have not observed the same differences for Japanese populations; however; thinner corneas are consistently observed in the African-American populations [10,13].
Genome-wide association studies (GWAS) on different human populations have identified a number of human loci/genes [16][17][18] associated with CCT. Several of the loci contain genes associated with CCT that are risk factors for human diseases, including ZNF496, which is associated with brittle cornea syndrome [19], and COL8A2, causing Fuch's endothelial corneal dystrophy [20,21]. A recent meta-analysis of over 20,000 individuals identified 16 additional loci associated with CCT [22]. Six of these loci conferred a significant risk for keratoconus, a disease characterized by an extremely thin cornea. This study of a large population successfully identified many different genes that contribute to the heritability of CCT and not surprisingly implicated collagens and extracellular matrix pathways in the regulation of CCT. Although this represents a significant increase in the number of genes involved in CCT, all of the loci only account for approximately 8% of the additive variance of CCT in the European population [22]. Taken together, these data reveal that this highly heritable trait, CCT, is a complex trait. In fact, it could be so complex that the contribution of many individual genomic elements would be difficult to prove, especially in a genetically heterogeneous species like humans.
Attempting to define the link between CCT and POAG in a human GWAS is complicated by the fact that one is comparing two complex traits. In this case, the effect size has to be very large. Many of the early studies specifically looked for a genetic/molecular association between CCT and POAG. The large GWAS meta-analysis of human CCT [22] identified many new loci associated with CCT. One of the loci that conferred a significant risk for keratoconus (FNDC3B), was also associated with POAG. Given the power of the study (20,0000 individuals) and the small effect of individual SNPs on CCT, it is clear that the association of individual SNPs with CCT and POAG will require an extremely larger (potentially unrealistic) sample size. Leveraging the unique genotype of recombinant inbred mouse strains can help to simplify such analysis.
The approach taken in the present study was to identify genomic loci modulating CCT using the largest recombinant inbred mouse strain set, the BXD strain set. The BXD RI strains were produced by crossing the C57BL/6J mouse with the DBA/2J mouse. The progeny were inbred (brother-sister matings) to produce over 80 inbred strains. The first 42 of these strains are from the Taylor series of BXD strains generated at the Jackson Laboratory by Benjamin Taylor [23]. BXD43 and higher were bred by the Williams group at the University of Tennessee (Peirce et al. 2004). The BXD RI strains offer a powerful tool to accurately identify genomic loci modulating phenotypes such as CCT. Among the strain set, there are over 7,000 breakpoints in the genome. All of the strains are fully mapped and the parental strains are fully sequenced. Thus, SNPs as well as insertions or deletions are known for every one of the BXD strains. With the aid of the bioinformatic tools present on GeneNetwork (Genenetwork.org), we were able to identify genomic loci and candidate genes modulating CCT in the mouse. We then used these data to interrogate human GWAS studies of corneal datasets [22] and glaucoma datasets [24,25] to relate our findings in the mouse to the normal human cornea and disease states [26]. Ultimately, we examine genes that modify CCT in the mouse that may also be risk factors for human glaucoma.

Results
The BXD RI strains were used to define genomic loci that modulate CCT. The CCT was measured in a total of 818 mice across 61 members of the BXD strain set (Fig 1). For the purpose of the present study we examined the total thickness of the cornea and no attempt was made to measure the separate layers of the cornea (epithelium, stroma or endothelium). The mean and standard error for each strain are shown in Fig 2. The mean central corneal thickness measured across 61 BXD strains was 100.9 μm with a standard deviation of 7.4 μm. The strain with the thinnest cornea was BXD 44 with an average corneal thickness of 85.3μm. The strain with the thickest cornea was BXD 74 with an average thickness of 124.6μm. The CCT for the parental strains is an intermediate value with the DBA2/J (D2) mouse having a CCT of 103.0μm and the C57BL/6J (B6) mouse having an average CCT of 93.1μm. Thus, there is considerable genetic transgression of CCT across the BXD RI strains with some BXD strains having corneas thinner than the parental strains and other strains having corneas thicker than the parental strains. This distribution of the CCT phenotype indicates that CCT is a complex trait, with multiple genomic loci segregating across the BXD RI strain set.
These data also demonstrate that CCT is a heritable trait. Fig 2 shows that there is considerable variability in the CCT from strain to strain and the standard error for each strain is rather small, which suggests that the genetic variability has a greater effect than the environmental variability. These data can be used to calculate the heritability of CCT in the BXD RI strains. Heritability (H 2 ) is the genetic variance (Vg) of the trait divided by the sum of genetic variance plus the environmental variance (Vg +Ve). The genetic variance can be estimated by calculating the standard deviation of the mean of the CCT for each strain (Vg = 10.054). The environmental variance can be estimated by taking the mean of the standard deviation across the   Using an unbiased forward genetic approach and the data collected from the cohort of 61 BXD strains, we performed a genome-wide scan in an effort to identify QTLs that modulate CCT. The genome-wide interval map (Fig 3) reveals one significant peak on chromosome 13 (13 to 19 Mb) and a suggestive peak on chromosome 7 (42 to 57 Mb). An expanded view of the peak on chromosome 13 is illustrated in Fig 4 along with a haplotype map of the BXD strains. The strains with the thicker corneas in general tend to have the D2 allele; while, strains with thinner corneas have the B6 allele. This is reflected in the genome wide map by the presence of the green line in the peak on chromosome 13 that indicates that higher phenotypic values are associated with the D2 allele (Fig 3). The QTL on Chr. 13 (Fig 4) rises above the significance level of p<0.05 as indicated by the pink line. The significant portion of the peak extends from 13Mb to 19Mb over Chr. 13. Genomic elements modulating CCT are located in this region. To identify candidates for modulating CCT in the BXD RI strains, we examined this 7Mb long region. The candidate genes can either be genomic elements with cis-QTLs or they can be genes with nonsynonymous SNPs changing protein sequence. Within this region are 29 traditional genes and one microRNA. In the Whole Eye Database (Eye M430V2 (Sep08) RMA) hosted on GeneNetwork.org we found two genes within this locus with cis-QTLs: Cdk13 and Mplkip. We examined the expression of both of these genes in the cornea and the eye using quantitative PCR (S1 Data, S1 Appendix). Neither of the two genes were expressed at a high enough level in the cornea to be monitored by microarrays of the whole eye. Thus, both genes were discounted as potential candidates for modulating CCT. Within the significant QTL on Chr 13, only one gene (Pou6f2) had nonsynonymous SNPs. Pou6f2 contained five separate nonsynonymous SNPs (rs29821949, rs52634762, wt37-13-18331131, rs29234524 and rs29250924). Three of these SNPs result in a change in the amino acid proline, which could cause a change in secondary structure of the protein and an alteration in protein function. On the y-axis is the linkage related score (LRS). Notice one significant QTL peak on Chr13 (above the pink Line, p = 0.05) and additional suggestive peaks (above the gray Line). https://doi.org/10.1371/journal.pgen.1007145.g003 Analyzing these changes using SIFT to predict functional changes in the protein structure [27], one SNP (rs29234524) was predicted to potentially alter the function of the protein.
Based on this examination of the locus, there are no valid cis-QTLs and only one gene, Pou6f2, with nonsynonymous SNPs that could have biological effects.
The genomic elements within significant QTLs on mouse chromosome 13 (13 to 19 Mb) were examined to determine if there were specific associations in two human datasets: the International Glaucoma Genetics Consortium dataset for association with CCT [22], and the  13. The haplotype map for the 61 strains in the CCT dataset is shown in the top panel. As indicated at the key on the right, red represents the B6 alleles, green defines the D2 alleles, blue represents regions of the DNA that are heterozygotic and gray is unmapped. The genomic markers used in the mapping process are listed at the bottom of the haplotype map. At the far right is a list of the specific BXD RI strains and the associated CCT measurements going from thick corneas at the top to thin corneas at the bottom. QTL map of CCT on Chr. 13 is shown below the haplotype map. An LRS of above the pink line is statistically significant (p<0.05). A positive additive coefficient (green line) indicates that D2 alleles are associated with higher trait values. The light fill under the peak defines the genomic locus associated with the modulation of CCT. The significant peak for CCT is from 13 to 19 Mb (gray shaded area). Vertical orange lines at the bottom of the plot show the SNPs on Chr. 13.
https://doi.org/10.1371/journal.pgen.1007145.g004 NEIGHBORHOOD glaucoma dataset [28]. The syntenic regions on human chromosome 1 (235 to 238 Mb) and chromosome 7 (38 to 43 Mb) were queried for associations of genetic effects relative to human CCT and with genetic risk for human glaucoma. In the International Glaucoma Genetics Consortium dataset for CCT [22], there were a number of nominally significant associations; however, none of the markers remained significant after corrections for multiple testing. When we looked specifically for the candidate genes from the mouse, weak associations were found in both European populations (peak marker rs4723833; P = 4.34x10 -3 β = -0.062μm) and Asian populations (peak marker rs17619647; P = 3.72x10 -3 ; β = -0.066μm), both of these SNPs (rs4723833 and rs17619647) overlap POU6F2. However, the associations are not significant after correction for multiple testing of 482 SNPs in Europeans and 265 SNPs in Asians. The β values from this linear regression analysis reflect the quantitative change in CCT per associated allele. The mouse locus was also used to interrogate the NEIGHBOR-HOOD human POAG glaucoma dataset [28]. In the NEIGHBORHOOD dataset, the top 50 SNPs associated with glaucoma within the syntenic region were all found within the POU6F2 locus. The most significant SNP (rs76319873) had a genome-wide p-value of 5.34 −6 and the next top 3 SNPs had p-values less than 10 −5 ( Fig 5, Table 1). None of these SNPs reached genome-wide significance when corrected for multiple testing.

Distribution of POU6F2 in the mouse eye
Our data suggests that Pou6f2 modulates central corneal thickness and it may also modulate risk for glaucoma in humans. To determine if Pou6f2 could be a molecular link between CCT and glaucoma, we examined the distribution of POU6F2 protein and mRNA in the retina and cornea in mice. The first approach was to examine its absolute expression levels in adult retina and cornea using digital droplet PCR. We first confirmed the specificity of tissue dissection using the corneal marker Uroplakin 1 [29], which was highly abundant in cornea but not retina (5760±260 copies/μL vs. 3.3±0.9 copies/μL). The expression levels of Pou6f2 were relatively high in the retina (76.3±1.5 copies/μL); however, levels of Pou6f2 message from cornea were just above the detection threshold of digital droplet PCR (0.5±0.2 copies/μL), suggesting that only few if any corneal cells express Pou6f2. These data indicate that Pou6f2 was barely expressed in the adult cornea, while it was readily detectable in the adult retina.
Previous studies [30] demonstrated that POU6F2 is expressed in neuroblasts in the future ganglion cell layer of the developing eye and in retinal ganglion cells in the adult mouse, cat and monkey. To confirm these findings and to examine the potential for a link between RGCs and cornea, we immunostained sections from the embryonic mouse eye, the adult mouse eye and flat-mounts of mouse retina. In the embryonic eye, there is strong staining of the neuroblasts destined to become retinal ganglion cells (Fig 6). There is also notable staining of the epithelium of the developing cornea and what appears to be corneal stem cells (Fig 6). The distribution of POU6F2 in the developing eye indicates a clear association between the development of retinal ganglion cells and the cornea.
In sections through the adult mouse retina (C57BL/6J), POU6F2 is found to label the nuclei of many cells in the retinal ganglion cell layer. To fully explore the distribution of POU6F2 in the retinal ganglion cell layer, we examined flat mounts of the mouse retina (Fig 7) counterstained with RBPMS (a ganglion cell marker [31]) and a nuclear stain. We observed that POU6F2 labels only a small percentage of the total number of cells in the ganglion cell layer. When we quantified the number of cells labeled with POU6F2, only 17.4% of the total cells were heavily labeled with POU6F2. Virtually all of the POU6F2-positive cells were ganglion cells positive for RBPMS. There were a few cells in the inner nuclear layer that were positive for POU6F2 and these cells were also positive for RBPMS, identifying these cells as displaced ganglion cells. These data demonstrate that POU6F2 is expressed in a small subset of retinal ganglion cells. To provide additional evidence that the POU6F2 positive cells are ganglion cells, we crushed the optic nerve unilaterally in three mice and allowed them to survive for 28 days. The retinas were then examined for POU6F2 staining. No labeling was observed in the RGC layer, indicating that all POU6F2 cells were gone (Fig 1, S1 Appendix). Taken together,   these data reveal that POU6F2 is expressed in a subset of retinal ganglion cells in the ganglion cell layer. Finally, as previously observed by Zhou et al. [30] there were a few cells labeled with POU6F2 on the inner surface of the inner nuclear layer. These cells appeared to be displaced ganglion cells (Fig 2, S1 Appendix). Thus, POU62 labels a subset of RGCs within the mouse retina.
To determine if the POU6F2-positive RGCs are selectively sensitive to glaucoma, we double stained 4 young DBA/2J (2 months old) mice for POU6F2 and RBPMS and compared these values to 4 aged (8-month-old) DBA/2J mice (Fig 8). When we counted the cells in the 4 young mice there were an average of 460 (±81, SEM) RBPMS-labeled RGCs per 40X field. In the same sections, we observed an average of 74 POU6F2-labeled RGCs (±15). When these results were compared to 4 aged DBA/2J mice (8 months old) there was a significant decrease (p<0.001, student t test) in the total number of RGCs to an average of 361 RGCs (±21) per 40 X field. This represented a total decrease in RBPMS labeled ganglion cells of 22%. When we examined POU6F2-positive cells in the aged DBA/2J retinas there was a 73% loss of cells, a significant difference relative to young DBA/2J retinas (p<0.0001). The average number of POU6F2 cells in the young retina was 74 (±5), while in the aged retina this number decreased to 21 (±3). These data demonstrate that the POU6F2 RGC subtype is very sensitive to early phases of glaucoma in the DBA/2J model.
To examine the expression of POU6F2 in the cornea, we first examined sections through the cornea and limbus extending down into the sclera. There was no staining in the cornea or sclera. Occasionally, we observed individual cells stained in the limbal area. However, it was difficult to unequivocally identify these cells as limbal stem cells. As an alternative approach, we stained the surface of six eyes using extended antibody incubation times. When examining the junction between the cornea and sclera, a line of cells is stained at the junction between the cornea and sclera in the limbus (Fig 9). This line of cells was also positive for the stem cell There was a 22% loss of RBPMS-labeled RGCs in aged DBA/2J mice (8 months of age, Aged D2) as compared to young DBA/2J mice (2 months of age, Young D2). There was a dramatic loss of 73% of the POU6F2-positive cells comparing the Young D2 mice to the Aged D2 mice. These data demonstrate the sensitivity of the POU6F2 RGC subtype to glaucoma.
https://doi.org/10.1371/journal.pgen.1007145.g008 marker ABCB5 (Fig 9E). Thus, POU6F2 is not only a marker for limbal stem cells but it may be responsible in part for the maintenance of the corneal integrity in the adult mouse eye.
To determine if POU6F2 is directly involved in corneal development and maintenance affecting CCT, we examined the corneas of 8 Pou6f2-null mice and compared their CCT to that from 10 wild-type littermates (Fig 10). There was a significant difference (p<0.01, Wilcoxon signed-rank test) in CCT between the Pou6f2-null mice and the wild-type littermates. In the Pou6f2-null mice the mean CCT was 102.8μm (±1.8μm). The wild-type littermates had a mean CCT of 109.9μm (±1.8μm). Thus, abolishing Pou6f2 expression has a significant effect on CCT in the adult mouse.

Discussion
The BXD RI strain set is particularly well suited for a systems genetic approach to study CCT. This relatively new branch of quantitative genetics has the goal of understanding networks of interactions across multiple levels that link DNA variation to phenotype [32]. The present approach involves an analysis of sets of causal interactions among classic traits such as CCT, together with networks of gene variants, and developmental/environmental factors. The main challenge is the need for a comparatively large sample size and the use of more advanced statistical and computational methods and models. The BXD strain set is sufficiently large to have adequate power for this approach [33,34]. The genetic reference panel used in this study consisted of 62 strains of mice. The novel aspect of our current approach is the combination of phenotypic mouse data, mRNA expression data (baseline and experimental) from the same BXD RI strains and the interrogation of GWAS studies for glaucoma (NEIGHBORHOOD Study [24]) and for CCT [22]. Data that we generated throughout this experiment and data already deposited in GeneNetwork (www.genenetwork.org) for gene expression in the eye [33] and retina [35][36][37] were used to test specific mechanisms and predictions.
In the present study, we have identified a significant QTL on Chr 13 (13)(14)(15)(16)(17)(18)(19) Mb) that modulates CCT in the BXD RI strain set. One suggestive QTL was also found on chromosome 7 (41 to 57 Mb). Previous studies from the Anderson group [38,39] have found significant QTLs associated with CCT. Using a F2 cross between C57BLKS/J and SL/J F2 mice, Lively et al. [38] found a significant QTL on Chr. 7 at 105Mb that modulated CCT. This differs from the suggestive QTL identified on Chr. 7 in the present study. The QTL identified with the C57BLKS/J and SL/J F2 cross is located at 105Mb on Chr7, while the suggestive QTL from the present study was at 41 to 57Mb on Chr7. A second study by Koehn et al. [39], identified two significant QTLs associated with CCT in an F2 cross between BXD24/TyJ and CAST/EiJ. One was on Chr 3 at 104Mb, and the second was on Chr11 at 88Mb [39]. Neither of these loci demonstrated even suggestive associations with CCT in the BXD RI strains used in the present study.
In humans, CCT is a highly heritable trait affected by many genomic loci [16][17][18][19][20][21][22]. Comparing these human loci to those identified in the mouse reveals several genomic regions that are similar. For example, the mouse locus on Chr7 at 105Mb [38] could include previously identified human loci near several known genes: TJP1, CHSY1, LRRK1 and AKAP13 [22,40]. Other than this region in the mouse, there does not appear to be any overlap between loci modulating CCT in the mouse and those identified in the human.
CCT is a phenotypic risk factor for glaucoma [1-3]. Glaucoma is a diverse set of diseases with heterogeneous phenotypic presentations associated with different risk factors. Untreated, glaucoma leads to permanent damage of axons in the optic nerve and visual field loss. Millions of people worldwide are affected [41,42] and glaucoma is the second leading cause of blindness in the United States [43]. Adult-onset glaucoma is a complex collection of diseases. The severity of glaucoma appears to be dependent on the interaction of multiple gene variants, age, and environmental factors [44]. These complex genetic risk factors are the focus of recent genome-wide association studies that are not only defining genetic risk factors but also aiding in our understanding of disease mechanisms [45][46][47]. These would include the identification of genes that could influence cerebral spinal fluid pressure [48] and ultimately the differential pressure observed across the optic nerve head [49]. Recent findings suggest that RGC number is also a risk factor for glaucoma. Several polymorphisms in the genes SIX6 and ATOH7 were identified that are responsible for thinner retinal nerve fiber layer due to fewer RGCs, suggesting that humans with fewer RGCs are at increased risk of developing glaucoma [25,50,51]. In addition, linkage analysis in families affected by glaucoma has led to an extensive list of genomic loci linked to POAG. The genes underlying some of the adult-onset glaucoma loci have been identified (MYOC, OPTN, TBK1, CDKN2BAS [52] [53] [54,55].
Many studies in human and mouse have looked for a genetic link between CCT and glaucoma. Studies in the human cornea have found several SNPs that modulate corneal thickness. However, none of the genes in these studies appear to relate CCT to the risk of developing glaucoma with the exception of FNDC3B [22]. Several studies have examined the genomic loci in the mouse that modulate CCT and again none of the loci controlling CCT thickness relate to glaucoma risk. In the present study, we have identified a well-defined QTL associated with CCT in the mouse. Within this locus, there is a limited number of candidate genes. One gene, Pou6f2, is uniquely poised to modulate CCT, for we have shown that it is expressed in the developing cornea and in limbal stem cells. In addition, Pou6f2 may also be associated with glaucoma risk. Interestingly, Zhou et al. [30] demonstrated that Pou6f2 is present in the developing and adult retinal ganglion cells. This gene product is also associated with the regeneration of neurons in zebrafish [56].
Based on the distribution of POU6F2 in the cornea and retinal ganglion cells in the developing and adult mouse eye, it is tempting to suggest that Pou6f2 may be a molecular link between CCT and the potential for RGC loss in glaucoma. A search of the normal mouse retina (DoD Normal Retina Affy MoGene 2.0 ST [36]) on GeneNetwork (genenetwork.org) reveals that Pou6f2 forms a genetic network with retinal ganglion cells markers, including Thy1, Pou4f1 (Brn3a) and Rbfox3 (NeuN) [57]. This association suggests that Pou6f2 could be commonly regulated in a specific subset of retinal ganglion cells. These RGC subtypes may be selectively sensitive to loss in mouse models of glaucoma. In the present study, we have also shown that Pou6f2 is expressed in the developing cornea and corneal stem cells. When we examined the NEIGHBORHOOD database, POU6F2 has an interesting association with POAG. However, for the International Glaucoma Genetics Consortium dataset for CCT [22] the P value was p = 0.0037, and when corrected for multiple testing it did not reach a significant level. Thus, we are not able to definitively say that Pou6F2 is a molecular link between CCT and glaucoma. As discussed above, CCT is a highly heritable trait and it is the second highest risk factor for glaucoma, following IOP. If this is the case, then why do we not know the molecular link between CCT and POAG? Some phenotypic traits are very complex or even hyper-complex. A well-known example of a complex trait is height. It is estimated that there are over 2,000 genomic loci that may contribute to height and that the heritability is relatively high as much as 80% [58], although others estimate that the heritability is lower due to dependence on family-based studies indicating that the true heritability is between 60% to 70% [59]. Nonetheless, the weight of any one genomic locus influencing height in the human population is low and genetic components known to be associated with heritability can only be accounted for a portion of the heritability. That being said, CCT is a glaucoma risk factor and a highly heritable trait. The identification of molecular links connecting CCT to POAG may be complicated by the fact that both are complex traits, influenced by multiple genomic as well as environmental factors. Using a mouse model like the BXD RI strain set (a restricted genetic reference panel) is an effective approach to identify genomic elements modulating phenotypes such as CCT. These data can now be used to design experiments that will attempt to link CCT to potential glaucoma risk.
POU6F2 appears to mark a subpopulation of RGCs that are selectively sensitive to injury in the DBA/2J mouse model of glaucoma. To understand the potential role of POU6F2 we examined data from a ChIP-chip experiment [60] in HEK293 cells. In this human cell line, POU6F2 bound to a number of different targets: CHD5, DACH1, ELAVL3, GFRA2, GRIN1, LHX1, MTSS1, NAP1L3, NEFH, NPTX1, NR4A2, NTNG2, PCDH7, ROBO2, SLC30A3, SLC7A8, SORCS2, UNC5A and VGF. If we examine these targets in the DoD Normal Retina Database [36] on GeneNetwork, all of the targets are expressed at relatively high levels in the retina, at least 2-fold above the mean expression of mRNA. In addition, their expression levels are highly correlated across the BXD strain set, with most probes having a Pearson's r value above 0.7. One of these down-stream targets, Slc30a3, encodes a zinc transporter (ZNT3). Recent studies from the Benowitz lab [61], reveal that ZNT3 plays an important role in injured ganglion cells and axon regeneration. Thus, it is tempting to hypothesize that POU6F2 exerts its effects on retinal ganglion cell survival in an injury and chronic neurodegeneration (glaucoma) context by modulation of Slc30a3.
In conclusion, we have shown that POU6F2 is involved in corneal development and modulates CCT. In retinal ganglion cells, POU6F2 is a transcription factor marking a subtype that is selectively sensitive to injury. POU6F2 is also known to be upstream of genes that play a critical role in ganglion cell death following injury and it is a potential glaucoma risk factor in humans.

Ethical statement
All of the procedures involving mice were approved by IACUC at Emory University and the University of Tennessee Health Science Center. The study adhered to the ARVO Statement for the Use of Animals in Research. Mice were anesthetized via intraperitoneal injection of ketamine 100mg/kg and xylazine 15mg/kg. For euthanasia, mice were anesthetized via intraperitoneal injection of ketamine 100mg/kg and xylazine 15mg/kg, and perfused through the heart with saline followed by 4% paraformaldehyde in phosphate buffer.

Mice
All of the mice in this study were between 60 to 110 days of age, which is long before any significant elevation in IOP due to two gene mutations (Tyrp1 and Gpnmb) carried by selected BXD strains originating from the DBA/2J mouse. We examined equivalent numbers of males and females for each of the strains. Mice were housed in a pathogen-free facility at UTHSC or at Emory University, maintained on a 12 hr light/dark cycle, and provided with food and water ad libitum. The B6.129-Pou6f2 tm1Nat /J mice were cryorecovered at Jackson Laboratory (Bar Harbor, ME) and bred at Emory. These mice were originally created in Jeremy Nathans' lab by replacing an approx. 3.5 kb long fragment of Pou6f2 that encompasses exon and flanking intronic regions with a PGK-Neo positive selection cassette in 129 ES cells. Animals were back-crossed to the C57BL/6J strain for more than 10 generations. All mice were genotyped in-house and only homozygote null or wild-type littermates were used in the present study. All CCT measurements on the Pou6f2-null and wild-type mice were performed in a double-blind manner. The mice were genotyped after all CCT measurements were made. To our knowledge, this is the first publication using the mice generated by Dr. Nathans.

Measuring CCT
CCT was measured using two different Ocular Coherence Tomography (OCT) systems: a Bioptigen SD-OCT system (Morrisville, NC) and a Phoenix Micron IV OCT Imaging system (Pleasanton, CA). Mice were anesthetized using ketamine 100mg/kg and xylazine 15mg/kg. The eye was positioned in front of the lens. The entire anterior chamber was imaged. The corneal scans were saved to a portable hard drive for subsequent analysis. CCT was then measured three times for each eye using the Mouse Retina Program, InVivoVue Clinic, in the Bioptigen Software or the Micron OCT software. These data were stored and entered into an Excel spreadsheet. Repetitive measures were possible with the BXD RI strains since each strain has the identical genetic background allowing us to sample the developmental consequences of the same genetic background for all independent measures. For the BXD strains the corneas were measured between 60 and 90 days. In the Pou6f2-null mice the corneal thickness was measured at 30 days of age.

Optic nerve crush
Optic nerve crush was performed as described in Templeton and Geisert [62]. Briefly, three C57BL/6J mice were anesthetized using ketamine 100mg/kg and xylazine 15mg/kg. Under the binocular operating scope a small incision was made in the conjunctiva. With micro-forceps (Dumont #5/45 Forceps, Roboz, cat. #RS-5005, Gaithersburg, MD), the edge of the conjunctiva was grasped next to the globe. The globe was rotated nasally to allow visualization of the posterior aspect of the globe and optic nerve. The exposed optic nerve is then clamped 2 mm from the optic nerve head with Dumont #N7 self-closing forceps (Roboz, cat. #RS-5027) for 10 seconds. At the end of the procedure, a drop of 0.5% proparacaine hydrochloride ophthalmic solution (Falcon Pharmaceuticals, Fort Worth TX) was administered for pain control and a small amount of surgical lube (KY Jelly, McNeil-PPC, Skillman NJ) was applied to the eye to protect it from drying. The mouse was allowed to wake up on a heating pad and monitored until fully recovered.

Interval mapping for the traditional phenotypes
CCT data was subjected to conventional QTL analysis using simple and composite interval mapping and pair-scans for epistatic interactions. Genotypes were regressed against each trait using the Haley-Knott equations implemented in the WebQTL module of GeneNetwork [63,64]. Empirical significance thresholds of linkage were determined by permutations [65]. We correlated phenotypes with expression data for whole eye and retina generated by Geisert, Lu and colleagues [33,36]. Based on recent work with the enlarged set of BXDs, we expected all significant (p < 0.05) QTLs to be mapped with a precision of ± 2 Mb [34,[66][67][68]. To identify loci, and also to nominate candidate genes, we used the following approaches: interval mapping for the traditional phenotypes, candidate gene selection within the QTL region, cis-eQTL analysis of gene expression, trans-eQTL analysis, multi-trait and complex analysis of molecular, clinical, and laboratory traits.

NEIGHBORHOOD analysis
The peak area of association on mouse chromosome 13 was examined for associations in human datasets. Syntenic regions on human chromosomes 1 and 7 were queried for associations with POAG in the NEIGHBORHOOD [24] dataset. Subsequently, the POU6F2 locus was queried in the International Glaucoma Genetics Consortium dataset for association with CCT [22].

RNA isolation and digital PCR
In order to quantify Pou6f2 mRNA expression, C57BL/6J mice (n = 5) were deeply anesthetized with a mixture of 15 mg/kg of xylazine and 100 mg/kg of ketamine and sacrificed at 9am. Corneas and retinas were dissected separately and stored in Hank's Balanced Salt Solution and RiboLock (Thermo Scientific, Waltham MA) at -80˚C. RNA was isolated on a Qiacube with the RNeasy mini kit (both Qiagen, Hilden, Germany) according to the manufacturer's instructions with additional on-column DNase1 treatment. RNA integrity was assessed using an Agilent Bioanalyzer 2100 and RIN scores for both pooled tissues were >9.5. Takara PrimeScript (Clontech, Mountain View, CA) was used to retrotranscribe equal amounts of RNA for both tissues. Digital Droplet PCR was then carried out using 30ng of total RNA in 20μL reactions of EvaGreen ddPCR supermix supplying 2μM Mg 2+ . Primers were designed using NCBI Primer-Blast to work at a combined annealing/extension temperature step of 60˚C and the sequences were as follows: Upk1b fwd (5-CAGGCAGCCGGTCTTTTAGAAA-3), Upk1b rev (5-ATCA TTGTTGGTGGCTTCGAGA-3), Pou6f2 fwd (5-CCCTCAATCAGCCAATCCTCAT-3), Pou6f2 rev (5-GTTCAGGGATGAGGTAGCTTGT-3). A combination of Ppia and Gapdh (Qiagen Quantitect primer assays) was used for ddPCR normalization.

Immunohistochemistry
Twelve C57BL/6J mice were deeply anesthetized with a mixture of 15 mg/kg of xylazine and 100 mg/kg of ketamine and perfused through the heart with saline followed by 4% paraformaldehyde in phosphate buffer (pH 7.3). The eyes were removed and embedded in paraffin and sectioned at 5μm. For the embryonic eyes, two timed-pregnant female mice were deeply anesthetized with tribromoethanol, decapitated and the E16 embryos were removed. The embryo heads were placed in 4% paraformaldehyde on ice for 1 hour, rinsed with PBS and dehydrated in stages of ethanol followed by xylenes, and then embedded in paraffin. The 5μm sections of retina and E16 mice had the paraffin removed and the sections were blocked with 5% normal donkey serum and stained with a rabbit antiserum directed against POU6F2 (MyBiosource, Cat. # MBS9402684) at 1:500 for 2 hours at room temperature (S2 Data, S1 Appendix). To demonstrate the specificity of the POU6F2 antibody, retinas from Pou6f2-null mice were stained and no nuclear labeling was observed. The sections were rinsed three times and transferred to secondary antibodies (Alexa Fluor 488 AffiniPure Donkey Anti-Rabbit, Jackson Immunoresearch, Cat. #715-545-152; at 1:1000 for two hours at room temperature. After three washes of 15 minutes each, To-PRO-3 Iodide (Molecular Probes, Cat. # T3605) was applied 1:1000 as a nuclear counterstain. After a final wash in PBS, coverslips were placed over the sections using Fluoromount -G (Southern Biotech, Cat. #0100-01) as a mounting medium. For the retinal flat mounts, the retinas were removed from the globe and rinsed in PBS, blocked in 5% normal donkey serum and placed in primary antibodies in primary rabbit antibody POU6F2 (MyBiosource, Cat. # MBS9402684) at 1:500 and to label RGCs, a primary guinea pig antiserum against RBPMS (Millipore, Cat # ABN1376) at 1:1000 overnight. The retinas were rinsed and placed in secondary antibodies, (Alexa Fluor 488 AffiniPure Donkey Anti-Rabbit, Jackson Immunoresearch, Cat. #715-545-152 and Alexa Fluor 594 AffiniPure Donkey Anti-Guinea pig, Jackson Immunoresearch, Cat. #706-585-148) at 1:1000 for two hours at room temperature. After three washes of 15 minutes each, To-PRO-3 Iodide (Molecular Probes, Cat. # T3605) was applied 1:1000 as a nuclear counterstain. After To-PRO3 has stained the section for 30 minutes, slides were washed twice with PBS for 10 minutes followed by a 5-minute wash with distilled water. Autofluorescence was decreased by treating the sections with cupric sulfate [69]. The 10mM cupric sulfate pH5 was applied to the section for 5 minutes. Slides were washed for 5 minutes with distilled water, then twice with PBS for 10 minutes and the slides were coverslipped.

Staining the surface of the limbus
Adult C57BL/6J mice were anesthetized with a mixture of 15 mg/kg of xylazine and 100 mg/kg of ketamine and perfused with saline followed by 4% paraformaldehyde in phosphate buffer (pH 7.4). The eyes were removed and post-fixed in 4% PFA for one hour at room temperature. The sclera along the ministry equator line was removed along with the posterior sclera and whole retina, vitreous body and lens. The remaining portion of the eye was rinsed three times in PBS containing 2% Triton X-100 for 10 minutes each. Dissected blocks the anterior hemisphere were placed in 10% donkey serum and 4% BSA overnight 4˚C. The tissue was incubated for two days in rabbit primary antiserum to POU6F2 (MBS9402684 MyBiosource, San Diego CA) and ABCB5 (Product ab126 Abcam, Cambridge MA) at 1:500 at 4˚C. The dissected limbal area was washed 3 times with PBS. Finally, the tissue was placed in secondary antibodies that included Donkey anti-Rabbit labeled with Alexa 488 (#711-545-152, Jackson ImmunoResearch Laboratories, West Gove CA) and Donkey anti-Goat IgG labeled with Alexa 594 (#705-585, Jackson ImmunoResearch Laboratories, West Gove CA) at 4˚C overnight. The tissue was cut into 5 pieces, radially mounted on a glass slide and a coverslip was placed over the tissue.

Members of the NEIGHBORHOOD consortium
Genome Research Institute (NIH), Baltimore, MD, USA. 57 Department of Ophthalmology, National University of Singapore and National University Health System, Singapore. 58 Department of Ophthalmology and Visual Sciences, School of Medicine and Public Health, University of Wisconsin, Madison, Wisconsin 53705, USA. 59 Clinic for General and Interventional Cardiology, University Heart Center Hamburg, Hamburg, Germany.
Supporting information S1 Data. Expression of Mplkip and Cdk13 in the cornea and the eye in the mouse. We examined the cornea and the eye cup for the expression of the transcripts monitored by three of the probes designated as cis-QTLs in the whole eye dataset. The expression level of Mplkip and Cdk13 is considerably higher (3-to 5-fold) in the eye cup relative to the cornea per ng of cDNA. We next measured the RNA content of the cornea in six independent biological samples relative to the eye cup in four independent biological samples. The cornea contained approximately 1.92% of the total RNA in the eye. To determine the contribution of the RNA expression in the cornea to the signal monitored by microarrays in the whole eye dataset, we multiplied the signal by the relative amount of RNA. The total contribution of corneal signal for Mplkip and Cdk13 to that monitored in the dataset is less than 1%. Thus, if the cis-QTL is real, then it does not represent a corneal QTL for there just is not enough expression in the cornea to be monitored in a whole eye sample ( Figure A). (DOCX) S2 Data. To validate the specificity of rabbit antiserum directed against POU6F2 (MyBiosource, Cat. # MBS9402684) we used in the present study, we first ran an immunoblot and found one major band in the retina sample at approximately 73 kDa, the appropriate size for POU6F2 (lane A). There were no significant bands observed in blots of similar tissue that were stained with secondary antibody only (lane B). Next, we examined four different tissues by PCR, including: retina, brain, colon and salivary gland. We found significant levels of Pou6f2 message in retina and brain, with virtually no Pou6F2 mRNA detected in the other samples. When we examined protein samples from these tissues using immunoblot methods, we observed a 74 kDa band in retina and brain and no bands in colon or salivary gland. This confirmed that the antibody was specific for POU6F2. (DOCX) S1 Fig. We conducted a series of RNA isolations and quantitative PCR to determine the relative expression of Mplkip and Cdk13 in the cornea and the eye. Based on ng of cDNA from the tissues, the expression within the cornea was less than 25% of the that in the eye. Almost no expression of the corneal marker Upk1b was found in the retina. The total amount of RNA from the cornea in a sample from the eye was only 1.92%. If we determine the relative contribution of corneal Mplkip and Cdk13 to the signal coming from whole eye samples it is less than 1% and thus negligible. The corneal expression of Mplkip and Cdk13 would not be reflected in the whole eye sample. For genes expressed at high levels in the cornea, this is not the case, as it can be seen that the majority of the signal from Upk1b in the whole eye sample originates form the cornea.