Further evaluation of differential expression of keratoconus candidate genes in human corneas

Background Keratoconus (KTCN) is a progressive eye disease, characterized by changes in the shape and thickness of the cornea that results in loss of visual acuity. While numerous KTCN candidate genes have been identified, the genetic etiology of the disease remains undetermined. To further investigate and verify the contribution of particular genetic factors to KTCN, we assessed 45 candidate genes previously indicated as involved in KTCN etiology based on transcriptomic and genomic data. Methods The RealTime ready Custom Panel, covering 45 KTCN candidate genes and two reference transcripts, has been designed. Then, the expression profiles have been assessed using the RT-qPCR assay in six KTCN and six non-KTCN human corneas, obtained from individuals undergoing a penetrating keratoplasty procedure. Results In total, 35 genes exhibiting differential expression between KTCN and non-KTCN corneas have been identified. Among these genes were ones linked to the extracellular matrix formation, including collagen synthesis or the TGF-β, Hippo, and Wnt signaling pathways. The most downregulated transcripts in KTCN corneas were CTGF, TGFB3, ZNF469, COL5A2, SMAD7, and SPARC, while TGFBI and SLC4A11 were the most upregulated ones. Hierarchical clustering of expression profiles demonstrated almost clear separation between KTCN and non-KTCN corneas. The gene expression levels determined using RT-qPCR showed a strong correlation with previous RNA sequencing (RNA-Seq) results. Conclusions A strong correlation between RT-qPCR and earlier RNA-Seq data confirms the possible involvement of genes from collagen synthesis and the TGF-β, Hippo, and Wnt signaling pathways in KTCN etiology. Our data also revealed altered expression of several genes, such as LOX, SPARC, and ZNF469, in which single nucleotide variants have been frequently identified in KTCN. These findings further highlight the heterogeneous nature of KTCN.


