Genetic variants and haplotypes in fibulin-5 (FBLN5) are associated with pseudoexfoliation glaucoma but not with pseudoexfoliation syndrome

Abstract Pseudoexfoliation (PEX) is a multifactorial age-related disease involving deposition of extracellular proteinaceous aggregates on anterior ocular tissues. The present study aims to identify functional variants in fibulin-5 (FBLN5) as risk factors for the development of PEX. Thirteen tag single-nucleotide polymorphisms (SNPs) in FBLN5 were genotyped using TaqMan SNP genotyping technology to identify association between SNPs of FBLN5 and PEX in an Indian cohort comprising 200 control and 273 PEX patients (169 PEXS and 104 PEXG). Functional analysis of risk variants was done through luciferase reporter assays and electrophoretic mobility shift assay (EMSA) using human lens epithelial cells. Genetic association and risk haplotype analysis showed a significant association of rs17732466:G>A (NC_000014.9:g.91913280G>A) and rs72705342:C>T (NC_000014.9:g.91890855C>T) within FBLN5 as risk factors with the advanced severe stage of the disease, pseudoexfoliation glaucoma (PEXG). Reporter assays showed allele-specific regulatory effect of rs72705342:C>T on gene expression, wherein, construct containing the risk allele showed a significant decrease in the reporter activity compared with the one with protective allele. EMSA further validated higher binding affinity of the risk variant to nuclear protein. In silico analysis predicted binding sites for two transcription factors, GR-α and TFII-I with risk allele at rs72705342:C>T, which were lost in the presence of protective allele. The EMSA showed probable binding of both these proteins to rs72705342. In conclusion, the present study identified the novel association of two genetic variants in FBLN5 with PEXG but not with PEXS, distinguishing between the early and the later forms of PEX. Further, rs72705342:C>T was found to be a functional variant.


Introduction
Pseudoexfoliation (PEX) is a progressive systemic disease of protein aggregation that involves the deposition of extracellular fibrillar material called PEX fibrils in the tissues of various organs, which manifests prominently in the eye. The initial stage of deposition termed as pseudoexfoliation syndrome (PEXS) progresses to the severe form, pseudoexfoliation glaucoma (PEXG) in almost half of PEXS-affected individuals [1]. The deposits impede the aqueous humour outflow pathways building up the intraocular pressure (IOP), which damages the retinal ganglion cell axons in the optic nerve leading to gradual blindness that are the hallmark features of PEXG [2,3]. PEXS accounts for 25-70% of all open-angle glaucoma cases depending on the country [4]. The prevalence of PEX is highly variable across the globe with 0.0% cases in Eskimos, 1.8% in the United States, 3.8-6.0% in India, 11.0-35.0% in East Africa, 5.0-25.0% in the Scandinavians and 38.0% in Navajo Indians [5][6][7][8].
The PEX fibrils or deposits comprise various extracellular proteins, membrane proteins and molecular chaperones, such as elastin, tropoelastin, fibrillin-1, amyloid P, lysyl oxidase like-1, vitronectin, fibulins, clusterin, latent transforming growth factor binding proteins and other components of elastic fibre system, and deregulation of these proteins has been observed at molecular levels in PEX individuals [9][10][11]. Previously, we reported the down-regulation of the extracellular matrix (ECM) protein, fibulin-5 (FBLN5) in the lens capsule of PEXS individuals [12]. FBLN5 encodes a secreted extracellular calcium-binding protein highly crucial for assembly of elastic fibres. FBLN5 deposits the cross-linking enzyme, lysyl oxidase-like 1 (LOXL1) in the ECM where it cross-links tropoelastin monomers to elastic fibres and is also involved in other essential cellular functions, such as cell matrix adhesion, integrin-dependent regulation of reactive oxygen species (ROS) in the ECM, regulation of cell receptor signalling, endothelial-to-mesenchymal transition and modulation of matrix proteases [13].
We have previously reported a novel genetic association of two variants, rs7149187 and rs929608 residing in the 5 -UTR and 10th intron, respectively, in FBLN5 as risk for PEX. However, these variants did not display any functional regulatory role in reporter assays. The main purpose of the present study was to identify putative regulatory variants in FBLN5 which could be explored further to find out whether they are functional variants playing pivotal role in PEX pathogenesis. In the present study, using a Tag SNP genotyping approach, we identified novel genetic association of two intronic variants in FBLN5 as risk factors for PEXG, one of which was found to be a functional variant.

