Leptin receptor is a key gene involved in the immunopathogenesis of thyroid‐associated ophthalmopathy

Abstract Thyroid‐associated ophthalmopathy (TAO), the most common and severe manifestation of Graves' disease (GD), is a disfiguring and potentially blinding autoimmune disease. The high relapse rate (up to 20%) and substantial side effects of glucocorticoid treatment further decrease the life quality of TAO patients. To develop novel therapies, we amid to explore the immunopathogenesis of TAO. To identify the key immune‐related genes (IRGs) in TAO, we integrated the IRG expression profiles in thyrocytes from a GD patient set (GD vs healthy control) and a TAO patient set (TAO vs GD). Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), protein‐protein interaction (PPI) and receiver operating characteristic (ROC) curve analyses identified the leptin receptor (LEPR) gene as the key IRG in TAO immunopathogenesis. Gene set enrichment analysis (GSEA) suggested enrichment of the antigen presentation pathway in TAO patients with higher LEPR. Increased LEPR expression was validated in TAO orbital tissues, and weighted gene co‐expression network analysis (WGCNA) showed that cell adhesion processes were positively correlated with LEPR. Our study revealed that LEPR is a key gene in TAO immunopathogenesis and plays different roles in thyrocytes and orbital tissues. Our findings provide new insights into diagnostic and therapeutic biomarkers for TAO.


| INTRODUC TI ON
Graves' disease (GD) is a common autoimmune thyroid disease caused by the breakthrough of tolerance to the thyroid-stimulating hormone receptor (TSHR). With an annual incidence of approximately 20 to 50 cases per 100 000 persons, 1 GD is associated with several complications, the most frequent and severe one of which is thyroid-associated ophthalmopathy (TAO). 2 Simultaneously or within 18 months of developing GD, almost half of GD patients develop TAO, whose symptoms include a dry and gritty ocular sensation, photophobia and excessive tearing. 3 Signs of TAO are swollen extraocular muscles, expansion of orbital fat and connective tissue and, in severe cases, corneal ulceration and decreased visual acuity. 4 Although glucocorticoids are a preferred first-line treatment for TAO patients, approximately 20%-30% of patients do not respond to glucocorticoids, and 10%-20% experience relapse after treatment withdrawal. 3,5 Considering the unsatisfactory therapeutic effects and decreased quality of life associated with TAO, the development of novel therapies targeting its pathogenic mechanisms has always been recognized as imperative. 6 The pathogenesis of TAO is still unclear but appears to involve both adaptive and innate immune responses. Regarding adaptive immune responses, single-nucleotide polymorphisms (SNPs) in the IL-21 and CTLA4 genes have been associated with TAO. 7,8 Both CD4+ and CD8+ T cells, as well as B cells, were found to be present in the majority of examined orbits from patients with TAO, and a 2018 study suggested that the level of infiltration correlates with disease activity. 9 Additionally, orbital fibroblasts, considered key target cells in TAO, exhibit up-regulation of CD40, which causes their activation by CD40L on T cells; moreover, TSHR-expressing T cells might further stimulate adipogenesis via cyclooxygenase up-regulation. 10,11 Regarding innate immune responses, macrophage infiltration in orbital fat from patients with TAO is higher than that in controls, 12 possibly promoting a profibrotic phenotype in orbital fibroblasts through increased hyaluronic acid production and cell contractility. 13 Based on these previous observations regarding the immune mechanism of TAO, advanced treatments, such as mycophenolate mofetil (an inhibitor of T and B cell proliferation), rituximab (an anti-CD20 monoclonal antibody) and tocilizumab (an anti-IL6 receptor antibody), have been proposed. 14 However, these treatments have certain side effects. Thus, exploring the underlying pathogenesis and further developing novel therapies for TAO from the perspective of immune responses is a great need.
In this study, we first identified that the leptin receptor (LEPR) gene is the key gene involved in the immunopathogenesis of TAO LEPR expression was validated not only in thyrocytes but also in orbital tissues from TAO patients. Bioinformatic analyses revealed that up-regulation of LEPR promotes antigen presentation in thyrocytes and mediates protein secretion, immune cell adhesion and adipogenesis pathway activity in orbital tissues. Our findings suggest that LEPR contributes to the immunopathogenesis of TAO and is a promising diagnostic and therapeutic target.

