Atg16L1 as a Novel Biomarker and Autophagy Gene for Diabetic Retinopathy

Objective Accumulating evidence suggests the critical role of autophagy in the pathogenesis of diabetic retinopathy (DR). In the current study, we aim to identify autophagy genes involved in DR via microarray analyses. Methods Gene microarrays were performed to identify differentially expressed lncRNAs/mRNAs between normal and DR retinas. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses of lncRNA-coexpressed mRNAs were used to determine the related pathological pathways and biological modules. Real-time polymerase chain reactions (PCR) were conducted to validate the microarray analyses. Results A total of 2474 significantly dysregulated lncRNAs and 959 differentially expressed mRNAs were identified in the retina of DR. Based upon Signalnet analysis, Bcl2, Gabarapl2, Atg4c, and Atg16L1 participated the process of cell death in DR. Moreover, real-time PCR revealed significant upregulation of Atg16L1. Conclusion This study indicated the importance and potential role of Atg16L1, one of the autophagy genes, as a biomarker in DR development and progression.


Introduction
Diabetic retinopathy (DR) is a major contributor to vision loss in patients with diabetes mellitus [1]. DR incidence has been increasing rapidly in recent years, from 127 million in 2010 to a projected 191 million by 2030 [2]. The underlying key biochemical pathways may include genetic and epigenetic factors, polyol pathway activation, production of advanced glycation endproducts (AGEs), protein kinase C (PKC) activation, hexosamine pathway activation, and poly (ADP-ribose) polymerase upregulation. However, DR pathogenesis is complex and remains incompletely understood [3].
Autophagy is the primary intracellular catabolic mechanism mediating degradation and recycling of proteins and organelles. Due to its essential role in development, aging, starvation, cellular differentiation, and cell death, autophagy has attracted marked attention in recent years. Dysregulation of autophagy and lysosomal pathways is the hallmark of many diseases, from diabetes to neurodegenerative disorders and lysosomal storage diseases [4][5][6]. Moreover, growing data have suggested the crucial role of autophagy in DR pathophysiology [7,8]. Heretofore, there remains limited understanding regarding the exact autophagy genes involved in DR development and progression. The advent of microarray technology has facilitated detection of the comprehensive pattern of simultaneous transcript expression [9]. In this study, we aim to identify the autophagy genes involved in DR by microarray analyses. Experimental mice were randomized to receive HFD (60 kcal%) (Research Diets Inc. D12492i) or normal diet control (ND, D12450Bi). Eight months after ND or HFD, mice were anesthetized with 2% isoflurane. Blood glucose concentrations were measured 48 hours after STZ injection and weekly thereafter. Only animals with blood glucose levels exceeding 250 mg/dL were considered diabetic. After 10-12 weeks, the animals were sacrificed by pentobarbital overdose. The retinas were quickly removed, placed in liquid nitrogen, and stored at −80°C for biochemical measurement.

Microarray Analysis.
Total RNAs were isolated from the retinas of diabetic mice and age-matched controls using TRIzol reagent (Life Technologies, Carlsbad, CA, USA) and purified with an RNeasy mini kit (Qiagen, Valencia, CA, USA) per manufacturer's protocol. Microarray profiling was performed by mouse Clariom™ D Assay (Affymetrix Gene-Chip®, USA, an assay containing 65956 gene-level probe sets). Raw data were normalized at the transcript level by the TAC software (Transcriptome Analysis Console; version: 4.0.1) using Affymetrix default analysis settings (Robust Multichip Analysis workflow). The median summarization of transcript expressions was calculated.
Based upon differentially expressed lncRNAs, hierarchical clustering was performed by R package heatmap (version:  1.0.12). Gene Ontology (GO) analysis determined the main function of the differentially expressed genes, yielding the likely gene regulatory network on the basis of biological processes and molecular function. Specifically, a two-sided Fisher's exact test and chi-square test were used to classify the GO category. The false discovery rate (FDR) was calculated to correct the P value (the smaller the FDR, the small the error in judging the P value). Pathway analysis was conducted upon the differential genes identified, per Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Fisher's exact test was performed to select the significant pathway. The threshold of significance was considered P < 0:05. Based the interactions of genes in the KEGG database, global signal transduction network (Signalnet) was generated to demonstrate the interaction between the differentially expressed genes in treated groups. The visualization of network was built by software Cytoscape (version: 3.6.0).

Real-Time Quantitative PCR.
Quantitative real-time PCR (qRT-PCR) was applied to validate the selected genes from the microarray analyses. Total RNAs were extracted using TRIzol reagent (Invitrogen, Carlsbad, Canada) and reverse-transcribed per manufacturer's instructions. qRT-PCR reaction was monitored by the ABI Prism 7500 Sequence Detection System (Applied Biosystems, Foster City, CA) and run in duplicate for each sample. The PCR reaction mixture (20 μl) contained 2 μl of cDNA template, 0.6 μl forward and reverse primers, and 10 μl of 2×SYBR-Green PCR Mix (Takara). The level of mRNA expression was calculated from the fluorescence intensity (b-actin served as internal control). Primers targeting Atg4c (F: 5 ′ -GATG AAAGCAAGATGTTGCCTG-3′ and R: 5′-TCTTCCCTG TAGGTCAGCCAT-3′) and Atg16L1 (F: 5′-CAGAGCAGC TACTAAGCGACT-3′ and R: 5′-AAAAGGGGAGATTC GGACAGA-3 ′ ) were used for real-time RT-PCR amplification. The relative gene expression was calculated by the ΔΔ threshold cycle (Ct) method. Real-time PCR reaction was run in biological triplicates for each sample. Melting curve analysis was used to verify the product purity at the end of the PCR run.