Study subjects' selection and recruitment
The following study was approved by the Institutional Biosafety and Human Ethics committee of National Institute of Science Education and Research and adhered to the tenets of Declaration of Helsinki. Participants were recruited at Sri Sri Borda Hospital, Bhubaneswar, India by our clinical collaborator Dr Pranjya Paramita Mohanty. All the study participants underwent a detailed ocular examination including slit lamp microscopy, ocular biometry, Goldman applanation tonometry, +90D biomicroscopic fundus evaluation and four-mirror gonioscopy. Cataract patients aged above 40 years with clinically evident PEX like material over lens capsule (LC) and pupillary ruff having untreated IOP < 21 mmHg without any visual field defects were included under PEXS group and those with untreated IOP > 21 mmHg, glaucomatous nerve head damage with repeatable field defects corresponding to disc damage were included under PEXG group. Patients with corneal or retinal pathology precluding reliable visual field were excluded from PEXG group. Patients with systemic diseases, such as diabetes, were excluded from the study. Cataract patients aged above 40 years without PEXS or PEXG, with untreated IOP < 21 mmHg having normal discs and visual field were included as controls. Written informed consent was obtained from all participants. Age-sex-matched subjects were chosen for the study.

SNP selection and TaqMan genotyping
For higher genetic coverage of the FBLN5 gene, tag SNP genotyping approach was chosen. Thirteen tag SNPs within FBLN5 were chosen based on 1000 genomes HapMap South Asian GIH dataset with a pair-wise tagging of r 2 > 0.9 and minor allele frequency (MAF) > 0.1. The selected 13 single-nucleotide polymorphisms (SNPs) spanned the FBLN5 gene and were genotyped in 169 PEXS, 104 PEXG and 200 age-sex matched control subjects (Supplementary Table S1). Peripheral blood was collected from the study subjects and genomic DNA was extracted through phenol-chloroform method. The SNPs were genotyped using TaqMan SNP genotyping assays (Applied Biosystems, Carlsbad, CA, U.S.A.) on Quantstudio 7 (Thermo Fisher Scientific, U.S.A.). Each 5 μl PCR reaction mix consisted of 20 ng DNA sample, 2.5 μl of TaqMan master mix 2X, 1.25 μl of TaqMan SNP assay, and the volume was made up with nuclease-free water. The PCR conditions were initial denaturation at 95 • C followed by 40 cycles of denaturation at 95 • C for 15 s and annealing and extension at 60 • C for 1 min. The data were analysed using the instrument software and TaqMan Genotyper software.

Plasmid construction and Luciferase reporter assays
To test the functional effect of the SNPs, two reporter vectors, pGL4.23 with minimal promoter of the reporter and pGL3 containing the FBLN5 core promoter region were used. The region from −675 bp to ATG start codon of FBLN5 gene was amplified using Phusion High Fidelity DNA polymerase (NEB) and cloned into pGL3 basic luciferase reporter vector using BglII and NcoI (NEB). DNA fragments (29 bp long) surrounding the SNPs (Supplementary Table  S2) with either of the alleles at the centre were cloned into the reporter vectors using KpnI HF and XhoI (NEB). HLE B-3 cells were seeded in a 24-well plate and at 80% confluency, the cells were transiently transfected with 500 ng of the constructs along with 5 ng of Renilla vector (pGL4.74) using lipofectamine (Thermo Fisher Scientific, U.S.A.). After 24-h post transfection, cells were harvested and luciferase activity assayed using Dual-Luciferase Reporter assay system (Promega, U.S.A.). The Firefly luciferase activity from each construct was normalized to Renilla luciferase activity, and the ratio has been plotted as per cent luciferase activity relative to that of empty vector (taken as 100%). Data represent at least three independent experiments.