| Study design and data collection
Three microarray data sets (E-MEXP-2612, GSE9340 and GSE58331) of GD or TAO patients were obtained from the NCBI Gene Expression Omnibus (GEO) and ArrayExpress in accordance with our selection criteria until 1 April 2020 ( Figure 1). E-MEXP-2612 was used to analyse gene expression profiles of GD in thyrocytes (2 healthy controls [HCs] and 6 patients with GD), and GSE9340 was used to analyse thyrocytes from TAO patients (8 GD patients without TAO and 10 GD patients with TAO). The GSE58331 data set was used to validate the findings regarding thyrocytes in orbital tissues from TAO patients (28 patients with non-specific orbital inflammation (NSOI), 24 patients with TAO and 21 HCs). The E-MEXP-2612, GSE9340 and GSE58331 data sets were generated utilizing the GPL571, GPL6014 and GPL570 platforms, respectively. The normalization and quality control of data in E-MEXP-2612 were carried out with the 'affy' R package (version 1.64.0), 15 and the remaining data sets were processed with the 'limma' R package (version 3.42.2). 16 In addition, a list of immune-related genes (IRGs), which were identified to participate in immune activity, was obtained from the Immunology Database and Analysis Portal (ImmPort) database. 17

| Screening of differentially expressed genes and principal component analysis
Data analysis was conducted with the 'limma' R package to detect differentially expressed genes (DEGs) between the GD and HC data or the TAO and GD data after normalization. The following cut-off criteria were used: |log 2 (fold change; FC)| > 0.1 and P value <.05.
Next, the differentially expressed IRGs were extracted from the complete set of DEGs. Two features were extracted from the differentially expressed IRGs or complete gene set of each group via an unsupervised principal component analysis (PCA) method. Heatmaps and volcano plots were generated using the 'pheatmap' (version 1.0.12) and 'ggplot2' (version 3.3.0) packages.

