Microarray Data of Lacrimal Gland Implicates Dysregulated Protein Processing in Endoplasmic Reticulum in Graves’ Ophthalmopathy

Graves’ ophthalmopathy (GO) has become one of the most common orbital diseases. Although some evidences announced the potential mechanism of pathological changes in extraocular muscle and orbital adipose tissue, little is known about that in lacrimal enlargement of GO patients. Thus, gene expression profiles of lacrimal gland derived from GO patients and normal controls were investigated using the microarray datasets of GSE105149 and GSE58331. The raw data and annotation files of GSE105149 and GSE58331 were downloaded from Gene Expression Omnibus (GEO) database. Bioinformatics including differentially expressed genes (DEGs), Gene Ontology, Kyoto Encyclopedia of Gene and Genome (KEGG) pathway, protein-protein interaction (PPI) network construction, hub gene identification, and gene set variation analysis (GSVA) were successively performed. A total of 173 overlapping DEGs in GSE105149 and GSE58331 were screened out, including 20 up-regulated and 153 down-regulated genes. Gene Ontology, KEGG and GSVA analyses of these DEGs showed that the most significant mechanism was closely associated with endoplasmic reticulum (ER). Moreover, we identified 40 module genes and 13 hub genes which were also enriched in the ER-associated terms and pathways. Among the hub genes, five genes including HSP90AA1, HSP90B1, DNAJC10, HSPA5, and CANX may be involved in the dysfunction of protein processing in ER. Taken together, our observations revealed a dysregulated gene network which is essential for protein processing in ER in GO patients. These findings provided a potential mechanism in the progression of lacrimal enlargement in GO patients, as a new insight into GO pathogenesis.


INTRODUCTION
Graves' ophthalmopathy (GO), also known as thyroid eye disease (TED), is an orbital disease that is uniquely linked to Graves' disease (GD), generally present in 25%-50% of GD patients (1). The main clinical features of GO include proptosis, lacrimal gland enlargement, eyelid retraction, diplopia, and exposure keratopathy, even developing into irreversible vision loss in severe cases (2). It is generally agreed that GO is an autoimmune disease that results in orbital remodeling, enlargement and fibrosis (3). These pathological changes involves three distinct but related immune processes: inflammation, adipogenesis and glycosaminoglycan accumulation (4). Most investigators have hypothesized that orbital fibroblasts are the main targets of inflammatory cytokines released by infiltrated immune cells in the extraocular muscle and orbital adipose tissue (5). Moreover, activated orbital fibroblasts could secrete many cytokines, such as interleukins 1b, 6,8,16, TNF-a, RANTES and CD154 (6). These cytokines promote orbital trafficking of monocytes and macrophages, facilitate differentiation of orbital fibroblasts and stimulate accumulation of hyaluronic acid-rich stroma, leading to uncontrolled immune responses (4)(5)(6). Although several mechanisms have been proposed for the development of GO, the exact pathogenesis of GO has not yet been illustrated.
Recently, high-throughput technologies provided ample evidence on the potential biomarkers and molecular mechanisms of GO. For examples, microarray analysis revealed the potential role of CASQ2 to trigger autoimmunity events (7), and identified several adipogenesis-related genes and some genes involved in Wnt and IGF-1 signaling as being potentially implicated in pathogenesis (8)(9)(10)(11). RNA-Seq also identified some meaningful genes (such as PTX3, HOXB2, HOXB3, CCL2, and SERPINA1) as potential biomarkers of GO (12)(13)(14). In addition, microRNA and protein sequencing found some blood circulating biomarkers that have the potential to diagnose GD, predict GO disease status and optimize patient management, such as hsa-miR-27a-3p, hsa-miR-22-3p, zonulin, haptoglobin, and lumican (15). All these accumulating findings improve our knowledge of the mechanisms underlying pathological changes in extraocular muscle and orbital adipose tissue of GO. Nevertheless, our understanding of the molecular mechanisms of lacrimal gland enlargement in GO is still limited.
In the present study, microarray analysis was performed to systematically investigate the key genes and pathways in lacrimal gland by mining the microarray datasets (GSE105149 and GSE58331) of GO and normal lacrimal samples. GSE105149 and GSE58331 were obtained from the public database GEO. These results would help to uncover a potential mechanism of lacrimal enlargement and offer new information for the pathogenesis of GO.