INTRODUCTION
Keratoconus (KTCN) is a progressive eye disorder characterized by thinning and conically shaped protrusion of the cornea, resulting in severe visual impairment (Rabinowitz, 1998). The first symptoms of KTCN usually appear during the second or third decade of life (Rabinowitz, 1998). However, individuals with early-onset KTCN have also been reported (Léoni-Mesplié et al., 2012). The estimated prevalence of KTCN varies between 1 in 375 and 1 in 2,000 individuals and depends on geographic location and ethnic origin (Godefrooij et al., 2016).
Keratoconus is considered a multifactorial disease with a substantial genetic contribution, supported by evidence from twin, family, and population studies (Tuft et al., 2012;Naderan et al., 2016;Lu et al., 2013). Because of the complex etiology and the genetic heterogeneity of the disease, the identification of specific risk factors for KTCN is difficult. However, the use of various approaches has already allowed detection of many candidate genes that could be implicated in KTCN as have been previously reviewed Abu-Amero, Al-Muammar & Kondkar, 2014;Bykhovskaya, Margines & Rabinowitz, 2016;Karolak & Gajecka, 2017;Valgaeren, Koppen & Van Camp, 2018;Loukovitis et al., 2018;Lucas & Burdon, 2020).
Gene expression studies are also a valuable source of information about functional genes and pathways important for KTCN pathogenesis. A high-throughput RNA sequencing (RNA-Seq) showed a serious alteration of numerous genes in KTCN corneas (Kabza et al., 2017). The most significant downregulation was observed among genes belonging to the TGF-β, Hippo, and Wnt signaling and collagen synthesis and maturation networks that are responsible for proper corneal organization and regulation of corneal ECM remodeling (Kabza et al., 2017). Abnormalities within genes involved in ECM have also been reported in subsequent RNA-and DNA-based studies (Khaled et al., 2018;You et al., 2018;Sharif et al., 2019;. Thus, disruptions within these molecular cascades could be potentially responsible for corneal changes underlying the development of KTCN.
While the genes mentioned above have been hypothesized to play a role in KTCN, some of these identified factors account for a limited number of individuals only being found in single families or particular studied populations. Also, in some cases, the results of the initial findings are contradictory. Thus, in this study, we have taken advantage of the previous RNA-and DNA-based KTCN analyses to select 45 candidate genes for further investigation aiming to validate their contribution to KTCN. The possible involvement of these genes in disease etiology was re-evaluated using the RealTime ready Custom Panels in the RT-qPCR experiments, performed in KTCN and non-KTCN corneas.

Subjects
Six KTCN patients and six non-KTCN individuals were enrolled in this study and underwent a complete ophthalmic evaluation. The KTCN diagnosis was made based on the previously described criteria (Karolak et al., 2016b;Kabza et al., 2017). The clinical characteristics of all enrolled individuals are shown in Table S1. All individuals provided written informed consent for participation in the study, in accordance with the Declaration of Helsinki. The research protocol was approved by the Institutional Review Board at Poznan University of Medical Sciences (453/14 and 755/19).

Candidate gene assay design
Forty-five genes were selected for further evaluation based on both corneal RNA-Seq results and previous genomic data. The first KTCN candidate gene set included 26 genes that have been differentially expressed in our earlier studies performed in KTCN and non-KTCN corneas. These genes were selected from the top molecular pathways overrepresented across deregulated genes and encode core elements of collagen synthesis and maturation pathways, the TGF-β, Hippo, and Wnt signaling pathways, as well as their potential regulators (Kabza et al., 2017). The second set consists of 19 genes that have been previously reported as involved in KTCN based on the detection of putative variants within the gene, localization within the linkage region, or proximity to variants associated with KTCN. The selected genes are involved in inflammation, intracellular signaling, ubiquitination, and other processes reported as crucial for the functioning of the eye or disease development. The complete list of studied genes is presented in Table S2.

RT-qPCR analysis
To analyze gene expression profiles, total RNA samples from six KTCN to six non-KTCN corneal tissues were reverse transcribed to cDNA with the Transcriptor First Strand cDNA Synthesis Kit (Roche Diagnostics, Penzberg, Germany), according to the manufacturer's procedure. Expression levels of 45 candidate KTCN genes and two housekeeping transcripts (Table S2) were assessed using the RealTime ready Custom Panel 96 (Roche Diagnostics, Penzberg, Germany), containing pre-plated primer pairs and probes. RT-qPCR analyses were conducted using the LightCycler 96 System (Roche Diagnostics, Penzberg, Germany) in a total volume of 20 ml containing the FastStart Essential DNA Probe Master (Roche Diagnostics, Penzberg, Germany), and 0.5 ng of cDNA. Each reaction was performed in triplicate.

RT-qPCR data analysis
Relative quantification of the gene expression was normalized to the level of the GAPDH and IPO8 transcripts with the comparative C T method. The R environment was used to analyze differential expression and to generate a heat map showing hierarchical clustering (Ward linkage) of KTCN and non-KTCN samples based on the expression profiles (R Development Core Team, 2008). The statistical significance of differential gene expression was evaluated using Student's t-test. The adjusted P-values were calculated using the Benjamini-Hochberg (FDR) procedure.

Comparison of RT-qPCR data with RNA-Seq data
The log2 transformed fold change (FC) values of gene expression levels between KTCN and non-KTCN samples were calculated for RT-qPCR and previously reported RNA-Seq data, GSE77938 (Kabza et al., 2017). To calculate the log2 FC values based on RNA-Seq results, gene expression data from both discovery and replication RNA-Seq experiments were analyzed together. Briefly, Illumina adapter sequences, poor quality regions (average Phred quality score below five), and sequences matching human rRNAs were removed from the sequenced short reads using the BBDuk2 program from the BBTools package (http://jgi.doe.gov/data-and-tools/bbtools). Filtered reads shorter than 50 bp were discarded. The expression values of known genes and transcripts from GENCODE 25 annotations (Ensembl 87) were estimated using Salmon with sequence-specific and GC content bias correction enabled. RNA-Seq FC values were obtained from limma package using a previously published protocol (Ritchie et al., 2015;Patro, Duggal & Kingsford, 2015;Law et al., 2016;Kabza et al., 2017). Pearson correlations between log2 FC values for selected genes obtained from RT-qPCR and RNA-Seq data were evaluated using the R package (R Development Core Team, 2008).
The analysis of RT-qPCR data showed that 23 genes had decreased expression level in the KTCN patients compared to the non-KTCN individuals, including six genes (CTGF, ZNF469, TGFB3, COL5A2, SMAD7, and SPARC) with log2 FC values lower than −2.0 and six genes (LOX, TEAD2, TEAD4, TGFB1, TGFBR2, and ZEB1) with log2 FC values between −1.0 and −1.99 (Table 1). Twelve genes presented with increased expression, including SLC4A11 and TGFBI with log2 FC values greater than 1.0 (Table 1). No difference between the level of CTNNB1 transcripts was observed for KTCN and non-KTCN individuals. Relative gene expression is presented in Table S3.
Sample clustering performed based on all differentially expressed genes showed almost clear separation of KTCN and non-KTCN clusters. Only one non-KTCN sample (KR49) was misclassified in the KTCN cluster ( Fig. 1).

Comparison of RT-qPCR data with RNA-seq data
The log2 FC values for 36 genes obtained during the RT-qPCR and RNA-Seq showed a strong correlation between the two sets of results. The Pearson correlation coefficient (r) was 0.91 (Fig. 2).

DISCUSSION
Differences in transcriptomic profiles between various cells and tissues play an essential role in the normal functioning of the human body (Frith, Pheasant & Mattick, 2005;Cummings et al., 2017). During the last few years, numerous transcriptome studies have been performed to understand how the changes in gene expression may influence different phenotypes, including KTCN (Nielsen et al., 2003;Macé et al., 2011;Joseph, Srivastava & Pfister, 2016;Kabza et al., 2017). The results of RNA-Seq-based transcriptomic profiling of human corneas suggested that downregulation of collagens and other genes encoding ECM components might be involved in KTCN etiology as it can lead to the reduction of the mechanical stability of affected tissues (Kabza et al., 2017). Moreover, it was postulated that the previously observed decreased level of collagen genes in KTCN corneas might be a result of downregulation of connective tissue growth factor (CTGF) and its modulators, which are involved in the TGF-β, Hippo, and Wnt signaling pathways (Kabza et al., 2017). Interestingly, changed expression patterns of the core elements of these three biological pathways in KTCN corneas were confirmed in the present RT-qPCR study. Three transforming growth factor beta (TGF-β) isoforms, TGFB1, TGFB2, and TGFB3, through binding to their receptors (TGFBR1, TGFBR2), initiate signaling that regulates the assembly of ECM. While TGFB1 and TGFB2 act as stimulators of a profibrotic response in the injured cornea, TGFB3 is known for its antifibrotic effect (Priyadarsini et al., 2015). TGF-β signaling stimulates the downstream synthesis of CTGF. As a result, corneal cells are activated and begin to produce different types of collagen and other ECM components (Blalock et al., 2003). The levels of TGFB1, TGFB3, TGFBR2, COL5A2, and CTGF were decreased in this study confirming that previously observed low expression of genes encoding major ECM components, including collagens, may be an effect of TGF-β signaling alteration in KTCN (Joseph, Srivastava & Pfister, 2016;Kabza et al., 2017). The TGF-β pathway is controlled by SMAD proteins, including SMAD family member 7, which through binding to TGF-β receptors, leads to an inhibition in TGF-β signaling (Yan, Liu & Chen, 2009). Interestingly, TGF-β transcriptional responses could also be regulated by non-Smad pathways (Zhang, 2009). That could explain decreased SMAD7 level and increased expression of TGFBI, a downstream effector of TGF-β, at the same time. The TGFBI is a significant component of human corneal stroma involved in cell adhesion and migration, which is induced by its interaction with several ECM elements, including collagens (Runager, Enghild & Klintworth, 2008). TGFBI variants have been frequently identified in patients with corneal dystrophies (Chao-Shern et al., 2019), as well as KTCN (Guan et al., 2012). While a significant upregulation of TGFBI in KTCN patients compared to control individuals was also detected previously (Bykhovskaya et al., 2016a), in most transcriptomic or proteomic studies the level of TGFBI in KTCN tissues was decreased (Takács et al., 1999;Zhao et al., 2002). Since different material was used in these mentioned analyses (whole corneas, corneal buttons, or corneal stroma), the observed results could vary because of the influence of other cells on the tissue expression. More studies are needed to verify this data.
It is known that TGF-β is integrated into higher-order networks and crosstalk between TGF-β, Wnt and Hippo signaling regulates the outcome of signaling activity (Attisano & Wrana, 2013). The Hippo molecular cascade is involved in controlling the eye size, and TEA domain transcription factors, encoded by the TEAD2 and TEAD4 genes, play an essential role in its regulation (Yu & Guan, 2013). Moreover, TEAD2 and TEAD4 transcription factors integrate Hippo and Wnt pathways at the nucleus through cooperation with YAP and TAZ molecules (Attisano & Wrana, 2013). Decreased level of TEAD2 and TEAD4, observed in this study in KTCN corneas, supports our previous data. Also, it indicates that abnormal expression of particular elements of the TGF-β, Hippo, and Wnt pathways might alter signaling crosstalk between these cascades in KTCN. Of note, our recent ES findings showed the accumulation of variants in several genes from Wnt signaling and/or focal adhesion pathways, deregulated in KTCN, further supporting possible involvement of these genes in disease etiology .
Apart from genes involved in the pathways mentioned above, selected as candidates for KTCN based on solely transcriptome study, there are numerous putative KTCN genes, which are revealed using linkage analyses, GWAS, or ES (Karolak & Gajecka, 2017). The level of these genes expression was also evaluated using RNA-Seq (Kabza et al., 2017). However, their significance in KTCN could not be clearly defined due to contradictory findings in genomic studies. Thus, in this study, 19 candidate KTCN genes were selected for re-evaluation using RT-qPCR. In KTCN corneas compared to non-KTCN corneas, we observed altered expression profiles of a few of them, including ZNF469, SPARC, and LOX.
The ZNF469 gene was initially associated with central corneal thickness, which is abnormal in KTCN or corneal dystrophies (Lu et al., 2010(Lu et al., , 2013Vitart et al., 2010) and now is one of the most discussed KTCN genes. Another KTCN candidate, the SPARC gene, is positioned in the 5q32-q33 region that shows a suggestive linkage with KTCN in previous studies (Bisceglia et al., 2009). A lot of single nucleotide variants in ZNF469 (Vincent et al., 2014;Lechner et al., 2014;Davidson et al., 2015;Karolak et al., 2016b;Kalantan et al., 2017;Lucas et al., 2017) and SPARC (De Bonis et al., 2011) have been identified in KTCN patients. However, the role of both genes in determining the disease has not been definitively clarified based on mutational studies and it should be further evaluated in the KTCN context. The proteins encoded by both ZNF469 and SPARC regulate ECM production and remodeling and participate in collagen homeostasis in the human cornea (Bradshaw et al., 2003;Burkitt Wright et al., 2011). Thus, the observed decreased expression of both genes in KTCN corneas might also partially explain previously identified downregulation of genes encoding ECM components and might shed more light on the role of these genes in KTCN.
The LOX gene encodes an enzyme responsible for collagen cross-linking in the cornea and is one of the well-studied genes in KTCN. It is mapped within 5q21 region linked to KTCN in Italian population (Bisceglia et al., 2009). While the first mutational screening of LOX has not confirmed the involvement of LOX variants in KTCN, evidence of a genetic association between KTCN and common SNVs located in the LOX gene has been revealed in Caucasian (Bykhovskaya et al., 2012), Czech (Dudakova et al., 2015, Iranian (Hasanian-Langroudi et al., 2015), Brazilian (Gadelha et al., 2020), and Han Chinese (Xu et al., 2020) KTCN patients. In our study, the level of LOX was decreased in KTCN corneas. Previously, the reduced expression or activity of LOX have been found in the cultured keratoconic fibroblasts (Dudakova et al., 2012), KTCN corneal epithelium and tears (Shetty et al., 2015), the cone apex of KTCN patients (Pahuja et al., 2016), as well as in the whole KTCN corneas (Kabza et al., 2017). The observed decreased level of LOX might affect the reduction of collagen cross-linking in the corneal stroma, leading to the thinning of the cornea (Dudakova et al., 2012;Pahuja et al., 2016).
Among the genes with detected increased expression in KTCN, the most interesting candidates are TGFBI, discussed above, and SLC4A11. Sodium bicarbonate transporter-like protein 11 (SLC4A11), functioning as an electrogenic Na + -coupled borate co-transporter, is involved in the stimulation of cell growth and proliferation via borate-dependent mitogen-activated protein kinase activation (Romero, Fulton & Boron, 2004;Park et al., 2004). SLC4A11 has been reported as involved in KTCN and its observed deregulation supports previous findings (Guan et al., 2012;Nowak et al., 2013). However, the mechanism through which they influence KTCN development is unknown and further research should be performed to interpret the obtained data.
In our study, hierarchical cluster analyses of RT-qPCR gene expression data correctly grouped the majority of KTCN and non-KTCN corneas. However, one non-KTCN sample (KR49) was misclassified with the KTCN cluster. As was hypothesized previously, the observed misclassification might be the consequence of the high heterogeneity of the control group, which could be recognized as a study limitation (Kabza et al., 2017). As it was not possible to obtain "healthy" corneas from living donors, corneal tissues have been ascertained from control individuals with different ocular phenotypes, including corneal ulcer, bullous keratopathy, and history of keratitis and ocular trauma. All of these diseases can result in corneal scarring. The misclassified KR49 cornea was obtained from a 38-years-old man referred for corneal transplantation after ocular trauma and had no other eye diseases. The tissue homeostasis during corneal scarring could be altered due to activation or inhibition of several molecular pathways. Depending on the stage of wound healing that results in ocular scarring, the gene expression of particular genes could be abnormal. Thus, we hypothesize that this could lead to the misclassification of our patient in hierarchical cluster analyses based on gene expression profiles. In addition, KR49 patient was the youngest individual among control participants, which could also be a potential confounding factor.
Another aspect of our study was to validate the RNA-Seq experiment performed previously in the same (except for KR49) human KTCN and non-KTCN corneas. The comparison of RNA-Seq data with results obtained from current RT-qPCR analysis showed a strong correlation between the two sets of data, indicating that the selected 36 genes were not false positives and suggesting high concordance between both methods. However, protein-based assays would provide further validity to observed findings.

CONCLUSIONS
In conclusion, in the present study, we have assessed the expression of 36 candidate genes across KTCN and non-KTCN human corneas to further investigate their contribution in KTCN etiology. We confirmed differential expression of 35 genes, including core elements of recently postulated KTCN pathways: the ECM formation, Hippo, and Wnt signaling. We found that the deregulation of genes identified using RNA-Seq was not incidental, and selected candidate genes were not false positives. Our results also showed abnormal expression of several other genes, such as LOX, SPARC, ZNF469, or TGFBI, in which single nucleotide variants have been frequently identified in KTCN individuals. Although the obtained results partially explain the molecular basis of KTCN etiology, our understanding of this complex disease is still rudimentary. To address this gap, further studies in a larger group of patients should be performed.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This study was supported by the National Science Centre in Poland, Grants 2012/05/E/ NZ5/02127 and 2018/31/B/NZ5/03280 (Marzena Gajecka). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author Contributions
Justyna A. Karolak conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Barbara Ginter-Matuszewska performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Katarzyna Tomela performed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Michal Kabza conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Dorota M. Nowak-Malczewska conceived and designed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Malgorzata Rydzanicz performed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Piotr Polakowski analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Jacek P. Szaflik analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Marzena Gajecka conceived and designed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft.

Human Ethics
The following information was supplied relating to ethical approvals (i.e., approving body and any reference numbers): The study protocol was approved by the Institutional Review Board at Poznan University of Medical Sciences (453/14 and 755/19).

DNA Deposition
The following information was supplied regarding the deposition of DNA sequences: The RNA-seq data are available in the Gene Expression Omnibus database (GSE77938).

Data Availability
The following information was supplied regarding data availability: The raw data are available in Table S4.

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.9793#supplemental-information.