Differentially expressed genes in orbital adipose/connective tissue of thyroid-associated orbitopathy

Subjects and tissue samples

All human studies were conducted according to the Declaration of Helsinki principles and were approved by the Ethics Committee of the Affiliated Wuxi People’s Hospital of Nanjing Medical University (identifier, KY23013). We collected human orbital adipose/connective tissues from 43 to 80-year-old patients with TAO undergoing routine resection of prolapsed orbital fat in the Department of Ophthalmology, the Affiliated Wuxi People’s Hospital of Nanjing Medical University, from July 2021 to August 2022. The demographics of the patients are presented in Table 1 and Fig. S1. All TAO patients included in this study were diagnosed according to Bartley’s criteria, and tissues of control individuals obtained in plastic surgery were collected as control samples. All patients provided written informed consent.

Bulk RNA sequencing analysis (RNA-Seq)

The total RNA in tissues were extracted. To ensure the quality of the samples for transcriptome sequencing, the concentration and integrity of RNA samples were checked using a Nanodrop ND-2000 spectrophotometer and an Agilent Bioanalyzer 2100/4200, respectively. The qualified RNA samples were used for mRNA preparation and cDNA library construction. After library construction, the qualified libraries were sequenced using the Illumina NovaSeq 6000 using PE150 mode. Following an extracting and filtering quality control, we obtained high-quality, cleaned reads, and a follow-up analysis was then conducted (Table S1). All experiments were repeated three times with biological replicates. The statistical power of this experimental design, calculated in RNASeqPower is 0.96, based on a sequencing depth of 6 GB, CV of 0.4. We have uploaded the RNA-seq into the NCBI, the NCBI accession number is PRJNA971380.

DEGs and differential alternative splicing (DAS) analysis

We used FeatureCount (version 2.0.2) (Liao, Smyth & Shi, 2014) to quantify transcripts at the gene level. Differential expression analyses were performed with edgeR (version 3.3.3) according to the criteria of —log2 (FC)— > 1 and P value < 0.05.

Alternative splicing (AS) is the process by which different splice sites in precursor messenger RNA are selected to generate multiple mRNA isoforms, so AS is an important mechanism in creating proteome diversity and regulating gene expression in different tissues and developmental stages. To identify the number of different splicing events in TAO patients and controls, the software rMATS (version 4.0.2) was used (Shen et al., 2014), a new statistical method for robust and flexible detection of differential AS from replicate RNA-Seq data. Five main alternative splicing events, A3SS, A5SS, MXE, RI, and SE, were analyzed. A significance threshold of P value < 0.01 was used to define differential alternative splicing events.

Functional enrichment analysis

