Thyroid eye disease (TED) is one of the most common autoimmune inflammatory diseases [29]. In about 3–5% of patients, TED is accompanied by loss of vision and compressive optic neuropathy [4]. Exploring disease-related genes on the molecular level can contribute to the development of drugs for TED treatment. In this study, we used the differential expression analysis, WGCNA and lmQCM approaches to identify TED-associated modules and constructed a solid coexpression network for TED. Based on the overlap of important genes retrieved from the different approaches, 11 genes (ATP6V1A, PTGES3, PSMD12, PSMA4, METAP2, DNAJA1, PSMA1, UBQLN1, CCT2, VBP1 and NAA50) were identified as hub genes in the PPI network of the overlap coexpression genes. We also found that NFRKB, ZNF711, ZNF407 and MORC2 were key candidate TFs possibly regulating genes in the TED-associated coexpression network. Similarly, the genes in the coexpression network were mostly regulated by predicted miRs such as hsa-miR-144, hsa-miR-3662, hsa-miR-12136 and hsa-miR-3646. Functional enrichment analysis indicated that genes in the coexpression network were enriched in the biological processes related to proteasomal protein catabolic process, proteasome-mediated ubiquitin-dependent protein catabolic process and protein polyubiquitination and the pathways of endocytosis, ubiquitin-mediated proteolysis and amyotrophic lateral sclerosis. Furthermore, drug prediction based on the coexpression network indicated that most of genes were targets for enzyme inhibitors, chlorodiphenyl, finasteride and nefazodone. The study is the first, to the best of our knowledge, to construct such a stringent coexpression network for TED and to explore its regulatory factors and potential targeting drugs. The present findings are important in the understanding of the pathogenesis of TED, its diagnostics and the development of drugs for this disease.

Previous studies have revealed the clinical features and some key genes associated with TED. For instance, Grusha Ia et al. [30] found that the expression of HBD-2 may cause corneal damage in patients. Previously, several researchers have also identified hub genes by using bioinformatics analysis. Lee et al. [8] identified that 328 DEGs were associated with active TED through RNA sequencing. Many of these were associated with inflammation, cytokine signaling, adipogenesis, IGF-1 signaling and glycosaminoglycan binding. Gopinath et al. [3] demonstrated that cardiac calsequestrin gene CASQ2 was markedly up-regulated (~ 2.2-fold) in patients with the orbital disease and/or healthy people by microarray analysis. Zhu and Yang [7] indicated that POMC, IL-2, GNG3, CXCR4, TLR4, CSF1R, LPAR3 and CXCL8 were hub genes among many DEGs. It was also found that 71 DEGs were associated with TED in the comparative Toxicogenomics Database. Rosenbaum et al. [6] extracted RNA from the orbit of 83 volunteers and found that TED had fewer inflammatory marker characteristics than other ophthalmologic diseases through PCA. In the present study, by combining multiple bioinformatics approaches, we were able to construct, for the first time, a robust regulatory coexpression network. The regulatory interaction genes identified here might shed light in the understanding of the pathogenesis of TED.

There are various tools for weighted gene coexpression network mining, for example, quasi-clique merger (QCM), lmQCM, Markov clustering (ML) and WGCNA. Except for WGCNA, lmQCM is suitable for identifying small yet densely connected gene coexpression networks. A detailed description of each approach is provided in “Introduction” and “Methods” sections. Here, we performed lmQCM and WGCNA network and found that 11 and 13 modules were identified by WGCNA and lmQCM, respectively. The correlation between TED and both lmQCM modules and WGCNA modules was similar, ranging from − 0.31 to 0.39. However, we obtained a small number of genes in TED-specific modules through lmQCM compared with WGCNA. Moreover, 1488 genes of TED-specific modules were found in both WGCNA modules and lmQCM modules, suggesting that the genes obtained by both methods were the most credibly related to TED and that both methods are highly efficient in module discovery. Pathway enrichment analysis results showed that the pathways regulated by the genes in the TED-specific modules identified by lmQCM and WGCNA showed overlap in functional terms, indicating that both WGCNA and lmQCM were able to extract important genes involved in TED.

In eukaryotes, the binding of transcription factors and histone modification can accurately regulate gene expression. TFs can control the transcription rate of genetic information from DNA to messenger RNA by binding to specific DNA sequences. A number of studies have reported that transcription factors can participate in the activation and shutdown of gene functions. In order to explore the transcription factors regulating the expression of genes involved in TED, we used ToppGene to predict the TFs that regulate the expression of TED-specific genes. We found a total of 221 TFs that may be related to the expression of TED-associated coexpression genes. Among them, NFRKB, ZNF711, ZNF407 and MORC2 were the most prevalent, but these TFs and many of other identified TFs have not been reported in TED before. Thus, we can conclude that we identified new TED-related TFs which might be important in the research on the pathogenesis of TED. More relevant evidence needs to be provided in future research.