Electrophoretic mobility shift assays
A total of 29 base pairs sense (S) and antisense (A) oligonucleotides encompassing the rs72705342:C>T (NC 000014.9:g.91890855C>T) were synthesized for performing the EMSAs (Supplementary Table S2). The oligos were synthesized with their 5 -end labelled with biotin and unlabelled oligos were procured as well. The oligos were annealed by incubating the mix of complementary strands at 95 • C for 5 min followed by gradually cooling down the mix to room temperature. Nuclear extract from HLE B-3 cells was prepared using the NE-PER Nuclear and Cytoplasmic Extraction Reagents kit (Thermo Fisher Scientific, U.S.A.). The binding reaction included poly (dI.dC) as non-specific competitor DNA. For competition experiments, a 200-fold excess of unlabelled oligonucleotides was included in the pre-incubation mixture. For supershift assays, EMSA-specific antibodies for TFII I (sc-46670X, Santa Cruz, U.S.A.) and GR-α (PA1516, Invitrogen) were pre-incubated with the nuclear extract for 1 h before adding the final reaction mixture. The complexes after incubation were resolved on 6% native polyacrylamide gels and transferred to nylon membranes and developed. The EMSA was performed with the Lightshift Chemiluminescent EMSA kit (Thermo Fisher Scientific, U.S.A.). Detection was done using Fusion Solo S Chemi-Doc (Vilber Lourmat) and gel shifts were quantified with the Evolution Capt software (Vilber Lourmat Fusion Solo S).

In silico analysis
The PROMO software (http://alggen.lsi.upc.es/cgi-bin/promo v3/promo/promoinit.cgi?dirDB=TF 8.3 was used to identify transcription factor binding sites to the region surrounding rs72705342C>T. The sequences of genomic region flanking the 'T' allele or 'C' allele of the SNP rs72705342 (+ − 15 bp) were used as input.

Genetic and statistical analysis
Age-sex matched samples were taken for the experiments. The matching was done by performing the Student's t-test between the groups. No data were missing for the participants. The allelic association tests, Hardy-Weinberg equilibrium (HWE), and logistic regression analysis for covariates were done using PLINK. Haplotype analysis and linkage disequilibrium (LD) analysis were done using Haploview V4.2. Statistical significance of group-wise results was analysed using Student's t-test and P<0.05 was considered as statistically significant. The Bonferroni and Holm correction was applied for multiple pair-wise comparisons. All experiments were done at least three times independently. Data are presented as mean + − SEM.

Demographics of the study subjects
A total of 273 PEX (169 PEXS and 104 PEXG) and 200 age-sex matched control subjects participated in the present study. The demographics of the study subjects are shown in Table 1. The mean age in years + − SD of controls, PEXS, and PEXG were 70.17 + − 7.17, 71.09 + − 7.29 and 70.21 + − 7.39, respectively. The age-range of controls, PEXS, and PEXG was 60-90 years, 50-90 years and 60-92 years, respectively. Of the study participants, 36.1% were females. About 171 females (80 control, 64 PEXS and 27 PEXG) and 302 males (120 control, 105 PEXS and 77 PEXG) participated in the study.

Intronic variants, rs72705342 and rs17732466, within FBLN5 are genetically associated with PEXG
Thirteen tag SNPs ( Figure 1) within FBLN5 were genotyped in 169 PEXS, 104 PEXG and 200 age-and sex-matched   control subjects. All the studied SNPs passed the HWE test set at a default significance threshold of P≤0.001. Allele frequencies, odds ratio (OR) and statistical significance of the genotyped FBLN5 variants are presented in Table 2. Two variants, NC 000014.9:g.91913280G>A (rs17732466:G>A) and NC 000014.9:g.91890855C>T (rs72705342:C>T) located in the 4th and the 6th introns of FBLN5, respectively, were found to be significantly associated with PEXG with the risk alleles being 'G' (P=0.04) and 'C' (P=0.02), respectively. Risk analysis showed that the minor alleles ' A' at rs17732466 and 'T' at rs72705342 confer a protective effect with an OR of 0.66 (95% CI: 0.43-0.99) and 0.60 (95% CI: 0.39-0.93), respectively. However, none of the studied variants showed a significant association with PEXS. Genotypic distribution of the variants in controls, PEXS, and PEXG is presented in Table 3. None of the SNPs showed   any genotypic association with PEXS. However, individuals with the risk genotype 'CC' at rs72705342 showed higher susceptibility to having PEXG compared with individuals carrying the 'TT' genotype (P=0.04).

Genetic association of haplotypes was observed with PEXG
The LD pattern across the 13 studied FBLN5 SNPs is shown in Figure 2A. The confidence interval algorithm (Gabriel et al.) defined two LD blocks. The frequency of haplotype 'T-G' in Block 1 (rs72705342-rs2267995) was found to be significantly lower in PEXG (0.17, P=0.03) but not in PEXS (0.23, P=0.70) compared with control (0.24) ( Table 4). None of the haplotypes in Block 2 (rs917908-rs2267997-rs2498835) showed association with either PEXS or PEXG (Supplementary Table S3). Also, the haplotype analysis was done for the SNPs that were significantly associated with   Figure 2B). The frequency of haplotype 'C-G' at 'rs72705342-rs17732466' was significantly higher in PEXG (0.76, P=0.03) compared with controls (0.67) but was not associated with PEXS (0.69, P=0.6, Table 4).

