Long Non-Coding RNAs Gene Variants as Molecular Markers for Diabetic Retinopathy Risk and Response to Anti-VEGF Therapy

Introduction

Diabetes mellitus (DM) is a devastating health problem worldwide with an estimated global prevalence of 9.3% (463 million people) in 2019, predicted to rise to 10.2% (578 million) by 2030, according to the recent “International Diabetes Federation” report.¹ Diabetic retinopathy (DR) is one of the main microvasculopathy associated with long-term diabetes. Retinal inflammation/neovascularization, vascular hyperpermeability, and cell apoptosis play major roles in DR etiopathology.² Healthy endothelial cells are the basis of normal blood vessels, whereas endothelial cell dysfunction is a risk indicator of diabetic angiopathy.³ Accumulating evidence indicated that diabetes-related microvascular complications result from genetic and environmental interactions.^4,5

Long non-coding RNAs (lncRNAs), which are >200 nucleotides in length, are emerging regulatory non-coding RNAs, participate in epigenetic and transcriptional/post-transcriptional regulation of several cellular pathways in different diseases, including DM and its complications.^6–9 Several lncRNAs have been implicated in DR,^10–12 and their deregulation was proved to be associated with DR susceptibility and/or response to treatment^13,14 (Table 1).

Table 1 The Role of lncRNAs in Diabetes and Diabetic Complications

The lncRNA TUG 1 “Taurine UP-regulated 1” was reported as one of the identified genes that are up-regulated in response to amino acid taurine, which induces rod photoreceptor biogenesis.^15,16 This type of lncRNA is essential for photoreceptors biogenesis in the developing rodent retina.¹⁶ Knockdown of TUG1 resulted in a defect of migration of the developing rod photoreceptors into the outer nuclear layer and increased transfected cell apoptosis.¹⁷ Furthermore, TUG1 upregulation stimulates Wnt pathway-dependent proliferation and migration of endothelial cells, thereby promoting the occurrence and progression of diabetic atherosclerosis.¹⁸

Another lncRNA reported to be implicated in DR pathophysiology is MIAT “Myocardial Infarction Associated Transcript”, also known as Gomafu, or RNCR2 “retinal non-coding RNA 2,” which is expressed in cardiomyocytes and the nucleus of multiple retinal cells.¹⁹ It was found to be up-regulated in the retina of diabetic rats and the fibrovascular membrane of patients with diabetes induced by cellular hyperglycemia.²⁰

The highly conserved lncRNA MALAT1 “metastasis-associated lung adenocarcinoma transcript 1”, also named NEAT2 “nuclear-enriched abundant transcript 2”,²¹ showed significant overexpression in hyperglycemic cell model (the RF/6A), in samples from the aqueous humor, and fibrovascular membranes of patients with diabetes.¹¹ It can regulate retinal endothelial cell pathophysiology and microvascular growth under a hyperglycemic milieu via several cellular pathways, including the “p38/MAPK (mitogen-activated protein kinase)” signaling pathway.²² Silencing of MALAT1 significantly mitigates diabetes-induced retinal neovascularization, vascular hyperpermeability, and retinitis.²³

Lastly, the cytoplasmic lncRNA-SENCR “Smooth muscle and Endothelial cell-enriched migration/differentiation-associated long non-coding RNA” was found to be a vascular cell-enriched lncRNA that can regulate FoxO1 “fork-head box protein O1” and TRPC6 “Transient Receptor Potential Cation channel 6”, hence promoting the proliferation and migration of smooth muscle cells.²⁴ SENCR has been implicated in regulating pluripotent cells endothelial differentiation and human vascular endothelial cell angiogenic capacity.²⁵ Its silencing downregulates myocardin and smooth muscle contractile genes and up-regulates the cell migration-related genes.²⁶