| Gene Ontology (GO) functional enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and protein-protein interaction (PPI) network construction
To explore the underlying mechanism by which differentially expressed IRGs impact the progression of GD and TAO, we conducted functional enrichment analyses, including GO term and KEGG pathway analyses, with the 'clusterProfiler' R package (version 3.14.3) 18 and set a P value of <.05 as the threshold.
Biological processes (BP), cellular components (CC) and molecular functions (MF) were included in the GO term enrichment analysis. The R package 'GOplot' (version 1.0.2) 19 was adapted to graphically show the results of GO analyses. The PPI network was constructed based on data collected from the STRING online database (https://strin g-db.org/; with a confidence score of >0.4 set as the cut-off value) and visualized with Cytoscape software version 3.8.0 20 (with genes with a degree of ≥3 considered key genes). In this part, the shared IRGs in the above analyses served as common key genes.

| Identification of key genes and common pathways
The 'VennDiagram' R package (version 1.6.20) 21 was used to conduct Venn diagram analysis for differentially expressed IRGs and pathways, which were obtained from data sets via GO, KEGG and PPI analyses. The diagrams show the shared IRGs (key genes in TAO) together with processes and pathways.

| Evaluation of diagnostic efficacy
The Pearson correlation test was used to compare the expression of key genes with disease states. A receiver operating characteristic (ROC) curve was plotted, and the area under the curve (AUC) was calculated with the 'pROC' R package (version 1.16.2) 22 to evaluate the capability of key genes to distinguish TAO from GD or TAO from NSOI. ROC curves were considered significant with AUC value >0.8.

| Gene set enrichment analysis and gene set variation analysis
To further explore the potential function of the LEPR gene in TAO,   24 were selected as the reference gene set, and a P value of <.01 was considered the threshold.

| Molecular characteristics of the LEPR gene
To explore the regulatory mechanism of LEPR in GD and TAO, we focused on the genes in the leptin signalling pathway that may be targets for LEPR. The gene list was obtained from the Reactome Pathway Database (https://react ome.org). The expression profiles of these genes were visualized using the 'forestplot' R package (version 1.9). The genes that were differentially expressed and significantly correlated with LEPR expression were considered potential target genes for LEPR in the TAO set.

| Weighted gene co-expression network analysis
Gene co-expression network analysis was specifically performed on TAO patient data in the GSE58331 data set using the R package weighted gene co-expression network analysis (WGCNA; version 1.69). 25 The similarity matrix was constructed by Pearson correlation analysis. Next, to achieve a scale-free topology (scale-free

| Statistical analysis
The statistical significance of differences between two groups was analysed by a non-parametric test or t test based on the data distribution characteristics. All analyses were performed with R3.6.2 software. A P value of <.05 was considered statistically significant.

| Identification of differentially expressed IRGs in the GD set
We identified 298 up-regulated and 897 down-regulated genes in the GD set in the E-MEXP-2612 data set (thyrocytes; GD vs HC; Figure S1A). From this set of genes, we extracted 75 differentially expressed IRGs, 17 specifically, 27 up-regulated and 48 down-regulated genes ( Figure 2A). The heatmap showed that the differentially expressed IRGs could discriminate the GD and HC groups ( Figure 2B), and the PCA plot of the GD group did not overlap with the profile of the HC group ( Figure 2C). Gene functional enrichment analysis revealed that 'tube development', 'plasma membrane' and 'hormone binding' were the most enriched biological terms ( Figure 2D,E). The KEGG pathway cytokine-cytokine receptor interaction was the most highly enriched with differentially expressed IRGs ( Figure 2F). The PPI network of differentially expressed IRGs was constructed, and 50 genes were identified as potential key IRGs in GD ( Figure 2G).

| Identification of differentially expressed IRGs in the TAO set
The expression profile of IRGs in the TAO set in the GSE9340 data set (thyrocytes; TAO vs GD) was also analysed. A total of 935 DEGs were identified ( Figure S1B). Among these genes, 58 were IRGs; 30 were up-regulated and 28 were down-regulated ( Figure 3A). The heatmap and PCA plot showed that the differentially expressed IRGs can clearly distinguish the TAO and GD groups ( Figure 3B,C). GO functional enrichment analysis results for the 58 IRGs showed that the most enriched terms were 'regulation of mononuclear cell migration' in the BP category, 'clathrin-coated vesicle membrane' in the CC category and 'hydrolase activity' in the MF category ( Figure 3D,E, Figure S2). Similarly, the cytokine-cytokine receptor interaction pathway was the most enriched KEGG pathways ( Figure 3F), and PPI network analysis identified 11 potential key IRGs in TAO ( Figure 3G).

| Venn diagram identification of the TAOspecific key IRGs and pathways among GDrelated factors
A total of 5 potential TAO-specific key IRGs were identified by Venn diagram analysis; these genes were detected based on the differentially expressed IRGs in the two sets ( Figure 4A). Next, we utilized  Figure 4C).
Additionally, 26 GO processes and 26 KEGG pathways (including cytokine-cytokine receptor interaction and Th1 and Th2 cell differentiation) were shared by these two sets ( Figure 4D,E, Appendix S1 and S2). Venn diagram analysis assisted us in detecting potential key IRGs and pathways in TAO pathogenesis.

| Linear regression and ROC curve analysis of the two shared key IRGs
Subsequently, we conducted linear regression and ROC curve analyses of the two potential key IRGs (THRB and LEPR). In the TAO set, both THRB and LEPR genes were positively correlated with the TAO state (THRB: P = .0207, LEPR: P = .0011; Figure 4F). ROC curve analysis showed that the AUCs were 0.7750 and 0.900 for THRB and LEPR, respectively ( Figure 4G). Considering that the AUC for LEPR was >0.8 but that for THRB was <0.8, we selected LEPR as the key gene involved in the immunopathogenesis of TAO as well as a potential early-warning biomarker and interventional target to prevent TAO in GD patients.

| GSEA: Robust immune response-related processes and pathways are prominent in the group of TAO patients with high LEPR expression
To determine the possible functional pathways affected by the LEPR gene in the TAO state, we performed GSEA to map the BP ( Figure 5).
In the group of TAO patients with high LEPR expression, antigen presentation-associated BP, including MHC protein complex and MHC class II protein complex, were up-regulated ( Figure 5A). Similarly, the KEGG pathways antigen processing and presentation pathway, complement and coagulation cascade pathway and type 1 diabetes mellitus pathway were obviously enriched in TAO patients with LEPR gene up-regulation ( Figure 5B). These results suggested that high LEPR expression on thyrocytes may contribute to the pathogenesis of TAO by up-regulating the antigen presentation pathway.

| Leptin signalling pathway-related genes in the TAO set
To investigate the molecular mechanism of the LEPR gene in the pathogenesis of TAO, we focused on the genes in the leptin signalling  Figure 5D). Thus, STAT5B was the potential target molecule of LEPR involved in the TAO pathogenesis in thyrocytes.

| LEPR expression in orbital tissues from TAO patients
The GSE58331 data set included orbital tissues from TAO patients, NSOI patients and HCs. Because the PCA plots indicated that the discrimination between samples from TAO patients and those from HCs was not as good as that between samples from patients with TAO and samples from NSOI patients ( Figure 6A), we considered the NSOI group as a control group. Consistent with previous results in thyrocytes from TAO patients, LEPR was obviously increased in orbital tissues from TAO patients compared with those from NSOI patients (P = 6.049e-08; Figure 6B). In addition, the LEPR gene was positively correlated with the TAO state (P = 8.933e-08; Figure 6B), and the ROC curve showed that LEPR expression can distinguish patients with TAO from controls (AUC: 0.8958; Figure 6C). In addition, the levels of PTPN11 and IRS2 in the leptin signalling pathway were significantly increased (P =.03 and 6.591e-09; Figure 6D), and the expression of both of these genes was positively correlated with LEPR expression (P = 6.168e-06 and 6.09e-11; Figure 6E).  Figure S4). Module-trait correlation analyses showed that multiple modules were positively related to increased expression of LEPR ( Figure 7C). Among these modules, the brown and blue modules were the most significantly associated with high LEPR expression in TAO patients (brown: P < 1e-04, blue: P = .004; Figure 7D; Appendix S3 and S4). GO analysis showed that the genes in brown and blue module were the most highly enriched in the MF term 'cell adhesion molecule binding' and the BP term 'neutrophil activation', consistent with the GSVA results for BIOCARTA pathways ( Figure 7E). In addition, considering that LEPR was related to lipid metabolism, we explored the expression profile of genes with the top 50% mad values in the 'adipogenesis' and 'fatty acid metabolism' hallmark gene sets. The results suggested that the expression of lipid metabolism-related genes can discriminate TAO patients with high LEPR expression from those with low LEPR expression and that genes such as ACADM, COL15A1 and ADH1C are up-regulated in TAO patients with high LEPR expression. These findings implied that LEPR participates in TAO progression via the lipid metabolism pathway ( Figure 7F).