Ubiquitination is believed to be involved in the regulation of almost all life activities, including cell cycle, proliferation, apoptosis, gene expression, transcription regulation, inflammatory immunity, etc. Several studies have found that ubiquitination has an important regulatory role in autoimmune diseases [31,32,33]. Ubiquitination controls the development, activation and differentiation of T cells and maintains an effective adaptive immune response to pathogens and immune tolerance to self-tissues. Autoimmune diseases are not only directly affected by ubiquitination, but may also be regulated by ubiquitination-related enzymes. Wang and his colleagues found that mutations in the ubiquitin-conjugating enzyme E2L3 (UBE2L3) gene may be related to the high risk of Hashimoto’s thyroiditis (HT) in Han Chinese, indicating that ubiquitination may be involved in thyroid-related autoimmune diseases [34]. TED is an autoimmune disease, and its occurrence and development may be closely related to ubiquitination. In addition, a previous study reported that small ubiquitin-like modifier 4 (SUMO4) is a common autoimmune disease susceptibility gene in the Japanese population [35]. We found that the genes in the regulatory coexpression network were enriched in the biological processes related to proteasomal protein catabolic process, proteasome-mediated ubiquitin-dependent protein catabolic process and protein polyubiquitination and the pathway of ubiquitin-mediated proteolysis, indicating that ubiquitination plays a significant role in the pathogenesis of TED. The ubiquitin–proteasome pathway is responsible for the degradation of most intracellular proteins in eukaryotic cells [36]. Proteasome dysfunction is related to autoimmune diseases [37], neurodegeneration [38], heart [39] and other diseases. We found that genes in the regulatory network were also enriched in the proteasome pathway and the processes of proteasomal protein catabolic process, and proteasome-mediated ubiquitin-dependent protein catabolic process, which can explain a previous finding that the immune proteasome LMP2 can be used as a new target for the treatment of autoimmune hypothyroidism [40]. A previous study also demonstrated that TNF-α can promote mTOR-dependent proteasome-mediated PDCD4 degradation in orbital fibroblasts and that rapamycin can enhance TNF-α induction and IL6 secretion by inhibiting PDCD4 degradation in orbital fibroblasts [41]. Thus, the coexpression network constructed in this study may regulate the pathogenesis of TED via the ubiquitin–proteasome pathway (UPP). Although we only discussed some of the pathways, other pathways also deserve a more in-depth study.

In order to identify the interaction between proteins associated with the pathogenesis of TED and the related hub genes, we performed PPI analysis. We found 11 hub genes (ATP6V1A, PTGES3, PSMD12, PSMA4, METAP2, DNAJA1, PSMA1, UBQLN1, CCT2, VBP1 and NAA50) in the PPI network. Several studies reported that the PSMD12, PSMD14, PSMD5, PSMD6 were associated with the pathway of the proteasome [42,43,44,45]. PSMA1, PSMA4 and PSMB7 are the three subunits of the 20S proteasome [46]. The regulatory effects of PSMA1 and PSMD12 in the liver of thyroid hormone were previously reported [47], which is consistent with the down-regulation of PSMA1 and PSMD12 in TED as observed in the present work. The 11 hub genes may regulate TED by regulating the pathway of UPP, which is consistent with the findings in the functional enrichment analysis. At present, there is still little about the role of these hub genes in the regulation of thyroid diseases. Although we have discussed the regulatory function of PSMA1, PSMD12 in the thyroid, there is no report indicating that the 11 hub genes are directly related to the occurrence and progression of TED, and thus, more studies are needed to validate their regulatory role in TED.

The miRs are important regulators of gene expression. However, there are limited data on the involvement of miRs in the pathogenesis of TED. Up to date, only miR-130a, miR-146a and miR-155 and miR-1287-5p have been reported in TED [48,49,50]. In the present work, we predicted 4581 miRs targeting genes in the coexpression network involved in TED. The most prevalent miRs were hsa-miR-144, hsa-miR-3662, hsa-miR-12136 and hsa-miR-3646. Thus, our study provided a new set of miRs with relevant potential functions in the pathogenesis of TED. Further studies based on these miRs might enlighten our understanding on the pathogenesis of TED and open new avenues for the development of therapeutics for this disease.

To date, there is no effective drug for the treatment of TED. Few studies have indicated that some drugs including teprotumumab, tocilizumab, rituximab and mycophenolate can improve the outcomes of TED [51]. Herein, we predicted a set of drugs that could target the predicted gene coexpression network. Enzyme inhibitors and chlorodiphenyl (54% chlorine), finasteride, nefazodone, sodium arsenate, beta-methylcholine, propiconazole, pentachlorophenol, glucose, thimerosal, thapsigargin, succimer were the drugs with the highest number of target genes. The predicted genes may be efficient in the treatment of TED, which needs further experimental screening and validation.

Our study presents some limitations. Although TED-specific genes were found through the combined use of DEGs, lmQCM and WGCNA, due to the small sample size we used, there may be some deviations. In addition, our results were only based on public data, and more in vivo and in vitro experiments are still needed to verify our findings. We found 11 hub genes regulating the pathway of ubiquitin and proteasome through PPI analysis, but the role of these genes in the occurrence and development of TED still needs in-depth validation. In addition, we tried to identify the association of genes in the TED coexpression network to other diseases in ToppGene, but the analysis output revealed no disease association (no output data), which implied that the coexpression network might be specific to TED; however, tremendous work is needed to confirm this observation. We hereby call on more colleagues to experimentally verify our results.