rs72705342 shows allele-specific regulatory effect
To evaluate the putative regulatory effect of the regions containing these SNPs, luciferase reporter assays were performed. The 29 bp regions flanking the SNPs were cloned upstream of the minimal promoter in pGL4.23 and transiently transfected into HLE B-3 cells. The cells containing constructs with either 'G' (P=0.1) or ' A' (P=0.5) allele at rs17732466 did not show any differential luciferase activity compared with empty vector. Also, no significant changes in luciferase activity were observed between the alleles (P=0.4) ( Figure 3A). On the other hand, the alleles at rs72705342 element showed significant allele-specific changes in the luciferase activity. The presence of the protective allele 'T' at rs72705342 significantly increased the expression of luciferase compared with the construct with the risk allele 'C' (P=0.03) or the empty vector (P=0.001), the p-values remained significant at P=0.04 and P=0.007, respectively, after the correction for multiple pair-wise comparisons. No significant difference was observed between empty vector and construct with 'C' allele ( Figure 3B). Further, to assess the direct effect of rs72705342 element on FBLN5 promoter, the core promoter of FBLN5 [24] was cloned into pGL3 basic vector and the 29 bp rs72705342 loci with either allele 'T' or allele 'C' was cloned upstream of the FBLN5 promoter. These constructs were transiently transfected into HLE B-3 cells. The reporter activity showed that the alleles at rs72705342 showed an allele-specific effect on the FBLN5 promoter ( Figure 3C). Change in allele 'C' to 'T' at rs72705342 showed an increased luciferase activity (P=0.02, corrected P=0.03) implying an allele-specific regulatory effect of the rs72705342 element on FBLN5 promoter.