| D ISCUSS I ON
In this study, we screened two public microarray data sets containing thyroid samples (E-MEXP-2612 and GSE9340) to investigate key IRGs in the processes of TAO and used the GSE58331 data set, containing orbital tissues from TAO patients, to validate our findings. We firstly found that LEPR is not only a key gene involved in the immunopathogenesis of TAO but also a diagnostic biomarker for TAO among GD patients. In addition, bioinformatic analyses were conducted to explore the potential functional processes and potential targets of LEPR.  (F) Heatmap showed expression profile of genes with top 50% mad value in 'adipogenesis' and 'fatty acid metabolism' from hallmark set. HC, healthy control; TAO, thyroid-associated ophthalmopathy; NSOI, non-specific orbital inflammation; PCA, principal component analysis Furthermore, we detected an increase in LEPR expression in orbital tissues from TAO patients, which corresponds to findings from previous studies that showed increased leptin expression in orbital tissues from patients with TAO and that adipocytes derived from orbital preadipocyte fibroblasts stimulated with TSHR can produce 6-37 times more leptin than those from controls. 33,34 The main pathological alteration in orbital tissues from TAO appears to involve both the extraocular muscles and the orbital fat compartments. 35 The results of hallmark and KEGG analyses via GSVA in TAO patients with high LEPR expression suggest that LEPR is involved in protein secretion. Enlarged extraocular muscles from patients with TAO are widely separated by an amorphous accumulation of granular material consisting primarily of collagen fibrils and glycosaminoglycans. 36 Thus, it is reasonable to infer that high LEPR expression in orbital tissues may promote the production of collagen fibrils, which can further induce muscle body oedema and exacerbate the TAO state. 37 The other pathological alteration in extraocular muscles is diffuse and focal immune infiltration. 35 Similarly, the results of GSVA of BIOCARTA pathways and WGCNA suggested that adhesion molecules participate in the progression of TAO with increased LEPR expression, which has been confirmed to be correlated with the progression of TAO and subsequently induces immune cell infiltration. [38][39][40] In summary, LEPR presumably participates in the enlargement of extraocular muscles during TAO progression by promoting collagen fibril production and immune cell infiltration.
Considering that orbital fat compartments were also involved in TAO and that LEPR was related to lipid metabolism, we explored the correlation between these factors in TAO and found that genes in the lipid metabolism pathway could discriminate TAO patients with high LEPR expression from those with low LEPR expression. Previous results have revealed increased mRNA expression of the adipocytespecific genes leptin, adiponectin, fatty acid synthase, adipocyte fatty acid binding protein (AP2) and PPARγ in TAO-affected adipose tissue, implying that de novo adipogenesis occurs in orbital tissues in TAO. 35 Additionally, non-targeted metabolite profiling of blood from TAO patients exhibited different profiles in the GD group and HC group, and the combination of proline and 1,5-anhydroglucitol was identified as a GO-specific modulator. 41 Moreover, plasma cholesterol is considered a risk factor for TAO, and a significant relationship was found between the therapeutic efficacy of intravenous methylprednisolone and pre-treatment triglyceride levels in TAO patients. 42,43 Taken together, these findings indicate that LEPR may contribute to the adipogenesis and lipid metabolism alterations observed in patients with TAO and may even result in the unsatisfactory treatment effects of glucocorticoids.
However, the current study has limitations. On the one hand, further research concerning the exact mechanism of the LEPR gene in TAO was not performed. On the other hand, the three data sets that we analysed utilized distinct platforms for gene expression analysis and were assembled from quite different populations. Therefore, LEPR expression still needs to be explored in TAO patients of more different races.
In summary, our study revealed that LEPR is a key gene in the immunopathogenesis of TAO and that its up-regulation promotes antigen presentation by thyrocytes and mediates protein secretion, immune cell adhesion and adipogenesis pathway activity in orbital tissues.
These findings provide new insights into diagnostic and therapeutic biomarkers for TAO, although the specific molecular mechanism and functional pathway of LEPR in TAO needs further exploration.

CO N FLI C T O F I NTE R E S T
The authors confirm that there are no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets used and/or analysed during the current study are available from the corresponding author on reasonable request. The public data source: GEO database (https://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress database (https://www.ebi.ac.uk/array expre ss/).