3.1.
Overview of lncRNA-mRNA Microarray Analysis. To reveal a differential gene expression profile, hierarchical clustering analysis compared lncRNA-mRNA expression between diabetic and nondiabetic retinas (Figures 1(a) and 1(b)). Differentially expressed lncRNAs (with statistical significance) between the two groups were identified via volcano plot filtering (Figure 1(c)). The combined criteria of a P value <0.05 and fold change > 1:1 identified 2474 lncRNAs expressed differentially, including 1487 upregulated and 987 downregulated lncRNAs. 317 significantly increased and 642 decreased mRNAs were also identified. The top 10 differentially expressed lncRNAs and mRNAs between diabetic and nondiabetic retinas are listed in Tables 1 and 2, respectively.  The GO database revealed (of 21 upregulated GOs) photoreceptor cell morphogenesis, activation of MAPK activity, and autophagy were related to retinal function (Figure 2(a)). Of 11 downregulated GOs, only eye photoreceptor cell development was directly related to retinal function (Figure 2(b)). The summaries of those genes involved in the significant signaling pathways relevant to retinal function are listed in Table 3. KEGG analysis revealed highly enriched upregulated signaling pathways included phosphatidylinositol signaling system, glycosphingolipid biosynthesis, fatty acid degradation, PPAR signaling pathway, arginine and proline metabolism, and protein digestion and absorption (Figure 2(c)). Significantly downregulated pathways included biosynthesis of unsaturated fatty acids, protein export, nitrogen metabolism, and calcium signaling pathway (Figure 2(d)).

Construction of the lncRNA-mRNA Coexpression Network and Identification of Genes in the Process of Cell
Death. Based upon significant pathway and GO analysis, the Signalnet analysis screened the key genes associated with DR pathogenesis. A total of 369 key genes were identified in the transduction network. Specifically, Bcl2, Gabarapl2, Atg4c, and Atg16L1 were involved in the process of cell death. Figure 3 showed the Signalnet of cell death-related lncRNA-mRNA coexpression network. To evaluate the expression of these 4 specific genes in diabetic and nondiabetic retinas, hierarchical clustering analysis was further performed. The heatmap (Figure 4(a)) suggested the good classification of mRNA expression profile between those two groups.

Confirmation of Gene Expression in Autophagy by RT-PCR.
To further validate the microarray analysis results, we conducted qPCR assays for the genes involved in autophagy. qRT-PCR analysis revealed significantly upregulated expression of Atg16L1 compared to control, while Atg4c expression was unremarkable (Figure 4(b)). The potential role of Atg16L1 in DR pathogenesis was further supported.

Discussion
The pathogenesis of DR is hugely complex and likely implicates the dysregulation of many biochemical and molecular   [10,11], but the underlying responsible genes are unclear. To better understand the significance of autophagy in DR, we evaluated the associated genes by microarray analyses. Signalnet analysis suggests the involvement of Bcl2, Gabarapl2, Atg4c, and Atg16L1 in the process of cell death in DR. As expression of Atg16L1 was significantly increased, this particular molecule may be an important participant in DR development and progression. lncRNAs/mRNAs have garnered attention in recent years for their potential regulatory role in DR [12,13]. Additionally, the important roles of lncRNAs were investigated in PDR by harvesting fibrovascular membranes [14]. In the current study, we identified 2474 significantly dysregulated lncRNAs and 959 differentially expressed mRNAs, further confirmed by PCR analysis. In another prior study, lncRNAs were investigated with respect to DR pathogenesis in a mouse streptozotocin-induced diabetic model, utilizing microarray analyses. Only 303 aberrantly expressed lncRNAs were identified in the retinas of early DR [13]. Different animal models and protocol may lead to this discrepancy. Future bioinformatic analysis of the lncRNAs/mRNAs will be helpful to confirm these results.
Accumulating evidence suggests that autophagy plays an essential role and may act as a double-edged sword in DR [15,16]. Activation of autophagy results in cell survival under mild stress in DR, while dysregulated autophagy can lead to massive cell death during severe stress conditions. In high glucose concentrations, autophagy activation in cultured ARPE-19 cells decreases proinflammatory cytokine production [17]. High glucose upregulates autophagy in retinal Müller cells during early DR pathogenesis phases [10]. In agreement with these previous studies, we demonstrate (based on GO analysis) that autophagy is upregulated during DR. While the contribution of autophagy to DR development requires further study, our findings may open interesting perspectives for novel therapies.
Interestingly, two conserved mRNAs, including Atg16L1, were identified to be involved with autophagy in the current study. As a key player in early autophagy initiation, Atg16L1 also regulates subsequent steps of this pathway [18]. Moreover, recent progress has demonstrated that Atg16L1 may be involved in diabetic pathophysiology by regulating autophagy [19,20]. In consistent fashion, we report significantly increased expression of Atg16L1 in diabetic mouse retinas. Based upon these findings, we hypothesize that Atg16L1 may be an important component of the DR pathological process. Further studies are warranted to better understand the functional role of this novel mRNA in DR. In conclusion, we identified dysregulated autophagy genes involved in DR by microarray analyses. Atg16L1 may be a potential biomarker for the diagnosis and prognosis of DR. It may serve as a potential therapeutic target blocking and slowing DR progression. Future work is required to confirm these findings and elucidate the specific underlying molecular signaling mechanisms.

Data Availability
The data used to support the findings of this study are available from the corresponding authors upon request.

Disclosure
Our results have been accepted to be presented at the Association for Research in Vision and Ophthalmology (ARVO) Annual Meeting in 2020.  The results showed that, compared to control, Atg16L1 was significantly upregulated while Atg4c expression was unchanged in diabetic retinas.