Risk allele 'C' at rs72705342 showed greater protein binding affinity compared with 'T' allele
The electrophoretic mobility shift assays (EMSA) were performed to study specific DNA-protein interactions at rs72705342. The EMSA yielded specific shifted bands, and shifts could be competitively inhibited by excess of unlabelled oligonucleotides ( Figure 4A). Quantitative analysis of the shifted bands showed greater protein binding to the sequence containing the risk allele 'C' compared with that containing allele 'T' (P=0.04) implying a differential transcription factor binding at rs72705342 ( Figure 4B and Supplementary Figure S1). Competitive EMSA on the labelled rs72705342 'C' probe with increasing concentration of unlabelled rs72705342 'C' oligo (50-, 100-, 200and 400-fold excess) showed progressive reduction of the shifted band ( Figure 4C). In silico analysis using PROMO software predicted binding of ten transcription factors to the region flanking rs72705342 (+ − 14 bp), i.e., glucocorticoid receptor α (GR-α), activating enhancer binding protein 2α (AP-2αA), CCAAT/enhancer-binding protein β (C/EBPβ), c-Jun, nuclear factor-1 (NF-1), estrogen receptor α (ER-α), forkhead box P3 (FOXP3), retinoid X receptor α (RXR-α), retinoic acid receptor α (RAR-α), CAAT box transcription factor (CTF) and c-ETS-2. Binding of only one transcription factor, TFII I was predicted to be affected by variation at rs72705342. Change in allele from 'C' to 'T' at rs72705342 predicted a loss of binding site for TFII I transcription factor. Also, although GR-α had binding sites in the 29 bp sequence, an additional binding site for GR-α was created in the presence of 'C' allele ( Figure 4D). To check the binding of TFII I and GR-α to rs72705342 'C' , the EMSA was performed using antibodies specific to either of the transcription factors. On pre-incubating the nuclear extract with TFII I antibody, a decreased intensity of the shift was noted which did not happen in the presence of non-specific HSF1 antibody ( Figure 4E). Further, as shown in Figure 4F, a reduced intensity in the shift was observed in the presence of GR-α antibody as well but not in the presence of the non-specific HSF1 antibody. However, contradictory to the in silico prediction, even with rs720705342 'T' probe, a reduction in the intensity of shift was observed in the presence of TFII I. A decreased intensity in shift was observed with rs72705342 'T' probe in the presence of GR-α antibody as well (Supplementary Figure S2A and 2B). These findings suggest that the protein-DNA complexes binding to rs72705342 element might comprise the TFII I and GR-α proteins.