Several studies have reported that lncRNA-related single nucleotide polymorphisms (SNPs) may impact susceptibility to DR.^27,28 However, the TUG1 (rs7284767G/A), MIAT (rs1061540T/C), MALAT1 (rs3200401C/T), and SENCR (rs12420823C/T) variants which are selected based on specific criteria detailed in 2.4 section of this study, were not investigated yet in association with DR risk or response to “anti-vascular endothelial growth factor (VEGF)” aflibercept treatment. It is worth noting that anti-VEGF agents used in clinical practice, such as ranibizumab and aflibercept, are considerably different in terms of molecular interactions when they bind with VEGF.²⁹ In this sense, this study aimed to investigate whether the selected genetic variants of the specified lncRNAs are associated with DR susceptibility, clinic-laboratory data, or the short-term (after four weeks) response to aflibercept therapy in the hope to find new biomolecular markers that help in patient risk stratification and/or listed within the pharmacogenetics-related variants.

Patients and Methods

Study Population

This cross-sectional study recruited 126 consecutive adult patients with type 2 DM attending the “Ophthalmology Department, Suez Canal University Hospitals” and a private Clinic, Ismailia, Egypt. The “Early Treatment Diabetic Retinopathy Study” (ETDRS) report³⁰ was followed in patients’ evaluation and subclassification into DR group (n = 73) and non-DR group (n = 53). About 42.5% of DR and 22% of non-DR patients were on insulin therapy. Both groups had a comparable period of diabetes duration. Patients presented with other ophthalmic disorders, including neovascular glaucoma, age-related macular degeneration, or hematogenous retinal detachment, vasculopathy other than DR, history of ocular trauma or surgery, having other chronic diseases and/or severe comorbidities were excluded. The “institutional Research Ethics Committee of Faculty of Medicine, Suez Canal University” approved the study (approval No.4465). This study was executed according to the Declaration of Helsinki principles and its amendments. Written informed consent was obtained from the study population before taking part.

Clinical Assessment

A complete ophthalmic evaluation was done for all participants. Examination at the first visit to the ophthalmic clinic and throughout the follow-up visits included (1) the BCVA “best-corrected visual acuity” using logMAR “log of the Minimum Angle of Resolution”, (2) the anterior/posterior segment examination using slit-lamp biomicroscopy and +20 D/+90 D lenses, respectively, (3) colored fundus photography and (4) optical coherence tomography (OCT). Based on the initial FFA “fundus fluorescein angiography” assessment, maculopathy type and subclassification (focal vs diffuse) have been assigned. The ischemic type of DR was ruled out to avoid any bias in the results. The “non-proliferative DR (NPDR)” (subclassified into mild, moderate, or severe subtype) and “proliferative DR (PDR)” diagnoses were reached based on the ETDRS. Based on the OCT and the clinical findings, the diabetic macular edema (DME) diagnosis was assigned.³⁰

Optical Coherence Tomography Scans

Patients retinas were subjected to OCT scans as explained previously.³¹ Patients presented with DME underwent an intravitreal injection (IV) of aflibercept (2 mg in 0.05 mL), “Eylea 40 mg/mL, Bayer Pharma AG, Berlin, Germany”, session after blood sampling. All the specified percussions and standards were followed pre-, during, and post-injection, as shown previously.¹⁴ Following-up was scheduled on the first day/one week, and one month after the first injection session. Both BCVA and OCT were reperformed after one month to evaluate the CMT. The initial response to aflibercept (one dose) was assigned as

BCVA improvement more than two lines of the Snellen s chart (Converted to Log MAR units for the statistical analysis), and CST reduction more than 15% of the pre-treatment thickness.¹⁴

The repeated injection was arranged when appropriate. Both “BCVA” and “CMT” changes were updated and recorded during the follow-up visits for all treated patients.

Blood Sampling and lncRNA Variants Genotyping