Characteristics of the Microarray Data Sets
The Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih. gov/geo/) is a public functional genomics data repository. GEO has collected a large amount of omics data, such as highthroughput gene expression data, chips and microarrays. Microarray datasets of GSE105149 and GSE58331 were obtained by searching the keywords including Graves' ophthalmopathy, thyroid eye disease and thyroid-associated ophthalmopathy in GEO. The raw data and annotation files of GSE105149 and GSE58331 were downloaded from the GEO database. Rosenbaum JT et al. contributed to the generation of GSE105149 and GSE58331and provided many significant findings from different perspectives (16)(17)(18)(19)(20)(21). GSE105149 was the microarray data of various lacrimal gland tissues, including sarcoidosis, granulomatosis with polyangiitis, thyroid eye disease, nonspecific orbital inflammation and healthy controls (16). Normal lacrimal samples from healthy individuals were obtained at the time of cosmetic surgery or blepharoplasty (16). GSE58331 was the microarray data of anterior orbit and lacrimal gland tissues, which were obtained from subjects with infl ammatory diseases such as nonspecifi c orbital inflammation, sarcoidosis, granulomatosis with polyangiitis, IgG4-associated disease and thyroid eye disease as wells as from normal controls (17)(18)(19)(20)(21). Normal control tissues were obtained during surgeries such as blepharoplasty and enucleation on eyes with non-inflamed orbits (17)(18)(19)(20)(21). The diagnoses of all samples were based on the clinical and histopathological information submitted by orbital disease specialists and ocular pathologists (16)(17)(18)(19)(20)(21). On account of GO or normal lacrimal gland as the inclusion criteria, the following samples were chosen in the present study: four TED lacrimal samples (GSM2823348, GSM2823349, GSM2823350, GSM2823351) and seven normal lacrimal samples (GSM2823306, GSM2823307, GSM2823308, GSM2823309, GSM2823310, GSM2823311, GSM2823312) were derived from GSE105149, and no sample was repeated; eight TED lacrimal samples (GSM1407195, GSM1407196, GSM1407199, GSM1407200, GSM1407202, GSM1407203, GSM1407204, GSM1407206) and seven normal lacrimal samples (GSM1407233, GSM1407235, GSM1407236, GSM1407237, GSM1407238, GSM1407240, GSM1407243) were obtained from GSE58331, in which eight TED samples were from four patients. Both microarrays have the same platform of GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array.

Screening of the Differentially Expressed Genes
According to the raw data and annotation files, microarray probes were converted into corresponding gene symbols, and the expression value of one gene symbol corresponding to multiple probes was calculated by taking the average value. In addition, Z-score was chosen to normalize the converted data to extend the range of differential expression. The Limma R package (http://www.bioconductor.org/packages/release/bioc/ html/limma.html) was used to detect differentially expressed genes (DEGs) between TED and normal lacrimal group. The cut-off criterion of DEGs was P< 0.05 and |log (fold change) | ≥1. The false discovery rate (FDR) was not used to identify DEGs in the present study. Nonetheless, we calculated FDR and published it alongside the p-values to give a better idea of the robustness of DEGs.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Analyses
Using the Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/), DEGs were mapped to Gene Ontology analysis to detect their possible functional processes from three terms of biological process (BP), molecular function (MF) and cellular component (CC). Moreover, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs was carried out to investigate their potential biological pathways by the KEGG Orthology-based Annotation System 2.0 (KOBAS 2.0, http://kobas.cbi.pku.edu. cn). P < 0.05 was considered as statistically significant difference.

Construction of Protein-Protein Interaction Network
DEGs were entered into the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/) to construct a functional protein-protein interaction (PPI) network. An interaction with the composite score > 0.4 was considered statistically significant. Then, PPI network complex was visualized with Cytoscape software. The central node genes were screened out with the cut-off criteria of connectivity degree ≥10. Whereafter, the Molecular Complex Detection (MCODE) in Cytoscape software was utilized to identify densely connected regions from PPI network complex. The connected region with MCODE score >3 and node number >4 was considered to be the significant module.