GO enrichment analyses for both the upregulated and downregulated genes were carried out using the R package topGO (The Gene Ontology Consortium, 2021) and the results were visualized using the REVIGO tool (http://revigo.irb.hr) (Supek et al., 2011). KEGG Orthology Based Annotation System (KOBAS) v3.0 (Bu et al., 2021) was used to perform the functional enrichment analysis. GSEA was carried out using the R package ‘clusterProfiler’ (Yu et al., 2012). The results are indicated in the appropriate figure legend and text.

The protein–protein interaction (PPI) network and hub gene identification

Construction of a PPI network was conducted using STRING (https://string-db.org/). We uploaded DEGs to STRING and obtained high-resolution bitmaps. By calculating the degree of connectivity, the hub genes in the PPI network were identified via cytoHubba, which is a plugin in Cytoscape software (version v3.9.1) (Shannon et al., 2003).

RNA quantification

Total RNA was extracted using the RNAiso Plus (Takara, Kyoto, Japan), according to the manufacturer’s instructions. Final RNA pellets were resuspended in nuclease-free H2O and then determine the purity and concentration by measuring the optical density at 260 nm and 280 nm (NanoDrop 2000c; Thermo Fisher Scientific, Waltham, MA, USA). Reverse transcription of total isolated RNA was performed using the PrimeScript RT master mix kit (Takara, Kyoto, Japan). Gene expression was measured by qRT-PCR. The data were analysed using the 2^−ΔΔCT method and normalized to the endogenous control GAPDH mRNA (for humans), and the amount of target gene mRNA expression in each sample was expressed relative to that of the control. Primer sequences for qRT-PCR were designed using Primer Express Software (Thermo Fisher Scientific, Waltham, MA, USA; Table S2).

Histological and immunohistochemical analysis

Human orbital adipose/connective tissues were obtained during orbital decompression and fixed overnight in 4% PFA (w/v) at 4 °C. The adipose sample was dehydrated through graded ethanol, and paraffin embedded. Histological sections of 5 µm were taken along the vertical meridian. Specimens were stained with H&E staining and observed under an Olympus BX-51 light microscope (Olympus, Tokyo, Japan). Standard immunohistochemical analysis with citrate antigen retrieval was performed with the antibodies against CD45 (#70257S; Cell Signaling, Danvers, MA, USA), Fibronectin (FN, #15613-1-AP; Proteintech, Chicago, IL, USA), and intercellular adhesion molecule 1 (ICAM1, #ab282575; Abcam, Cambridge, UK) to localize expression. Standard immunofluorescence analysis was performed to indicate F4/80 (#ab6640; Abcam, Cambridge, UK) expression, followed by Goat anti-Rabbit IgG (H+L) Cross-Adsorbed Secondary Antibody, Alexa Fluor™ 488 (#A-11008; Thermo Fisher Scientific, Waltham, MA, USA), and Goat anti-rat IgG (H+L) Cross-Adsorbed Secondary Antibody, Alexa Fluor™ 488 (#A-11006; Thermo Fisher Scientific, Waltham, MA, USA).

Statistical analysis

The results are expressed as the mean ± SD. Significance was established between the two groups using Student’s t test (paired t test). Age was compared using the t-test, and gender was compared using chi-squared tests. The data were analysed using GraphPad Prism 5 statistical software (Prism v5.0; GraphPad Software, La Jolla, CA, USA). A P value < 0.05 was considered statistically significant.

DEGs in orbital adipose/connective tissue samples of TAO patients

Deep sequencing identified 183 DEGs with the conditions of —log2(FC)— > 1 and P value < 0.05 between the orbital adipose/connective tissues of TAO patients and control individuals. Among these, 114 genes were upregulated, and 69 genes were downregulated. The fragments per kilobase million (FPKM) value of mRNAs shows that there is no abnormal expression in the three samples in each group (Figs. S2A–S2C). Principal component analysis (PCA) showed a significant separation between the two sets of samples (Fig. S2D). In our volcano plot and heatmap analysis of TAO-enriched genes, we showed the top 40 most DEGs in TAO samples compared to the controls (Figs. 1A, 1B). To identify and analyze the corresponding changes in these underlying functional DEGs, the enrichment analyses were employed.

DAS gene analysis

Alternative splicing (AS) refers to the process of selectively removing or retaining exons/introns during the initial transcription of DNA into RNA and further processing into mature mRNA, resulting in multiple transcripts of a gene. To learn the potential AS of TAO patients, five main types of AS events were analyzed using rMATS, including exon skipping (SE), intron retention (RI), alternative 5′splice site (A5SS), alternative 3′splice site (A3SS), and mutually exclusive exons (MXE) (Fig. 2A). We selected the DAS genes with a threshold of P value < 0.01. The numbers of A3SS, A5SS, MXE, RI, and SE events were 65, 57, 22, 18, and 477, respectively. SE was the most prevalent AS event in TAO patients, whereas RI was the least prevalent (Figs. 2B, 2C). This data suggests that an abnormal splicing process leads to specific splicing isoforms, which may have a close relationship with the occurrence and development of TAO.

Enrichment analyses of DEGs

To explore the functions of DEGs, functional enrichment analysis was performed on DEGs by linking them with biological phenomena and their underlying mechanisms.

GO annotation analyses

GO analysis is a common useful method for large-scale functional enrichment research, which can significantly distribute DEGs into the biological process (BP), molecular function (MF), and the cellular component (CC). The most significant GO terms of upregulated and downregulated DEGs are shown in Figs. 3A–3C, and detailed information is listed in Table 2.

In the GO terms of TAO samples, inflammation response was the main BP category, including inflammatory response, regulation of inflammatory response, acute inflammatory response, regulation of acute inflammatory response, and myeloid leukocyte migration (Fig. 3A). This suggests that the pathogenesis of TAO is closely related to the aberrant activation of inflammatory responses, which play a key role in the activation of orbital adipogenesis. The MF category was abundant in glycosaminoglycan binding, G protein-coupled receptor binding, signaling receptor binding, and extracellular matrix structural constituent (Fig. 3B). In addition, CC mainly displayed extracellular region, extracellular space, and cell surface (Fig. 3C).

GSEA

GSEA is a promising, widely used software package that derives gene sets to determine different biological functions between two groups. By GSEA, we identified that cytokine–cytokine receptor interaction, cytokine–cytokine receptor interaction, NF-kappa B signaling pathway, rheumatoid arthritis, TNF signaling pathway, and viral protein interaction with cytokine and cytokine receptor were the top five enriched pathways (Fig. 3E). In summary, the biological processes from the enriched GO terms, KEGG pathways, and GSEA for the DEGs were mainly involved in the regulation of inflammatory response, glycosaminoglycan binding and hyaluronic acid binding.

Cross with gene expression omnibus (GEO) database

We downloaded the microarray data of GSE185952 from the GEO database (Yue et al., 2021). This dataset contains six samples, including three TAO patients who underwent orbital decompression for proptosis correction and three control groups obtained from patients who underwent plastic surgery. We screened out the DEGs on the cut-off criteria with —log2 (FC)— > 1 and P value < 0.01. Intersection analysis was performed on the DEGs of the two independent samples. We obtained six co-upregulated genes, cartilage oligomeric matrix protein(COMP), interleukin-8(CXCL8), fmet-leu-phe receptor (FPR1), chemokine (c-c motif) ligand 4-like 1 (CCL4), protein s100-a9 (S100A9), and nf-kappa-b inhibitor zeta isoform x2 (NFKBIZ), and 14 co-downregulated genes, alx homeobox protein 1 (ALX1), protein kinase c-binding protein nell2 isoform x3 (NELL2), cadherin-4 isoform x1 (CDH4), puratrophin-1 isoform x1 (PLEKHG4), cochlin (COCH), low-quality protein: bcl-2-modifying factor (BMF), probable carboxypeptidase x1 (CPXM1), prolyl endopeptidase fap isoform x1 (FAP), low quality protein: hemicentin-2 (HMCN2), melanophilin isoform x1(MLPH), collagen alpha-6 chain (COL6A6), epidermal growth factor-like protein 6 isoform x1 (EGFL6), shc-transforming protein 3 (SHC3), and periplakin (PPL) (Fig. 4A).

The protein–protein interaction (PPI) network and hub gene

To better understand the molecular mechanism of TAO, we visualized the importance of the relationship between proteins of DEGs using Cytoscape software (Table S3, Fig. 4B). Moreover, we identified the top 10 HUB genes according to node degree among these target genes via the Cytoscape plug-in cytoHubba (Fig. 4C). Collectively, these results suggest that the core proteins CXCL8, Toll-like receptor-2 (TLR2), CCL4 and angiotensinogen (AGT) in the PPI network may be involved in the regulation of TAO pathogenesis.

Validation of the expression of DEGs

We confirmed the DEGs by qRT-PCR in orbital adipose/connective tissues to confirm the results of RNA-seq. In orbital adipose tissues, alkaline tissue-nonspecific isozyme isoform x1 (ALPL), ceruloplasmin isoform x3 (CP), and AGT were significantly upregulated, and protein mab-21-like 1 (MAB21L1), phosphoinositide 3-kinase gamma-subunit (PIK3C2G), and clavesin-2 (CLVS2) were significantly downregulated compared with controls (Figs. 5A, 5B). In orbital muscle, glutathione peroxidase 3 (GPX3) and alpha-1-antichymotrypsin isoform x1 (SERPINA3) were upregulated, while MAB21L1 and PIK3C2G were downregulated, but the difference did not achieve statistical significance (Figs. 5C, 5D). Among these genes, only R-spondin 1 (RSPO1) was significantly downregulated in orbital muscle tissues. Thus, the validation results by qRT-PCR are consistent with the RNA sequencing results.

Histology and inflammation in the orbital adipose/connective tissues of TAO patients and control individuals

H&E staining showed the morphology of orbital adipose tissue, and consistently indicated an increased level of the inflammatory cell infiltration (black arrows) in TAO patients compared with the control individuals (Fig. 6A). Meanwhile, we immunohistochemically stained sections for the CD45, a protein expressed on all leukocytes, and found that CD45 expression (black arrows) also increased in the TAO patients compared with controls (Fig. 6B). To identify and quantitate macrophages within adipose tissue, we detected the expression of F4/80 antigen, a marker specific for mature macrophages in the orbital adipose tissues. Indeed, the TAO group had significantly increased amounts of F4/80-positive macrophages, compared with the control groups.

Besides, in the orbital muscle tissues, the inflammatory markers, CD45 and ICAM1 also increased surrounding the myofibrils in patients with TAO compared with the control groups (Figs. S3A and S3B). Moreover, there was a potent increase in the expression of fibrotic proteins, including α-SMA and FN (Figs. S3C and S3D), indicating the orbital fibrosis or myositis in individuals with TAO.

Collectively, our results showed that there were enhanced inflammatory responses in orbital adipose/connective tissue and increased levels of fibrosis in the extraocular muscles among TAO patients.