A total of five milliliters of venous blood samples were withdrawn from all participants after overnight fasting under aseptic conditions on EDTA tubes. Buffy coat genomic DNA was isolated using “QIAamp DNA Blood Mini kit (Catalog # 51104; Qiagen GmbH)” according to the vendor’s guidelines. The nucleic acid concentration and purity were evaluated using “NanoDrop ND-1000 (NanoDrop Technologies, Inc. Wilmington, DE, USA)”.

The studied variants selection was based on (1) dbSNP (www.ncbi.nlm.nih.gov) search for the minor allele frequency (MAF) ≥0.1 in the selected lncRNAs to get adequate statistical power, (2) mining in previous literature which showed evidence of the functional significance of these polymorphisms, and/or (3) no previous literature relates these SNPs with DR at least in our population. Genotyping was done on “StepOne Real-Time Polymerase Chain Reaction system (Applied Biosystems, USA)”. Polymerase chain reaction (PCR) was run in duplicates in a total volume of 25 µL containing 20 ng of extracted genomic DNA diluted to 11.25 µL with nuclease-free water, 12.5µL TaqMan genotyping PCR Master Mix, and 1.25 µL predesigned primer/probe sets TaqMan SNP Genotyping Assay (20x) working stock (assay IDs: C___2566592_10, C___2467719_1_, C___3246069_10, and C__11783392_10), for TUG1 (rs7284767G/A), MIAT (rs1061540T/C), MALAT1 (rs3200401C/T), and SENCR (rs12420823C/T), respectively. “No-template negative controls” (NTCs) were applied in each PCR run to confirm free reaction contamination. The PCR program was set as followed: 10 min (95°C) for one cycle, followed by 15 sec (90°C) and 1 min (60°C) for 40 two-stage cycles. The genotyping call was done by Applied Biosystems software with a 99.2% call rate and 100% concordance rate.

Statistical Analysis

The SPSS for Windows (version 27) and R (version 3.5.3) were used for the statistical analysis. G*Power (version 3.1.9.2.) was applied for sample size calculation. At a 95% significance level (alpha 0.05) and an effect size of 0.37 with a minimal sample size required to reject the null hypothesis (n= 126), the calculated study power was 91%. Data distribution and normality were checked by the Kolmogorov–Smirnov test. “Hardy-Weinberg equilibrium” (HWE) and allele/genotype frequencies of the studied variants were estimated within each group, as mentioned previously.³² Multivariable regression models were employed. Adjusted odds ratio (OR) and 95% confidence interval (CI) for multiple genetic association models were calculated using SNPStats (https://www.snpstats.net/start.htm). One-way ANOVA test was applied for parametric attributes, and Mann Whitney or Kruskal–Wallis tests were executed for non-parametric data. A Chi-square test was applied for qualitative variables. Significance was set at P<0.05.

Results

Baseline Characteristics of the Study Population

A total of 126 diabetic cohorts were included in the current analysis. Of these, 73 patients who developed diabetic retinopathy were compared to those who did not progress to retinopathy phenotype. Their mean ages were 59.8 ± 9.5 years and 62.5 ± 7.5 years, respectively (P = 0.09). Females accounted for 69.9% in DR cohorts and 79.2% in DM group (P = 0.30). Hypertension comorbidity was found in 71.2% and 62.3% of patients, with no significant difference between the group (P = 0.33). In contrast, patients with DR had more prevalent insulin intake (41.1% versus 17.0%, P = 0.006), and higher glycosylated hemoglobin (HbA1c) levels (8.92% ± 2.2% vs 5.45% ± 1.01%, P = 0.025) (Table 2). These variables were taken into consideration in the downstream multivariable regression analysis.

Table 2 Baseline Characteristics of the Study Population

Clinical Assessment of Diabetic Retinopathy Patients

DR patients were categorized according to the disease stage into mild, moderate, and severe non-proliferative and proliferative diseases. As depicted in Table 3, patients with PDR were significantly older (65.8 ± 10.9 years, p = 0.018), more likely to have prolonged disease duration (26.4 ± 7.1 years, P <0.001), and higher frequency of insulin intake (76.9%, P = 0.001). There were no significant inter-group differences regarding their therapeutic response (P = 0.89) nor their best-corrected visual acuity (P = 0.17).

Table 3 Clinical and Ophthalmologic Assessment of Patients with Diabetic Retinopathy of Different Grades

Allelic Discrimination Analysis

Apart from TUG1 polymorphism (p <0.001), genotype frequencies in diabetic controls were in accordance with Hardy–Weinberg equilibrium (MIAT: p = 0.29, MALAT1: P = 0.79, SENCR: P = 0.78). For TUG1 (rs7284767), MIAT (rs1061540), MALAT1 (rs3200401), and SENCR (rs12420823), overall minor allele frequencies were 0.37 (A), 0.26 (C), 0.47 (T), and 0.42 (T), respectively. On comparison between the two study groups, MIAT C variant was more frequent among DR patients (C allele: 44% in DR group versus 26% in DM group, P = 0.004) (Figure 1A). Similarly, higher proportions of MIATT/C and C/C genotypes were observed in DR group (C/C: 27% versus 13%, C/T: 33% versus 26%, P = 0.049) (Figure 1B). Comparison with other ethnic populations from 1000Genome Project is shown in Figure 1C.

Figure 1 Genotype and allele frequencies of the studied genetic variants. (A) Genotype frequencies of polymorphisms. (B) Allele frequencies of polymorphisms. A Chi-square test was applied. Statistical significance was set at P < 0.05. (C) Allele frequencies of TUG1 (rs7284767), MIAT (rs1061540), MALAT1 (rs3200401), and SENCR (rs12420823) in 1000Genome Project. Bold values indicate significance at P-value < 0.05.

Association of lncRNA Variants with Disease Risk

As seen in Table 4, TUG1 A/G was 3 times more likely to develop DR under heterozygote comparison (OR = 3.15, 95% CI = 1.15 to 8.64) and over-dominant model (OR = 3.33, 95% CI = 1.28 to 8.67). MIAT T/C and C/C were associated with three and eight times more odds for developing disease [heterozygote model: OR = 2.97, 95% CI = 1.11 to 7.93, homozygote model: OR = 8.43, 95% CI = 2.40 to 29.59, dominant model: OR = 4.31, 95% CI = 1.78 to 10.47, recessive model: OR = 5.23, 95% CI = 1.66 −16.45, and log-additive model: OR = 2.92, 95% CI = 1.60 to 5.33]. In contrast, MALAT1 C/T rendered protection with 60% decreased susceptibility to develop DR disease under heterozygote comparison (OR = 0.40, 95% CI = 0.16 to 0.99) and dominant model (OR = 0.41, 95% CI = 0.18 to 0.96). Gene–gene interaction analysis revealed that carriers for GTCT and GCCC genotype combinations had 17- and 20-times higher disease risk, respectively (P = 0.012) (Table 5).

Table 4 Genetic Association Models for the Study Long Non-Coding RNAs and Disease Risk

Table 5 Combined Genotype Association with Disease Risk

Association of lncRNA Variants with Disease Severity

SENCR T variant was associated with better pre-treatment BCVA level (0.63 ± 0.19 in C/C compared to 0.45 ± 0.27 in C/T and 0.46 ± 0.26 in T/T carriers, P = 0.021). In contrast, MIAT T/T homozygosity was associated with higher grade (PDR: 33.3% in T/T compared to 10% in C/C and 4.2% in C/T carriers, P = 0.012) (Table 6).

Table 6 Association Between Genotypes and Patient Characteristics

Association of lncRNA Variants with Drug Response After Anti-VEGF Treatment

TUG1 A and MIAT C alleles were associated with non-early response to aflibercept treatment (P = 0.020 and 0.014, respectively). Carriers for TUG1 A allele showed the lowest median change in CMT (−11.5 in A/A compared to −90 in A/G and −54 in G/G, P = 0.016) and BCVA (0 in A/A compared to −0.15 in A/G and −0.1 in G/G, P = 0.005) (Table 6). The same allele (TUG1 A) was significantly associated with treatment failure under homozygote comparison (OR = 25.10, 95% CI = 2.53 to 249.44), dominant model (OR = 8.18, 95% CI = 1.48 to 45.29), recessive model (OR = 8.81, 95% CI = 1.49–52.06), and log-additive model (OR = 5.02, 95% CI = 1.60 to 15.76). In addition, MIAT T/C genotype was associated with poor therapeutic response under heterozygote comparison (OR = 10.23, 95% CI: 1.51 to 69.15), dominant model (OR = 5.95, 95% CI = 1.04 to 34.23), and over-dominant model (OR: 6.26, 95% CI: 1.41–27.80) (Table 7).

Table 7 Association Between Genotypes and Drug Response Phenotype

Discussion

The current treatment strategies for DR are mainly used in the advanced stages of the disease and may be related to several side effects. In this sense, it is necessary to explore a new landscape of susceptibility genetic variants associated with disease development, severity, or response to treatment to identify the individuals most likely to develop DR and predicting treatment effects. To this end, we genotyped for the first time four variants belonging to lncRNAs which proved previously to be implicated in one or more pathophysiological mechanisms of diabetes or its complications (Table 1). Our results show a variable degree of disease risk, severity, and/or early response to aflibercept treatment.

We found that carriers of lncRNA TUG1 A/G of rs7284767 variant were three times more likely to develop DR under heterozygote and over-dominant models. Also, carriers of the TUG1 A allele under several genetic models were associated with non-early response to aflibercept treatment as they showed the lowest median change in CMT and BCVA.

The rs7284767 variant was found to be located in one of “novel biologically relevant; NBR” target sites on the TUG1 and was associated with unusual patterns of long-range haplotype conservation in the human genome³³ and exhibited significantly high integrated haplotype score value; “a gnomically standardized measure for conservation of long-range haplotype associated with a given SNP in a population”.³⁴ This variant was predicted to disrupt the NBR target site for miR-20 in humans.³³ Interestingly, miR-20 was associated with diabetes and DR in previous studies^35,36 and proved to have diagnostic and prognostic roles in patients with DR.³⁶ Given the targeting role of this microRNA for the 3′-untranslated regions of hypoxia‐inducible factor‐1α and VEGFA genes, it is not surprising to be associated with angiogenesis-related disorders, including DR and with treatment response to anti-VEGF aflibercept treatment. A recent study by Duan et al has identified the association of rs7284767 allelic variation with increased plasma TUG1 levels in patients with knee osteoarthritis.³⁷ This finding and all the mechanisms mentioned above could prove that the rs7284767 variant may be functional, supporting its association with disease risk and response to anti-VEGF treatment in the present study population. Future functional assays and detailed genetic analyses will be required to determine the detailed biological significance of this variant and other studied ones in the future.

The present work also explored that MIAT rs1061540T (C/T) and (C/C) genotypes were associated with three and eight times more odds for developing DR, and T/T homozygosity was associated with a higher grade of PDR relative to other genotypes. Also, the C/T genotype was associated with poor therapeutic response under heterozygote, dominant, and over-dominant models. Although accumulating evidence confirmed the implication of MIAT in every stage of DR by several genetic/epigenetic mechanisms and its knockdown has a therapeutic advantage in neovascular-related disorders as DR,^{20,27,38–40} the specific impact of the rs1061540 variant on MIAT expression and association with DR is yet to be identified. In this sense, more studies are recommended to explore the precise molecular mechanism by which this variant could impact disease susceptibility/severity and/or treatment response.

In contrast to the SNPs mentioned above, the lncRNA MALAT1 rs3200401 (T/C) genotype in the current study population conferred protection with 60% decreased susceptibility to develop DR under heterozygote comparison and dominant models. Several in vivo and in vitro studies unraveled the association between MALAT1 level and diabetes/hyperglycemia. It was found to be up-regulated in diabetic mice retinas,¹¹ regulating the retinal endothelial cells function in terms of cell proliferation/migration and tube biogenesis, and the pathological growth of the retinal microvasculature under hyperglycemic conditions.^23,41 Its pharmacological inhibition can reduce the retinal endothelial cell proliferation phenotype, vasculature growth, inflammatory mediators, and the pro-inflammatory cytokines such as serum amyloid antigen 3, tumor necrosis factor-α, and interleukin-6.^39,42 It is worth noting that the lncRNA variants may exert their impact through lncRNA transcript alternative splicing or change in its secondary structure, resulting in loss or gain of function.⁴³ Interestingly, the rs3200401 variant is located in the binding site of “serine/arginine-rich splicing factor 2 (SRSF2)” that controls pre-mRNA alternative splicing of several targets,⁴⁴ including the VEGFA/VEGFR genes, which are considered the main stimulatory signal of angiogenesis in vivo.^45,46 Using the “lncRNASNP database” [10.1093/nar/gkx1004], Wang et al predicted that the rs3200401 variant could cause “1.62 kcal/mol minimal free energy change” which change MALAT1 structure, weaken the MALAT1/SRSF2 interaction and alter miRNA-MALAT1 binding (eg, impact binding to miR-3661 and miR-1324).⁴⁷ In this sense, it is biologically plausible to speculate that deregulation of the mentioned cellular mechanistic due to MALAT1 (rs3200401C/T) polymorphism may impact the “molecular sponging” function and stability of MALAT1 that result in a less pathophysiological derangement, mainly the angiogenesis process and its progress with lower susceptibility to DR.

Regarding the SENCR (rs12420823C/T) intronic variant, the present study showed that it was associated only with better pre-treatment best-corrected visual acuity levels assessed by the logarithm of the minimum angle of resolution. This vascular-enriched lncRNA gene has been found to overlap the “Friend Leukemia Integration virus 1; FLI1” gene, an essential regulator of endothelial development, and the lncRNA SENCR can contribute to this regulation and induce proliferation, migration, and angiogenesis of human umbilical endothelial cell.²⁵ As there were no previous publications that uncover the role of the studied SENCR rs12420823 in general, more future studies are warranted to confirm the relation of this variant with angiogenesis-related diseases, including DR.

Our gene–gene interaction analysis revealed that carriers for TUG1, MALAT1, MIAT, and SENCR (GTCT) and (GCCC) genotype combinations had 17- and 20-times higher disease risk, respectively. This finding is consistent with that the concurrent impact of multiple gene variants is strong enough to unleash the genetic association with disease susceptibility and/or response to treatment.⁴⁸ This multi-polymorphic model could be a potentially useful tool in patient risk stratification and future personalized medicine.

Although this study, up to the authors’ knowledge, is the first study to uncover the impact of four angiolncRNAs polymorphisms on DR susceptibility and early response to aflibercept, some limitations should be considered. First, the study cohort is a hospital-based population with difficulty in avoiding selection bias. Second, the limited sample size available in our hospital in the study period with an application of very stringent exclusion criteria. Third, the limited selected study variants with MAF ≥ 0.10 to be included in this study to achieve adequate statistical power. Hence, it is recommended to replicate the work in multi-center, larger-scale studies in different ethnicity populations to validate the current findings. Moreover, including other variants related to the studied lncRNAs, supported with functional analysis to uncover their molecular mechanisms underlying the disease susceptibility and/or response to treatment, will clarify the complete picture.

Conclusion

This study revealed that the genetic variation of lncRNAs TUG1 (rs7284767G/A) and MIAT (rs1061540T/C) could be molecular determinants for increased DR susceptibility and early response to aflibercept treatment. Otherwise, the MALAT1 rs3200401 (T/C) variant conferred DR protection. Also, the multi-polymorphic models (GTCT and GCCC) of the TUG1, MALAT1, MIAT, and SENCR, respectively, were associated with increased DR risk. Further large population-based and functional in vivo and in vitro studies are required to confirm the significant association of the studied variants with DR susceptibility and/or response to aflibercept treatment, particularly in other ethnicities, and to explore the molecular mechanism(s) by which these variants confer DR risk and/or response to anti-VEGF therapy.