Gene Set Variation Analysis
Gene set variation analysis (GSVA) is a nonparametric and unsupervised clustering approach to estimate the score of gene-set enriched pathways based on the expression profile of each sample (22). For GSVA analysis, we chose HALLMARK gene sets as the reference, t value > 2 and p value < 0.05 as the cut-off to identify significantly altered pathways. The GSVA R package (http://www.bioconductor.org/packages/release/bioc/ html/GSVA.html) was used to calculate a pathway enrichment score between TED and normal lacrimal group. GSVA analysis for GSE105149, GSE58331, the combined data of GSE105149 and GSE58331 were performed, respectively.

Screening of the Common Molecular Origins Associated with GO Lacrimal Enlargement
Using P< 0.05 and |log (fold change) | ≥1 as the cut-off criterion, data from GSE105149 and GSE58331 were separately analyzed by the limma R package to identify DEGs. Representative heat maps showed different expression of all DEGs in the TED lacrimal samples relative to normal lacrimal samples ( Figure  1A). 80 up-regulated and 440 down-regulated genes were present in the TED lacrimal samples of GSE105149, 868 up-regulated and 252 down-regulated genes were in the TED lacrimal samples of GSE58331. The overlap between the two microarray datasets contained 173 genes, including 20 up-regulated and 153 downregulated genes in the TED lacrimal glands ( Figure 1B, Table  S1). Collectively, the common DEGs in two different gene expression profiles provided evidence on the common molecular origins of lacrimal gland enlargement in GO patients.

Investigation of the Potential Functions and Pathways in GO Lacrimal Enlargement
To investigate the mechanism, the overlapping DEGs were further measured through Gene Ontology and KEGG pathway analyses to explore the possible functions and pathways in the pathogenesis of GO lacrimal enlargement. Based on DAVID software, we observed a total of 81 Gene Ontology assignments, including 31 BPs, 9 MFs, and 41 CCs. The top five significant terms of Gene Ontology classification were illustrated in Figures 2A-C. In the BP category, most DEGs were involved in RNA splicing, protein folding in ER, cell-cell adhesion, mRNA processing and ER to Golgi vesicle-mediated transport. In the MF category, most DEGs may play roles in RNA binding, unfolded protein binding, peptidyl-prolyl cis-trans isomerase activity, cis-trans isomerase activity and nucleotide binding. In the CC category, a large percentage of the DEGs were associated with intracellular organelle lumen, nuclear lumen, cytosol, Golgi apparatus and ER. Furthermore, the KEGG analysis showed the dysfunction of protein processing in ER, phagosome and antigen processing and presentation pathways in GO lacrimal samples ( Figure 2D). Taken together, these findings suggested that multiple functions and pathways were significantly associated with GO lacrimal enlargement.

Identification of the Most Significant Pathway and Hub Genes in GO Lacrimal Enlargement
To identify the most significant DEGs "hub genes" associated with GO lacrimal enlargement, we performed a PPI network construction of all the overlapping DEGs on the basis of the STRING database. Finally, 157 nodes and 282 edges were present in the PPI network complex ( Figure 3A), in which 16 central node genes were selected with the cut-off criteria of connectivity degree ≥ 10 ( Figure 3B). Subsequently, module analysis was performed by the plug-in MCODE in Cytoscape software in order to identify the significant modules from the PPI network complex. Then, 6 modules successfully emerged from the above PPI network complex, including 40 genes. According to the selection criteria of MCODE score > 3 and node number > 4, module a, module b, module c and module d were regarded as significant modules ( Figure 3C). Among them, module A with the highest MCODE score and the maximum node number was considered as the most important module. Therefore, the common genes present in 16 central nodes and module A were selected as hub genes, including HSP90AA1, DNAJC10, HNRNPM, HNRNPD, HNRNPC, CANX, HSPA5, HSP90B1, SRSF11, DDX5, RBM25, HNRNPA2B1, and HSPA9 ( Figure  3C, Table 1). Besides, function and pathway analyses of DEGs in the modules or hub genes were further carried out. The functional annotations of the genes involved in the modules were protein folding in ER, poly (A) RNA binding, focal adhesion ( Figures 4A-C). Pathway analysis showed that the module genes mainly had a close relationship with protein processing in ER, spliceosome and pathogenic Escherichia coli infection ( Figure 4D). Moreover, the 13 hub genes were significantly enriched in multiple functions such as protein folding in ER, poly (A) RNA binding and extracellular matrix ( Figures 5A-C). And KEGG revealed that these hub genes were enriched in several pathways, especially protein processing in ER ( Figure 5D).
According to Gene Ontology and KEGG pathway analyses for the entire DEGs set, 40 module genes and 13 hub genes, we observed that protein processing in ER was the most significant pathway, in which five hub genes including HSP90AA1, HSP90B1, DNAJC10, HSPA5 and CANX were involved.

Differences in Pathway Activities between GO and Normal Lacrimal Gland
KEGG analysis requires relatively arbitrary selection of DEGs as input, only a small number of pathways were enriched in GO. In order to correct this selection to some extent, we applied GSVA analysis to investigate the critical pathways of GO. Since GSVA transforms data from a gene to a gene set based on the sample matrix in the coordinate frame, the enrichment pathway of each sample can be evaluated. As a result, GSVA analysis of TED versus normal lacrimal samples significantly revealed four upregulated and five down-regulated signatures in the TED lacrimal samples of GSE105149 ( Figure 6A), five up-regulated and 19 down-regulated signatures in the TED lacrimal samples of GSE58331 ( Figure 6B). The common enriched signals in both GSE105149 and GSE58331 were KRAS_SIGNALING_DN, MYOGENESIS, MYC_TARGETS_V1, UNFOLDED_PROTEIN_ RESPONSE, MTORC1_SIGNALING, PROTEIN_SECRETION and ANDROGEN_RESPONSE. In addition, the GSVA analysis of the combined data of GSE105149 and GSE58331 further supported the above results ( Figure 6C). These results further supported the differences in pathway activities between GO and

DISCUSSION
Currently, our understanding of GO pathogenesis is mainly focused on the extraocular muscle and orbital adipose tissue. The molecular origins of lacrimal gland enlargement in GO has not yet been elucidated. In this work, two gene expression profiles from different cohort studies were analyzed to identify DEGs with common changes between lacrimal gland from patients with GO and that from normal groups. Further bioinformatics analysis was performed to investigate the significance of these DEGs, including Gene ontology analysis, KEGG pathway analysis, PPI network construction, hub gene identification and GSVA analysis.
We identified a total of 173 common DEGs which may be involved in the most significant function and pathway in ER, including BPs "Protein folding in ER" "ER to Golgi vesiclemediated transport", MFs "Peptidyl-prolyl cis-trans isomerase activity" "Cis-trans isomerase activity" and CCs "Golgi apparatus" "ER" and pathways "Protein processing in ER" "Phagosome". As we known, ER is like a special factory for eukaryotic cells, in which at least one-third of molecules are synthesized, folded and transported (23). It maintains proper physiological proteostasis that the correctly folded molecules are modified and transported into the Golgi complexes, while the misfolded or unfolded molecules are trapped in the ER together with the molecular chaperone complex. When a large number of misfolded and unfolded molecules are accumulated, ER stress occurs and stimulates the unfolded protein response (UPR). Significantly, UPR intersects with different networks of inflammation, stress signals, and oxidative stress, in which NF- kB and JNK-AP1 pathways are involved (24). Inflammation, as an important pathological change of GO, directly affects the severity and activity of GO. Clinically, glucocorticoids can effectively alleviate the condition of some active moderate to severe GO patients by controlling the inflammatory response (25). In addition, previous reports have confirmed the dysfunctions of NF-kB pathway and oxidative stress in the progression of GO (26)(27)(28)(29). In this study, GSVA analysis also showed the dysregulation of protein secretion, oxidative phosphorylation, UPR and inflammatory response in lacrimal samples of GO. Therefore, these above evidences supported the ER and related signaling networks as potential mechanisms of GO lacrimal enlargement. Protein processing in ER was identified as the most significant pathway associated with the pathogenesis of GO lacrimal enlargement, in which 5 hub genes including HSP90AA1, HSP90B1, DNAJC10, HSPA5, and CANX were involved. These proteins respectively encoding by HSP90B1, DNAJC10, HSPA5, and CANX are all located in the ER. Of them, HSP90AA1 has the highest connectivity in the PPI construction of DEGs, which    encodes the well-known cytoplasmic Hsp90a, a highly conserved molecular chaperone that facilitates folding, maturation, structural maintenance and appropriate modification of target molecules (known as clients) (30). These clients are enriched in vesicle-mediated transport, telomere maintenance, mitotic signal transduction, apoptosis, cell-cycle progression, innate immunity and targeted protein degradation (31). For example, the binding of Hsp90a to lipopolysaccharide (LPS) can mediate LPS-induced inflammation, including the secretion of TNF by monocytes (32). Interestingly, LPS-induced inflammation has been proved to affect the differentiation of orbital pre-adipocytes in GO (33). HSP90B1 encodes another Hsp90 isoform GRP94, an ERresident molecular chaperone, which plays a vital role in the processing and transportation of secreting proteins (34). GRP94 shows many biochemical characteristics with other HSP90 family members, but also exhibits different activities, such as calcium binding, which is required for conditions in the ER (35).
Studies have shown that GRP94 functions as a key molecular chaperone for various Toll-like receptors and integrins (36,37), and also regulates innate and adaptive immunity (38). The absence of GRP94 in the hematopoietic system of mice leads to thrombocytopenia, prolonged bleeding, and massive platelets, which are clinically difficult to distinguish from human Bernardsoulier Syndrome (39). In addition, DNACJ10 and HSPA5 respectively encode the ER-resident proteins ERdj5 and BiP, which are crucial parts of ER-associated degradation (ERAD) complex involved in recognizing and degrading misfolded or unfolded molecules (40). ERAD is an essential quality-control system for ER molecules, in which the correctly folded molecules are transported to the secretion pathway and finally misfolded or unfolded molecules are translocated into the cytosol for proteasome degradation (41). Ushioda et al. confirmed that ERdj5 acts as a reductase that breaks the disulfide bonds of misfolded or unfolded molecules and facilitates ERAD by the  physical and functional connections with BiP (42). Furthermore, studies suggested ERdj5 and BiP as the master regulator of calcium homeostasis and redox homeostasis in the ER (43). Similarly, Calnexin (CANX) is a key ER-resident chaperone to help folding and quality control, thus ensuring that only correctly folded and assembled molecules can further enter into the secretion pathway (44). All these findings further deepen our understanding of the important roles of HSP90AA1, HSP90B1, DNAJC10, HSPA5 and CANX in protein processing in ER. Meanwhile, the observation that the downregulated expression of these 5 hub genes in GO lacrimal samples, further indicated that ER dysfunction may contribute to the enlargement of GO lacrimal gland, and may extend to other pathological manifestations of GO. In summary, our study presented a comprehensive bioinformatics analysis to identify potential key genes and pathways that may help to uncover the molecular mechanisms of the lacrimal gland enlargement in GO. Our findings firstly revealed that down-regulated expression of HSP90AA1, HSP90B1, DNAJC10, HSPA5, and CANX might generate the dysfunction of protein processing in ER, providing a novel explanation for the etiology and molecular events of the lacrimal enlargement in GO. When further molecular experiments confirmed the findings of the identified candidate genes and pathways in GO, the underlying molecular mechanisms of lacrimal gland enlargement in GO would be elucidated. All of these would contribute to a novel insight into GO pathogenesis.

DATA AVAILABILITY STATEMENT
The data sets (GSE105149 and GSE58331) for this study can be found in Gene Expression Omnibus (http://www.ncbi.nlm.nih. gov/geo/).

AUTHOR CONTRIBUTIONS
WT and YS conceived the experiments. JY, ZM, and XJ analyzed the data. WT drafted the manuscript. All authors contributed to the article and approved the submitted version.