Abstract
Photoreceptor cells (PRCs) are the light-detecting cells of the retina. Such cells can be non-invasively imaged using optical coherence tomography (OCT) which is used in clinical settings to diagnose and monitor ocular diseases. Here we present the largest genome-wide association study of PRC morphology to date utilising quantitative phenotypes extracted from OCT images within the UK Biobank. We discovered 111 loci associated with the thickness of one or more of the PRC layers, many of which had prior associations to ocular phenotypes and pathologies. We further identified 10 genes associated with PRC thickness through gene burden testing using exome data. In both cases there was a significant enrichment for genes involved in rare eye pathologies, in particular retinitis pigmentosa. There was evidence for an interaction effect between common genetic variants, VSX2 involved in eye development and PRPH2 known to be involved in retinal dystrophies. We further identified a number of genetic variants with a differential effect across the macular spatial field. Our results suggest a continuum between normal and rare variation which impacts retinal structure, sometimes leading to disease.
Introduction
The retina is a layered structure at the back of the eye, responsible for receiving light stimulus and converting it to neurological signal which is interpreted by the visual system of the brain. Each retinal layer comprises different cells that are responsible for particular stages of the signal conversion process. The central area of the retina is termed the macula, and has a valley like structure. At the centre of the macula, the bottom of the valley structure, is the fovea which is the area of the retina responsible for highest acuity vision [1]. The retina can be imaged using optical coherence tomography (OCT). OCT produces high resolution, three-dimensional representations of the retina in which the different retinal layers and macular structure can be identified (fig. 1A). OCT is a non-invasive imaging method commonly used in the clinic to diagnose and monitor a number of ophthalmic conditions [2].
(A) An optical coherence tomography (OCT) image of the cross-section of the retina. The layers are segmented, outlined in blue, by the Topcon advanced boundary software (TABS). On the left-hand side of the image, the PRC layers have been shaded: The outer nuclear layer (ONL) in yellow, the inner segment (IS) in orange and the outer segment (OS). (B) The early treatment of diabetic retinopathy study (ETDRS) segmentation grid super imposed on a fundus image of the retina. The ETDRS grid is composed of nine segments arranged in a bulls-eye. (C) A schematic of the data utilised to describe PRC morphology as an input phenotype in genome-wide association studies (GWAS). The ONL, IS and OS, combined, comprise the photoreceptor layer (PRL).
The photoreceptor cells (PRCs) are found towards the back of the eye within the retina, backed by the retinal pigment epithelium layer (RPE) (fig. 1A), and are responsible for detecting light photons and generating electrical response through the process of phototransduction. They consist of two main types, the rods and cones. Rods are responsible for vision in low light but offer lower visual acuity and lack colour definition. Cones are active in higher light levels, and confer colour and high acuity vision. The cones are densely populated at the centre of the retina, the fovea, where there are no rods. Conversely, rods are most densely populated more peripherally, where cone denisty is low. The PRC layer can be further divided into three component intracellular layers; the outer nuclear layer (ONL), which comprises the cell body of the PRC, including the cell nucleus; the inner segment (IS), which contains both the PRC’s mitochondria and ribosomes required for photopigment assembly; and the outer segment (OS), which interfaces with the RPE and contains the stacked membranes in which photopigments are stored.
There are numerous pathologies, both rare and common, that affect the morphology of the PRCs. Notably, age-related macular degeneration (AMD) is defined by the accumulation of lipid deposits, known as drusen, under the RPE. This causes damage to the RPE and subsequently degradation of the PRCs [3]. AMD causes progressive and irreversible loss of visual function, and is a major cause of blindness in older people. [4].
Many retinal dystrophies affect the photoreceptors [5]. These rare inherited disorders are due to highly penetrant mutations in genes expressed in PRCs or less commonly the RPE, resulting in retinal dysfunction with or without degeneration of the outer retinal layers. These can affect predominantly the central retina (macular dystrophies), the rod photoreceptors (rod-cone dystrophy or retinitis pigmentosa), or the cone photoreceptors (cone and cone-rod dystrophies). Though individually rare, taken together these diseases are a major cause of blindness in children and working age adults.
There have been previous studies on retinal morphology with known associations to age, sex, and ethnicity [6, 7]. Previous studies have further looked at the genetic variation underlying the overall retinal thickness [8] and the thickness of the inner retinal layers [9]. These studies identified loci with prior associations to general retinal development and ophthalmic disease.
In this paper, we utilise data within the UK Biobank to interrogate how genetic variation affects PRC morphology using quantitative phenotypes extracted from OCT images. The UK Biobank provides a large-scale rich resource of both genotypic and phenotypic data. Here we conduct genome-wide association studies (GWAS) of the thickness of the component PRC layers, the ONL, IS and OS. Following meta analysis of the three PRC layers, 111 loci are significantly associated with the thickness of one or more PRC layers. These loci include variants with prior associations to ophthalmic phenotypes, including ophthalmic diseases. We further explored the genetic variation of the morphology of the retina at a higher resolution by considering the retinal layer thickness in different concentric areas of the macula. This identified a number of loci that appear to differentially affect PRC morphology across the retinal field. Interestingly, several of these loci have prior associations to retinal diseases whose pathological action bears similarity with the spatial variation in genetic effect on morphology. Further analysis utilising the exome sequence data identified 10 genes associated with the thickness of one or more of the PRC layers. Several of these genes also had prior associations to ocular phenotypes. Additional analysis identified interaction between two of the common genetic variants, VSX2 and PRPH2, affecting outer retinal thickness. We further explored the spatial genetic effect of our discovered genetic variants to capitalise on more of the available dimensionality of the phenotype data.
Methods
Ethics Statement
The UK Biobank study was conducted with the approval of the North-West Research Ethics Committee (ref 06/MRE08/65), in accordance with the principles of the Declaration of Helsinki, and all participants gave written informed consent. The research presented here has been conducted using the UK Biobank Resource under Application Number 2112.
UK Biobank Data
The UK Biobank is a cohort comprising ∼500,000 individuals recruited through the NHS registry at age 40-69 from across the UK. Individuals were were not selected on the basis of having disease, resulting in a broad cross-section of the UK population. For all participants, a number of baseline physical measurements were taken and an extensive questionnaire completed comprising information on lifestyle, demographic and socioeconomic factors. Additionally, blood and urine samples were collected and further tests including a heel-bone ultrasound, bio-impedance, hand-grip strength, spirometry, blood pressure and several cognitive tests were performed. Each individual was further genotyped and exome sequenced. The participants also agreed to ongoing linkage of their medical records [10].
Ophthalmic Measurements and Optical Coherence Tomography
A subset of 132,041 participants completed an eye examination. The measures collected included best corrected visual acuity using a logarithm of the minimum angle of resolution (logMAR) chart (Precision Vision, LaSalle, Illinois, USA). Visual acuity was measured with participants wearing their distance glasses at 4m, or at 1m if a participant was unable to read letters at 4m. Participants were asked to read from the top of the chart downwards, with the test terminated when two or more letters were read incorrectly. Further, refractive error was measured using a Tomey RC-5000 auto refraktometer (Tomey Corp., Nagoya, Japan) [11]. For each eye, up to 10 measurements were taken and the most reliable measure was automatically recorded. Intraocular pressure (IOP) was measured using an Reichert Ocular response Analyzer [12] from which corneal compensated IOP was calculated that accounts for rigidity of the cornea [13].
Optical Coherence Tomography Imaging Data
A further subset of 67,321 individuals underwent Spectral Domain Optical Coherence Tomography (SD-OCT) imaging. OCT imaging was conducted using the the Topcon 3D OCT1000 Mark II machine using 3-dimensional 6 × 6mm2 macular volume scan mode (512 A scans per B scan; 128 horizontal B scans). The imaging was completed on undilated eyes in a dark room following other eye measurements [13]. The right eye was consistently imaged first, and the majority of individuals had imaging repeated in their left eyes [12].
The images were processed using the industry standard Topcon advanced boundary software (TABS) [14]. This software performs automated retinal layer segmentation. It calculated thickness values for each of the retinal layers averaged across the macular field, as well as in each of the subfields of the early treatment of diabetic retinopathy study (ETDRS) grid [15]. The ETDRS is a nine-segment grid arranged as a bullseye and is commonly used during clinical assessment of the macula including the PRCs (fig. 1B).
Quality Control of Genotypic and Phenotypic Data
Participants who were recommended for exclusion from genetic studies by the UK Biobank were removed from the dataset. The most densely populated well-mixed population was selected for using principal component analysis. Individuals within a defined euclidean distance of the mean of the selected population, as identified by comparison to the HapMap Phase III study [16] in PC1-PC2 space were selected. Individuals were further excluded if they were related to third degree or more [17].
In addition to the genotypic quality control, rigorous quality control was applied to the phenotypic data. Exclusion/inclusion criteria were applied to the OCT images and the quantitative measures derived from them utilising methods previously implemented in Patel et al. (2016) [18]. In line with this method all participants with an OCT image quality score <45 were removed from the study. Further, individuals with values within the poorest 20% of the population in each of the OCT segmentation indicators were removed. These segmentation indicators include: Inner limiting membrane (ILM) indicator, a measure of the minimum localised edge strength around the ILM boundary across the entire OCT scan. ILM indicator is indicative of blinks, severe signal fading, and segmentation errors; Valid count, used to identify significant clipping in the z-axis of the OCT scan; Minimum motion correlation, maximum motion delta and maximum motion factor, all of which utilise the nerve fibre layer and total retinal thickness to calculate Pearson correlation and absolute differences between the thickness values from each set of consecutive B-scans. The lowest correlation and highest absolute difference in a scan define the resulting indicator values. These values identify blinks, eye motion artefacts and segmentation failures. It should be noted that the image quality score and segmentation indicators are often correlated with one another. Finally individuals with outlier values of refractive error were removed from the study. Outlier refractive error scores were defined as values lying outside one standard deviation of 1.5 times the inter-quartile range from the median. The final dataset included 31,135 individuals.
Genome-wide Association Studies
GWAS were implemented using an additive linear model in BGENIE [17]. The mean thickness of each component PRC layer across the ETDRS grid was used as the input phenotype (fig. 1C). Eye-specific covariates including image quality measures obtained from the OCT machine (listed above) and refractive error (calculated as spherical error + 0.5 × cylindrical error) were regressed out of thickness values for the left and right eye separately before a mean was made across the two eyes. Covariates including age, weight, sex, height, OCT machine ID and the first 20 genotype PCs were used in the model. A genome-wide significance threshold of P <5 × 10−8 was used. LD-score regression was implemented using LD SCore v1.0 [19].
Associated variant discovery set
The summary statistics from the GWAS of the ONL, IS and OS were meta-analysed using MTAG [20]. The resulting lowest p-value of the three meta-analysed phenotypes was considered the meta p-value. The meta-analysed PRC results were filtered using GCTA Conditional and Joint Analysis (COJO) [21] to select for independent SNPs that were more than 10Mb apart. Further filtering was applied and loci within 1.5Mb of one another were labelled as being in the same locus as indicated in the results tables where loci considered to be within the same locus are shaded the same colour, alternating grey and white (table S3). A plot of a magnified locus of interest was generated using LocusZoom [22].
Exome analysis
We performed gene burden testing using the UK Biobank (UKBB) 200,000 whole exome release to search for genes that showed an increased level of rare loss of function (LoF) or rare missense variants that were predicated as deleterious and probably damaging. Briefly, we ran the ensemble variant effect predictor (VEP) [23] on the multi sample VCF from the 200K UKBB release to allow us to classify variants and assign their consequence. We defined high confidence LoF variants using the VEP plugin LOFTEE and selected missense variants using SIFT “deleterious” and Polyphen “probably damaging” flags, all variants were filtered to have a minor allele frequency (MAF) of less than 1% prior to collapsing into gene based genotypes for burden testing. To define gene level genotypes we used the standard collapsing method where we considered all variants within the gene for each sample separately and defined zero, one and two (0,1,2) genotypes based on: zero (genotype 0) no rare variant of the given type present in the gene; one (genotype 1) more than zero rare variants of the given type present in the gene but all are in the heterozygous state; and two (genotype 2) if any of the rare variants of the given type where in the homozygous state. It is worth noting that by using this approach we do not account for compound heterozygous mutations or the overall variant load within individual genes but it does allow us to interrogate potential rare variant associations to the trait of interest at the gene level. To perform the association testing we treated the gene level genotypes for rare LoF and potentially damaging missense variants as standard genotypes within a linear model with covariates against the OCT derived PRC layers. The covariates included within the models and the overall sample selection was exactly the same as for the common SNP GWAS analysis. A significance threshold of P <5 × 10−5 was used.
Annotation of variants and geneset enrichment
Significant variants were manually annotated with associated gene and any prior ophthalmic or general phenotypic associations using ENSEMBL [24] and Open Targets Genetics platform [25]. Geneset and tissue enrichment analysis was conducted using DEPICT [26].
Comparison of Genetic Effect on Retinal Sub-fields
The foveal, intermediate and peripheral fields of the ETDRS grid were considered separately (fig. 2). A mean thickness of each of the PRC layers was calculated across each of these fields. A GWAS was then performed for each of these phenotypes. Eye specific covariates were regressed out of the individual eye measurements prior to a mean being made across left and right eyes, as described above. The same covariates were used in these GWAS as were in the simple phenotype GWAS. The GWAS were performed using BEGENIE [17].
The macula divided into three concentric fields for comparison utilising the ETDRS grid. The peripheral field is shaded white, the intermediate field is shaded pink and the foveal field is shaded grey.
To compare the genetic effect of each SNP on the different subfields, each retinal layer was considered at a time. For each PRC layer, comparisons were made between the genetic effect on: the foveal field compared to the intermediate field; the foveal field compared to the peripheral field; and the intermediate field compared to the peripheral field. For each comparison, SNPs were selected that were significantly associated with either of the two phenotypes being considered. A z-score was calculated for each of the SNPs comparing the genetic effect size for the two phenotypes. The z-score was defined as:
P-values were calculated for each of the z-scores and SNPs selected that were significant following bonferroni correction. The list of SNPs was clumped into loci based on linkage disequilibrium values. Each variant was annotated with all SNPs with r2>0.05 in a 1000kb window using PLINK and the locus with the lowest p-value in each window selected.
Visualisation of spatial genetic effect
A linear model was assessed of the effect of each of the loci identified as differentially affecting the concentric retinal fields on the thickness of each of the PRC layers in each of the nine segments in the ETDRS grid. The model included the same covariates as the GWAS of average retinal thickness and the effect size was plotted back into the ETDRS grid space.
Interaction models
Linear models were used to explore the presence of genetic interaction effecting the thickness of the three PRC layers. The same covariates were used as those used in the GWAS. Significance values were corrected for multiple testing using Benjamini and Hochberg correction and a false discovery threshold of < 0.15 was used.
Results
Initially we conducted three genome wide association studies (GWAS), one for the thickness of each component PRC layer - the ONL, IS and OS - averaged across both the ETDRS grid and across both left and right eyes (fig. 1C). This is analogous to many clinical uses of macular layer thicknesses. Quality control was applied to the OCT image data consistent with previous work [9, 18] and we used a well established method to select the European-associated PCA cluster using genetic data, which aims to minimise variation in non-genetic factors and genetic factors (See Methods). The resulting subset of samples that passed both image and genetic filters comprised of 31,135 individuals (table S1). The three GWAS (one for each component layer of the PRC), each identified a number of associated loci using the established genome-wide significance level of P <5 × 10−8 threshold. Following calculation of the linkage disequilibrium score for each of the three GWAS, there was minimal evidence for inflation due to residual population structure (fig. S1, table S2). SNP heritability of the thickness of each of the PRC layers was calculated: ONL h2=0.39, IS h2=0.18, and OS h2=0.21. Unsurprisingly we observed numerous loci that were significantly associated with the thickness of more than one of the PRC layers, and so we meta-analysed results using MTAG [20]. Following selection of independent loci using COJO [21], there were 111 loci significantly associated with the meta-analysed PRC layers’ thickness (Figure fig. 3, table S3 and table S4).
Manhattan plot of the photoreceptor cell layer thickness phenotype GWAS p-values. These result from meta-analysis using MTAG [20] of the summary statistics from GWAS of the thickness of the outer nuclear layer (ONL), inner segment (IS) and outer segment (OS). Variants are considered significantly associated if they reach genome wide significance (P <5 × 10−8). Loci of interest, selected through manual labelling by ophthalmic experts, are circled and labelled.
We immediately noticed that many of these loci were close to well known genes associated with rare ocular diseases (fig. 4 and table S5). 17 of the 111 loci overlap with a known rare retinal disease locus, which is highly unlikely to occur by chance (P <2.2 × 10−16 Hypergeometric distribution).
Some genes are also expressed significantly in other cell types not shown (e.g. RLBP1 is expressed also in Muller cells). Also, some genes (including those not shown) have a wider role in development, but expression is more restricted in adults (e.g. VSX2 is expressed in retinal progenitor cells in development, but mainly in bipolar cells in the mature retina).
We further leveraged the initial release of the UK Biobank exome dataset of 247,000 individuals, complementing the common variation GWAS with a rare variation gene burden test on the same phenotypes. We considered both strict loss of function mutations and overall missense mutations as two separate tests applied to each layer. Although there were no strong rare variant hits, setting a more liberal significance threshold (P <5 × 10−5) - consistent with other exome wide burden tests - we found 10 significant loci (table S6,fig. S2). Three of these loci overlap with rare eye disease, again unlikely to have been observed by chance (P=2.65 × 10−6 by hypergeometric distribution).
Overall between rare and common discovery we have found over 100 loci associated with PRC thickness in some manner, and these continuous quantitative phenotype associations have significant overlap to known rare ocular disease genes.
Selected Loci of Interest
Of these genetic variants, the locus with the smallest p-value and with a reasonably large effect size for a single locus was VSX2. VSX2 is a homeobox transcription factor that has been found expressed in the developing retina of humans, mice and zebrafish [27]. The gene is known to have an important role in ocular size, with pathogenic variants leading to microphthalmia [27]. This locus showed allelic heterogeneity in the GWAS results with the most significant genetic variant being rs1972565 (P = 5.23 × 10−191) whilst COJO analysis retained 4 other SNPs across the VSX2 locus. In addition, there are SNPs in nearby genes which could be long range regulatory SNPs for VSX2 rather than the closest genes we have assigned this to (fig. 5). Interestingly in previous work we identified genetic variants in neighbouring loci that had a statistically significant effect on the thickness of the inner retinal layers [9]. All this evidence points towards the VSX2 locus having multiple alleles which affect different aspects of eye development.
Association of the locus on chromosome 14 centred on the genetic variant rs1972565 (highlighted in purple). Other genetic variants found significantly associated with PRC thickness are labelled in grey.
Further there were several loci found associated with PRC layer thickness that are known to be involved in the retinoid cycle. The retinoid or visual cycle refers to the recycling of chromophore: photons of light isomerise 11-cis retinal to all-trans retinal in the photoreceptor outer segments; to be converted back to 11-cis retinal, a multi-stage process is entailed, which includes transport to the retinal pigment epithelium and then back to the photoreceptor [28]. RDH5, RLBP1 and RBP3, all found associated with PRC layer thickness, each have a known role in the retinoid cycle and are associated with rare inherited retinal diseases. RDH5 (rs3138142, P = 5.58 × 10 -60) encodes 11-cis retinol dehydrogenase, an enzyme responsible for conversion of 11-cis retinol to 11-cis retinal. The specific SNP in RDH5 (rs3138142, P = 5.58 × 10 -60) has previous associations with a number of ophthalmic phenotypes including the age one started wearing glasses, cataract and myopia. Pathogenic variants in this gene are known to cause fundus albipunctatus, a retinal dystrophy and form of night blindness that primarily affects the rods [29, 30]. RLBP1 (rs3825991, P = 2.70 × 10−31) encodes the cellular retinaldehyde binding protein, a transporter of 11-cis retinal [31, 32]. Mutations in this gene lead to a number of different retinal dystrophies [33]. RBP3 (rs111245635, P = 8.41 × 10−10) encodes the interphotoreceptor retinoid-binding protein which has a role in transporting retinoids between the RPE and the photoreceptor cells [34]. Additionally, rare variants within RPB3 cause recessive retinal dystrophy and high myopia [35]. Much of the retinoid cycle occurs in the supporting RPE layer, for example RDH5 is almost exclusively expressed in this layer [36], showing that one does not need gene expression in a cell type to impact its structure.
Two loci, SAG (rs7564805, P = 6.16 × 10−29) and GRK1 (rs9796234, P = 4.69 × 10−18), were found associated with PRC layer thickness and these genes have known involvement in Oguchi disease. Oguchi disease is a rare autosomal recessive disease. Features include the Mizuo-Nakamura phenomenon, a golden retinal sheen and night-blindness [37]. SAG encodes arrestin and GRK1 encodes rhodopsin kinase; both proteins are involved in shut-off of activated rhodopsin, enabling termination of the light response. When these proteins are dysfunctional, rhodopsin continues to activate the phototransduction cascade, shutting off the cyclic nucleotide-gated current, hyperpolarising the rods, and rendering them electrically unresponsive to further light photons. Thus patients primarily experience night blindness in this disorder. Studies have also found changes in the microstructure of the retinal layers in Oguchi disease patients [38, 39]. The presence of common variation in these genes suggests that there might be functional differences in rod physiological due to genetics present at the population level.
Further, several variants were found in the large heterogenous NPLOC4-TSPAN10 locus, and both rare and common variants were discovered in the OCA2 locus. Variants in these loci have previously been found associated with foveal structure [9] suggesting that overall retinal development has an impact on the PRC layer. Additionally the NPLOC4-TSPAN10 locus is associated with strabismus, the abnormal alignment of the eyes [40]. Strabismus is known to be caused by defects in the extraocular muscles, cranial nerves, refractive error and visual cortex processing. Our finding suggests there may be a role of retinal structure in strabismus which warrants further investigation utilising appropriate clinical datasets.
We also note that several loci significantly associated with the PRC layer thicknesses had prior associations to AMD. These included CFHR2 (rs410895, P=8.19 × 1022), HTRA1 (rs60401382, P=2.43 × 10−16) and RAX2 (rs76076446, P=1.07 × 10−13). AMD is a disease characterised by the atrophy of PRCs caused by the accumulation of lipid deposits in the RPE. Extensive recent work has previously been completed to study the relationship between PRC thickness and AMD [41, 42].
The discovery of loci associated with PRC thickness at genes encoding proteins involved in phototransduction, generating the photoreceptor electrical response to light, is of interest. These include GNAT2, SAG, GRK1 and RHO. This suggests there is a shared importance of these genes between phototransduction and structure. It may warrant further investigation into the relationship between the structure and function of these cells and these variants might provide interesting genetic instruments for the integration of the physiological process of light transduction with other features.
To further explore the biological pathways underlying PRC morphology, we used DEPICT [26] to perform gene set and tissue enrichment analysis. When applied to the loci associated with PRC thickness, numerous genesets reached statistical significance following false discovery rate (FDR) correction (table S7). These include ocular-associated genesets such as abnormal lens morphology, microphthalmia, disorganized retinal layers and abnormal ocular fundus morphology amongst others. Additionally, tissue enrichment analysis of the loci associated with PRC thickness resulted in significant associations with retina and eye tissue (table S8).
Additionally, there were a number of genes identified in the gene-burden test analysis with interesting prior associations. These include OCA2, also significantly associated in our GWAS analysis, which has a prior association to oculocutaneous albinism. Individuals with this condition are characterised by foveal hypoplasia, the absence of the formation of the foveal pit. Further, three genes, ABCA4, MYO7A and NR2E3, are shown to have an association to PRC thickness and have a prior association with retinal dystrophies. ABCA4 encodes ATP-binding cassette sub-family A member 4, a protein expressed almost exclusively in the PRCs. It is responsible for mediation of all-trans-retinal aldehyde across the PRC membrane, clearing it from the retina following phototransduction. ABCA4 is the gene most frequently implicated in monogenic retinal dystrophies [5, 43]. Biallelic pathogenic variants are the most common cause of of juvenile macular degeneration (Stargardt disease) and can also give rise to a cone-rod dystrophy. MYO7A encodes myosin VIIA, a protein found in inner ear and retina. Within the PRCs, this protein helps regulate transport of opsins, light-sensitive proteins required for vision. Mutations in MYO7A are known to cause Type 1 Usher Syndrome, a condition characterised by deafness and retinal degeneration [44]. NR2E3 encodes a nuclear receptor transcription factor that activates rod development and represses cone development [45]. Mutations in this gene are also associated with retinitis pigmentosa in addition to enhanced S-cone syndrome [46].
Fine grained topography of genetic variation
In addition to the clinically utilised measure of the mean PRC thickness taken across the ETDRS grid, TABS outputs the thickness of each of the PRC layers within each of the nine ETDRS subsections. To capitalise on this increased dimensionality of data, we explored how genetic variation affects PRC morphology in the different areas of the ETDRS segmentation grid. We chose to consider the macular field split into three concentric fields which align with the topography of the valley-like macula (fig. 2). A mean thickness across each of these concentric fields - the foveal, intermediate and peripheral fields - was taken for each of the PRC layers. These phenotypes were used as an input to GWAS and the SNPs which passed the significance threshold of 5 × 10−8 were taken forward for further analysis. Z-score analysis was applied to the resulting effect size to identify loci with significantly different effect on thickness of the different concentric fields of each PRC layer. Numerous (32) loci were found in each layer that affect the thickness of the concentric fields differentially (fig. 6, fig. S3).
The figure shows the effect size of the SNPs (shown by a blue-red gradient) with differential effect on the concentric fields of the retina, on the nine different segments of the ETDRS grid across the three PRC layers listed from top to bottom: ONL, IS and OS. Six loci are visualised: RREB1 (rs75757892), TMEM106B (rs13237518), RDH5 (rs3138142), MERTK (rs869016), CDHR1 (rs55798570), and RAX2 (rs76076446). Effect sizes are scaled for each genetic locus individually.
The locus showing the most significant difference in effect on the thickness of the ONL across the concentric fields was RREB1 (rs75757892) which had a differential effect on the foveal field (effect size=-1.12) compared to both the intermediate field (effect size=-0.55) and peripheral field (effect size=-0.20). RREB1 encodes a zinc finger transcription factor, the Ras-responsive element binding protein 1 [47]. Genetic variants at RREB1 have prior associations to various ocular phenotypes including the age one starts wearing glasses and spherical power.
The locus showing the most significantly different genetic effect on the IS thickness across the macula was TMEM106B (rs13237518) which had a significant difference between the foveal field (effect size=-0.23), and the peripheral field (effect size=-0.02). TMEM106B is known to be involved in frontotemporal lobar degeneration, a neurodegenerative disease [48]. There is growing interest in the use of OCT derived retinal thickness measures as biomarkers for neurodegenerative diseases [49, 50] making this finding of potential interest.
The locus with the most significantly different genetic effect on the thickness of the OS was RDH5 also found in the initial meta-GWAS. RDH5 (rs3138142) shows a significantly different effect on the thickness in the peripheral field (effect size=0.81) compared to both the foveal field (effect size=1.37) and the intermediate field (effect size=1.28). As described above RDH5 encodes 11-cis retinol dehydrogenase, an enzyme that is part of the retinoid cycle. It has been previously associated with the retinal dystrophy, fundus albipunctatus.
Notably several further loci were discovered to differentially affect the concentric fields of the retina that had prior associations with retinal dystrophies. These include MERTK (rs869016) which shows a significant difference of effect across the concentric fields of the OS. It has a significantly larger effect size on the fovea (effect size=-0.31) as compared to the peripheral field (effect size=0.02). MERTK encodes MER proto-oncogene, tyrosine kinase, a transmembrane protein localised to the RPE. MERTK has a known role in the phagocytosis of the OS by the RPE. Mutations in MERTK have been associated with retinitis pigmentosa and early onset retinal dystrophies [51].
Additionally CDHR1 (rs55798570) shows a significantly different effect size on the thickness of the IS in the foveal field (effect size=0.13) compared to both the intermediate field (effect size=0.002) and the peripheral field (effect size=-0.06). CDHR1 encodes cadherin-related family member 1, a cell adhesion protein that localises to the junction of the IS and OS. Genetic variants in CDHR1 have been associated with a number of inherited retinal dystrophies including retinitis pigmentosa and cone-rod dystrophy [52].
Further RAX2 (rs76076446), shows a different effect on the thickness of the IS at the fovea (effect size=0.31) compared to the peripheral field (effect size=-0.09). RAX2 encodes retina and anterior neural fold homeobox 2, which has a known role in eye development. Pathological variants in this locus are known to cause cone and rod dysfunction and retinitis pigmentosa, a rod-cone disease [53].
Each of these retinal dystrophies have distinct and sometimes opposing spatio-temporal patterns in which they affect the rods and cones. Some dystrophies, including retinitis pigmentosa, affect the rods first, with the cones affected later. Patients have problems with night vision, then peripheral vision, and finally central vision. Fundus albipunctatus also affects rod function, whilst cones are usually spared. Others such as cone-rod dystrophy affect the cones first and the rods afterwards. Cones are found at the highest density at the fovea, whilst rod density peaks more peripherally. There appears to be a parallel between the spatial genetic effects of these loci on thickness of the PRC and the organisation of PRCs affected by the diseases associated with these loci.
Genetic Interactions
As there are both clear common biological pathways associated with multiple discovered variants, and overlaps with known rare disease loci we wanted to explore whether genetic interactions were affecting the retinal phenotype. Genetic interactions of common variation are challenging to discover and characterise due to the large space of possible interactions. To restrict our space we considered interactions in the following categories: (1) Interactions of the developmentally important VSX2 locus internally and with all other common genetic variants; (2) Interactions between the Retinoid cycle loci; Interactions of the Oguchi syndrome loci; (4) Interaction of all loci with genetically defined sex; and (5) Interactions between rare burden levels in genes and nearby common SNPs.
Overall we do not find widespread interaction terms that would pass multiple testing criteria. There was no statistically significant interaction term between the genetic variants with prior association to Oguchi disease, or between those with a prior association with the retinoid cycle. Additionally, there was no evidence of an interaction between any of the 111 identified loci and genetically defined sex. There was no evidence for interactions between rare burden levels in genes and nearby common SNPs; however this analysis was limited by the sample size of individuals with heterozygous alternative allelic states. Expansion of this analysis in a cohort containing a higher proportion of individuals with rare disease phenotypes would have more power to identify interactions.
However, following Benjamini and Hochberg correction for multiple testing, and using a False Discovery threshold level of < 0.15, four variants show an interaction with the lead VSX2 SNP, rs1972565 in relation to their effect on ONL thickness. Two of the four variants (rs112145470 and rs28488340) are within the same locus as VSX2, consistent with the already noted allelic hetreogenity at this locus. The other two loci are DYNLRB2 (rs7206532) and PRPH2 (rs375435). DYNLRB2 is expressed in retinal Muller cells and astrocytes, and encodes dyenin light chain roadblock type 2. Mutations in this gene are associated with deafness and hepatocellular cancer. PRPH2 encodes peripherin 2 which is localised to the photoreceptors. Mutations in this gene are known to result in forms of retinal dystrophy and macular degeneration [54]. Upon stratification of the population by their genotype at VSX2 and PRPH2, a change in ONL thickness can be observed, specifically between PRPH2 homozygous reference and homozygous alternative genotype in those with the VSX2 homozygous reference genotype (fig. 7,fig. S4). In individuals with PRPH2 homozygous reference, a thinner ONL is found in the superior segments of the macula, and thicker ONL in the inferior segments. In individuals with PRPH2 homozygous alternative, the inverse pattern is seen. This is of particular interest as there is often a difference in PRC degeneration across the macula eqautor in retinal dystrophies. PRPH2 has a doccumented specturm of penetrance and associated disease phenotypes [55, 56]. This finding suggests one of the possible mechanisms by which the effects of variants at the PRPH2 locus are modified is via the genotype at the VSX2 locus.
The mean thickness of the ONL at each of the nine segments of the ETDRS grid across individuals stratified by their genotype at VSX2 and PRPH2. The thickness is mean normalised per VSX2 genotype.
Discussion
Here we have performed the largest and most detailed GWAS and rare variation association of PRC layer thickness to date, leveraging the high-dimensional spatial resolution available within retinal OCT. We initially explored a simple phenotype of the mean thickness of each of the PRC layers - the ONL, IS, OS - across the entire macula field as defined by the ETDRS grid. These phenotypes are advantageous in that they are aligned with those used in clinical practice. Using these phenotypes, we identified 111 independent loci significantly associated with one or more of the component PRC layers following GWAS meta analysis and 10 loci associated with rare variation following gene burden testing. Many of the loci identified here had prior associations to ocular phenotypes and the loci collectively were enriched for known rare ocular diseases.
The collection of rare disease loci discovered in this study of common, physiologically healthy individuals without obvious eye diseases shows the continuum between common and rare variation, and rare variation which leads to disease. In this study we explicitly looked for interactions between common variation at different loci and common and rare variation at the same loci. We found several relatively weak interactions between a SNP near the well known eye development gene, VSX2, and two SNPs within this locus as well as two SNPs near PRPH2 and DYNLRB2. No evidence was found for an interaction between common and rare genetic variants but this was partially due to limitations of the dataset and its absence of adequate numbers of individuals with rare variation and subsequent lack of statistical power.
Genetic discovery also found sets of loci associated with PRC thickness with prior associations to different ocular phenotypes. Four loci (GNAT2, SAG, GRK1 and RHO) have a prior association to phototransduction, the generation of an electrical response to light in the PRC. Further three variants were found in or near genes (RDH5, RLBP1 and RBP3) with prior associations to the retinoid cycle, the process by which the chromophore is regenerated following photoisomerisation by light. In both cases, this finding suggests a relationship between the form and function of these cells which may warrant further investigation. Whilst it might be intuitive that genes encoding structural and matrix proteins are associated with photoreceptor layer thicknesses, it is interesting that genes encoding proteins involved in bringing about or shaping the photoreceptor electrical light response to light are also implicated, potentially suggesting a role for visual signalling in development or maintenance of outer retinal structure.
There were several loci associated with specific disease. For example both SAG and GRK1 in addition to an association with PRC thickness have a known involvement in Oguchi disease, a rare autosomal eye condition which causes visual impairment. As well as the overlap between rare eye genetic disease and these common variants affecting retinal structure, there are some interesting, less expected findings. For example, three loci - MYEOV (rs10737153), TNS1 (rs201030469) and FGFR2 (rs17102399) – have impacts on both PRC layer thickness and risk of diagnosis for breast cancer. Additionally TMEM106B, which had a particularly notable spatial genetic effect, has a known association to a form of dementia. In both cases, the link between morphology and disease suggests the possibility of using OCT derived PRC morphology measures as biomarkers for complex diseases in the future. Additional datasets enriched for individuals with these diseases with OCT data will be required to further explore these concepts.
This study has offered the ability to study the genetics underlying complex retinal structures. Whilst OCT does not directly offer single-cell resolution, the quantitative values extracted from the images in the photoreceptor layer describe a single cell layer and capture variation in sub-cellular structures. It is interesting to note that some of the strong hits have gene expression in a different cell type. For example RDH5 which is strongly associated with thickness of all three PRC layers is expressed in the RPE [57]. Therefore we are identifying genetic variation impacting sub-cellular morphology in a between cell type manner, highlighting the complex route between genetic variation, gene function and ultimate physiological differences. Imaging offers a methodology for interrogating such processes at a high-dimensional scale in a non-invasive manner and the protocol used here could be applied to other organs of the human body.
Data Availability
The genetic and phenotypic UK Biobank data are available upon application to the UK Biobank.
Data Availability
The genetic and phenotypic UK Biobank data are available upon application to the UK Biobank.
Author Contribution
HC conceived the analysis, performed the analysis and wrote the paper. TWF performed analysis and wrote the paper. PJP, APK and ARW provided advice on analysis. OAM performed analysis and oversaw the work. EB conceived the analysis, wrote the manuscript and oversaw the work. All authors reviewed the final manuscript.
Competing Interests
APK has acted as a consultant to Abbvie, Aerie, Google Health, Novartis, Reichert, Santen and Thea.
Supporting information
Results are presented as mean ± standard deviation.
Linkage score disequilibrium metrics for the three GWAS of the photoreceptor layers following analysis using LDSCore [19]. Metrics include lambda genomic control (GC), linkage disequilibrium score (LDSC) intercept and LDSC ratio.
Each SNP is annotated with the associated gene and any prior associations to ocular or non-ocular phenotypes. Variants considered to be representative of a single locus, examples of allelic heterogeneity, are highlighted in the same colour alternating white and grey. Full results including effect size, effect allele specification and standard error are available in Supplementary Table S4.
Extended results of GWAS of the three photoreceptor cell layers, the outer nuclear layer (ONL), inner segment (IS) and outer segment (OS) following meta-analysis as well as the results selected by meta-analysis using MTAG, labeled MTAG. Each locus found significantly associated following meta-analysis is listed alongside the chromosome (Chr) and location (BP). The reference and alternative alleles are listed alongside the allele frequency (AF) of the alternative allele (A1). The effect size refers to the effect of having an additional copy of the A1 allele and the standard error (SE) is further listed.
(A) The quantile-quantile plot (qq-plot) for the GWAS of ONL thickness prior to meta-analysis (Lambda GC = 1.15, Intercept = 1.04, Ratio = 0.13). (B) The qq-plot for the GWAS of IS thickness prior to meta-analysis (Lambda GC = 1.07, Intercept = 1, Ratio <0). (C) The qq-plot for the GWAS of OS thickness prior to meta-analysis (Lambda GC = 1.09, Intercept = 1.02, Ratio = 0.15).
Several genes are involved in bringing about or shaping the photoreceptors’ electrical response to light (RHO, GNAT2, SAG, GRK1) or in recycling of the light-sensitive chromophore (ABCA4, RLBP1, RBP3, RDH5). Others are involved in maintaining outer retinal structure or in retinal and eye development. Superscript letters denote the relevant analyses in which the association was found: aGWAS; bdifferential association with ONL across concentric fields; cdifferential association with IS across concentric fields; ddifferential association with OS across concentric fields; eexome analysis.
Manhattan plot of p-values resulting from exome wide loss of function burden testing analysis for the thickness of each of the retinal layers (A) ONL, missense model, (B) IS, missense model, (C) OS, missense model, (D) ONL, loss of function model, (E) IS, loss of function model, (F) IS, loss of function model. Variants are considered significantly associated if they reach a a p-value threshold of P <5 × 10−5.
List of genes associated with one of the three photoreceptor layers, the outer nulcear layer (ONL), inner segment (IS) and outer segment (OS) following gene burden testing. A significance threshold of P <5 × 10−5 is used. For each gene the gene name and chromosome is listed. Each gene is also annotated with any prior associations to ocular and non-ocular phenotypes. The number of individuals with the homozygous loss of function allele is listed as N. Two modes were tested, a missense (MS) model and a loss of function (LoF) model.
The results of geneset enrichment analysis using DEPICT applied to the GWAS results of the meta analysed PRC layers. All genesets which were significantly enriched following correction for false discovery rate are listed below.
The results of tissue enrichment analysis using DEPICT applied to the GWAS results of the meta analysed PRC layers. All Tissues which were nominally significantly enriched are listed below.
Plots depicting the residuals from models comparing the effect size of SNPs on the thickness of each retinal varinats in different concentric retinal fields. Variants with significantly different effect sizes on the two different areas tested, as determined by a z-score, are highlighted in pink. (A) ONL, fovea compared to intermediate; (B) ONL, fovea compared to peripheral; (C) ONL, intermediate compared to peripheral; (D) IS, fovea compared to intermediate; (E) IS, fovea compared to peripheral; (F) IS, intermediate compared to peripheral; (G) OS, fovea compared to intermediate; (H) OS, fovea compared to peripheral; (I) OS, intermediate compared to peripheral.
List of SNPs with a significant z-score describing the differential effect on the ONL thickness at the foveal (F), intermediate (I) and peripheral (P) fields. The field (F1 or F2) and corresponding effect size from GWAS of thickness in each field are listed alongside the p-value of the comparative z-score. Each genetic variant is also annotated with associated gene and any ocular and non-ocular phenotypes previously associated with it. The different concentric comparisons are separated by bold horizontal lines.
List of SNPs with a significant z-score describing the differential effect on the IS thickness at the foveal (F), intermediate (I) and peripheral (P) fields. The field (F1 or F2) and corresponding effect size from GWAS of thickness in each field are listed alongside the p-value of the comparative z-score. Each genetic variant is also annotated with associated gene and any ocular and non-ocular phenotypes previously associated with it. The different concentric comparisons are separated by bold horizontal lines.
List of SNPs with a significant z-score describing the differential effect on the OS thickness at the foveal (F), intermediate (I) and peripheral (P) fields. The field (F1 or F2) and corresponding effect size from GWAS of thickness in each field are listed alongside the Bonferroni adjusted p-value of the comparative z-score. Each genetic variant is also annotated with associated gene and any ocular and non-ocular phenotypes previously associated it. The different concentric comparisons are separated by bold horizontal lines.
Comparison of normalised outer segment (OS) thickness and outer nuclear layer thickness (ONL) in population subsetted by their genotype at VSX2 (rs1972565) and PRPH2 (rs375435).
Acknowledgments
The authors are extremely grateful for the selfless participation of individuals in the UK Biobank used in this study and the staff managing these cohorts. HC has received support from EMBL and is supported by Novo Nordisk Foundation (grants NNF17OC0027594 and NNF14CC0001). APK is supported by a UKRI Future Leaders Fellowship and an Alcon Young Investigator Award. TWF and EB are supported by EMBL. OAM is supported by the Wellcome Trust (206619/Z/17/Z). The UK Biobank Eye and Vision Consortium members are: Naomi Allen, Tariq Aslam, Denize Atan, Sarah Barman, Jenny Barrett, Paul Bishop, Graeme Black, Tasanee Braithwaite, Roxana Carare, Usha Chakravarthy, Michelle Chan, Sharon Chua, Alexander Day, Parul Desai, Bal Dhillon, Andrew Dick, Alexander Doney, Cathy Egan, Sarah Ennis, Paul Foster, Marcus Fruttiger, John Gallacher, David Garway-Heath, Jane Gibson, Jeremy Guggenheim, Chris Hammond, Alison Hardcastle, Simon Harding, Ruth Hogg, Pirro Hysi, Pearse Keane, Peng Tee Khaw, Anthony Khawaja, Gerassimos Lascaratos, Thomas Littlejohns, Andrew Lotery, Phil Luthert, Tom MacGillivray, Sarah Mackie, Bernadette McGuinness, Gareth McKay, Martin McKibbin, Tony Moore, James Morgan, Richard Oram, Eoin O’Sullivan, Chris Owen, Praveen Patel, Euan Paterson, Tunde Peto, Axel Petzold, Nikolas Pontikos, Jugnoo Rahi, Alicja Rudnicka, Naveed Sattar, Jay Self, Panagiotis Sergouniotis, Sobha Sivaprasad, David Steel, Irene Stratton, Nicholas Strouthidis, Cathie Sudlow, Zihan Sun, Robyn Tapp, Dhanes Thomas, Emanuele Trucco, Adnan Tufail, Ananth Viswanathan, Veronique Vitart, Mike Weedon, Katie Williams, Cathy Williams, Jayne Woodside, Max Yates, Jennifer Yip, Yalin Zheng.
References