Discussion
In the present study, we investigated the association of FBLN5 variants with PEX by Tag-SNP genotyping approach. PEX is a complex progressive multifactorial disorder of the ECM that manifests primarily in the ocular tissues. Impaired cross-linking of the ECM proteins and their subsequent aggregation is a hallmark of PEX. Association of numerous ECM proteins, such as elastin, tropoelastin, fibrillin-1, matrix metalloproteinase and their inhibitors, LOXL1 and FBLN5, with the disease pathology substantiates the debilitating effect of impaired ECM production and maintenance in development of PEX fibrils and subsequent PEX pathogenesis. FBLN5 is a matricellular scaffold protein with a crucial role to play in elastogenesis. FBLN5 interacts with integrins on the cell surface through its N-terminal domain and with LOXL1 through its C-terminal domain and brings the other ECM proteins into close proximity to facilitate elastogenesis [25]. Genetic variants in FBLN5 and its deregulation lead to various elastinopathies such as, cutis laxa, pelvic organ prolapse (POP), Charcot-Marie-Tooth disease and AMD [14,20,21,26]. Missense substitutions leading to improper secretion of FBLN5 and reduced interaction with elastin and fibrillin-1 have been reported in recessive cutis laxa [27]. Khadzhievaa et al. reported association of several tag SNPs with the advanced POP [22].
We and others have reported the dysregulation of FBLN5 in PEX patients which may result in impaired elastic fiber formation, degenerative tissue alterations, and subsequent ECM protein deposition [12,28]. We also identified novel genetic association of two variants, rs7149187 in the 5'-UTR and rs929608 in the 10th intron within FBLN5 with PEX which were, however, not found to be causal variants [12]. In the present study, we identified two intronic variants, NC 000014.9:g.91913280G>A (rs17732466:G>A) and NC 000014.9:g.91890855C>T (rs72705342:C>T) to be associated with PEXG as risk factors. However, these SNPs were not found to be associated with the early stage of PEX, PEXS. The minor alleles ' A' and 'T' at rs17732466 and rs72705342 were present in higher frequency in the controls and conferred a protective effect. The frequency of risk alleles at rs17732466 and rs72705342 was 78.0% and 82.0%, respectively, in PEXG. The 1000 Genomes Project data on Ensembl database recorded that the frequency of the binding at rs72705342 as predicted on PROMO software. The 'C' allele at rs72705342 is predicted to bind to GR-α and TFII I transcription factors, the binding of which is predicted to be lost on change in allele from 'C' to 'T'. (E) The EMSA shows that in the presence of TFII I specific antibody, the intensity of the shift reduced which did not happen in the presence of non-specific antibody (HSF1) taken as control. (F) The EMSA shows that in the presence of GR-α-specific antibody, the intensity of the shift reduced which did not happen in the presence of non-specific antibody taken as control. risk allele 'G' at rs17732466 is the highest in African population (88.0%), followed by American (87.0%), East Asian (84.0%), European (75.0%), and South Asian (73.0%) populations. The frequency of the risk allele 'C' at rs72705342 is the highest in African population (96.0%), followed by American (88.0%), East Asian (85.0%), European (77.0%) and South Asian (74.0%) populations. We found that the frequency of risk genotype 'CC' at rs72705342 was the highest in PEXG (64.0%) compared with 56.0% in control and 58.0% in PEXS. Further, the haplotypes among the Tag-SNPs across the gene showed significant association with PEXG but not with PEXS. This suggests that the underlying mechanism of pathogenesis of PEXS and PEXG could be different with novel risk factors contributing to the severity of PEX in its advanced stage. Many diseases show association of genetic variants with the severe forms of the disease compared with their early stages. Variants in CFH and ARMS2 were found to be more common with increasing severity of ARMD [29]. Tang et al. show that an intronic variant in the PAX6 gene is associated with extreme myopia but not with mild myopia [30]. Our findings suggest that these variants in FBLN5 might contribute to development and/or progression of PEX to a severe form rather than to the onset of the disease.
As causal variants for complex disorders are also found in regions outside the coding areas of protein coding genes, the functional effect of rs17732466 and rs72705342 was assessed [31,32]. Reporter assays showed that rs72705342 has an allele-specific regulatory effect on FBLN5 promoter. The risk allele 'C' at rs72705342 significantly reduced the reporter gene expression compared with allele 'T' . This finding is supported by the eQTL data from GTeX database which shows a significant increase in FBLN5 expression in tissues from individuals with the 'TT' genotype compared with those with 'CC' genotype (Supplementary Figure S3). Further, DNA-protein binding assay showed that the sequence with 'C' allele has more affinity to protein binding compared with that with the 'T' allele implying that the allele-specific differential regulation by rs72705342 could be due to differential transcription factor binding at the alleles. This finding was supported by the in silico analysis which showed differential binding of two transcription factors, GR-α and TFII I with the 'C' allele at rs72705342. Electrophoretic mobility shift assays including the antibodies against TFII I and GR-α showed that both these proteins might bind to rs72705342 C/T. However, the extent of binding affinity of these proteins to the alleles at rs72705342 needs to be verified further through the use of CRISPR-edited cells lines followed by chromatin immunoprecipitation experiments.
This study might have a limitation in terms of the sample size included and replication of this genetic association study at a larger scale and in different cohorts is needed to gain confidence in the association of these intronic variants with PEXG and their non-association with PEXS. Also, though rs72705342 showed regulatory effect on gene expression, we observed that FBLN5 is down-regulated in the lens capsule of PEXS individuals but not PEXG. However, the association of rs72705342 with only PEXG but not PEXS suggests that this SNP might have an unknown effect on PEXG progression. It is possible that this variant could influence the expression of distal genes contributing to PEXG pathology that needs to be studied further. In conclusion, the present study elucidated genetic association of 13 tag SNPs in FBLN5 with PEX. We identified several SNPs and haplotypes of FBLN5 to be associated with the advanced stage of PEX, PEXG but not with the early stage of PEXS. The intronic SNP rs72705342 showed a plausible regulatory effect on FBLN5 expression. Further, in vivo studies can help to understand the exact effect of these deep intronic variants and haplotypes on the progression of the disorder.

Data Availability
The data pertaining to this study are within the published article and its supplementary files. Any additional data will be made available by the corresponding author on reasonable request.

Ethics Approval
This study was approved by the Institutional Biosafety and Human Ethics committee of National Institute of Science Education and Research and adheres to the tenets of Declaration of Helsinki. Written informed consent has been obtained from all the participants for the study.