Aberrant Expressions of Co-stimulatory and Co-inhibitory Molecules in Autoimmune Diseases.

Co-signaling molecules include co-stimulatory and co-inhibitory molecules and play important roles in modulating immune responses. The roles of co-signaling molecules in autoimmune diseases have not been clearly defined. We assessed the expressions of co-stimulatory and co-inhibitory molecules in autoimmune diseases through a bioinformatics-based study. By using datasets of whole-genome transcriptome, the expressions of 54 co-stimulatory or co-inhibitory genes in common autoimmune diseases were analyzed using Robust rank aggregation (RRA) method. Nineteen array datasets and 6 RNA-seq datasets were included in the RRA discovery study and RRA validation study, respectively. Significant genes were further validated in several autoimmune diseases including Graves' disease (GD). RRA discovery study suggested that CD160 was the most significant gene aberrantly expressed in autoimmune diseases (Adjusted P = 5.9E-12), followed by CD58 (Adjusted P = 5.7E-06) and CD244 (Adjusted P = 9.5E-05). RRA validation study also identified CD160 as the most significant gene aberrantly expressed in autoimmune diseases (Adjusted P = 5.9E-09). We further found that the aberrant expression of CD160 was statistically significant in multiple autoimmune diseases including GD (P < 0.05), and CD160 had a moderate role in diagnosing those autoimmune diseases. Flow cytometry confirmed that CD160 was differentially expressed on the surface of CD8+ T cells between GD patients and healthy controls (P = 0.002), which proved the aberrant expression of CD160 in GD at the protein level. This study suggests that CD160 is the most significant co-signaling gene aberrantly expressed in autoimmune diseases. Treatment strategy targeting CD160-related pathway may be promising for the therapy of autoimmune diseases.

Introduction

Autoimmune diseases are complex diseases in which aberrant immunity cause immune attacks and serious damage to normal human tissues or organs (1, 2). Despite considerable progress in the studies of both the risk factors and treatment of autoimmune diseases, the mechanisms of most autoimmune diseases remain largely elusive (3–5). Current literature supports both genetic and environmental factors have important roles in the development of autoimmunity (6, 7). Genome-wide association studies have provided deeper insights into the genetic causes of autoimmune diseases (8). Additionally, some environmental factors are important risk factors of autoimmune diseases, such as vitamin D deficiency and infections (6, 9–11). Recent epigenetic studies have further proven the essential role of epigenetics in the pathogeneses of autoimmune diseases (12–14). Nevertheless, the molecular mechanisms underlying the development of autoimmunity are still not clearly defined, and additional studies are needed (15, 16).

It has been well accepted that abnormal interactions between immune cells have crucial roles in the development of autoimmunity (17–19). Co-signaling molecules include co-stimulatory and co-inhibitory molecules and play essential roles in modulating the interactions of immune cells (20, 21). The activation and differentiation of T cells are directed by both the antigen-specific signal and costimulation signal, and T-cell costimulation is a critical part during the induction of T cells-mediated immune response (22). During intracellular contacts, specific recognition between co-signaling molecules can trigger changes in the expressions of functions of downstream elements and thus regulate TCR signals (22). Co-stimulatory molecules can enhance TCR-mediated immune responses and are pivotal in the activation of T cells. Co-inhibitory molecules can inhibit TCR-mediated immune responses, and regulation of T cells-mediated immune responses can be achieved by the expressions of co-inhibitory molecules on B cells, antigen-presenting cells (APCs), or peripheral tissues (23). Therefore, co-stimulatory and co-inhibitory molecules are important controllers of T-cell responses, and a precise balance between these pathways is crucial in preventing the development of autoimmune diseases (24).

T cells have obligatory roles in the pathogeneses of most autoimmune diseases (25, 26). Because T-cell costimulation is a critical part during the induction of T cells-mediated immune response, the roles of co-signaling molecules in autoimmune diseases also have gained much attention (21). There is increasing evidence that some co-signaling molecules have a pivotal role in the pathogeneses of autoimmune diseases, which may contribute to promising therapeutic targets (23, 27). There are multiple co-stimulatory and co-inhibitory pathways, but no study systematically assesses the aberrant expressions of these genes in autoimmune diseases. We thus sought to assess the expressions of co-stimulatory and co-inhibitory molecules in autoimmune diseases through both a bioinformatics-based study and a validation study using clinical samples.

Methods

Datasets of Autoimmune Diseases

NIH Gene Expression Omnibus (GEO) database was systematically searched to identify datasets of whole-genome transcriptome of autoimmune diseases, such as systemic lupus erythematosus (SLE), rheumatoid arthritis (RA), inflammatory bowel disease (IBD), juvenile idiopathic arthritis (JIA), type 1 diabetes mellitus (T1DM), autoimmune thyroiditis, Graves' disease (GD), psoriatic arthritis (PsA), Sjögren's syndrome and ankylosing spondylitis. To be included in our study, datasets must meet the following eligibility criteria: (1) Array datasets must contain at least 20 cases and at least 20 controls, while RNA-seq datasets should have at least 5 cases and at least 5 controls; (2) Assessing transcriptome in whole blood or peripheral blood mononuclear cells (PBMCs); (3) Identifying at least 50 differentially expressed genes (DEGs); (4) Raw data or gene expression profiling were available in GEO; (5) Containing the expressions of at least 90% of total co-signaling molecules analyzed in this study. To exclude the possible impact of treatment on the expressions of those co-signaling molecules, we only used data of the first visit for patients with transcriptome data of multiple visits.

RRA Discovery Study

To integrate the outcomes from multiple datasets, robust rank aggregation (RRA) method was utilized, which is a well-designed tool to analyze data from multiple arrays (28, 29). We firstly performed a RRA discovery study by using data from all array datasets. Both the gene expression matrix and related annotation document for each array dataset were downloaded from GEO database, and microarray probes were then mapped to gene symbols using corresponding annotation document. If multiple probes mapped to the same symbol, the mean value was used. The expression values of 54 co-stimulatory or co-inhibitory genes were subsequently extracted (Supplementary Table 1). Samples of each dataset were categorized into two groups including patients with autoimmune diseases and healthy controls. We uniformly used the log-transformed expression data, and those datasets without logarithmic transformation were firstly log2-transformed before calculating DEGs. The data above were then normalized using the “limma” package for R and DEGs were finally determined.

Before RRA analysis, both the up-ranked and down-ranked gene lists in each dataset were generated by their fold changes between cases and controls. The ranked gene lists of all eligible datasets were then finally integrated using “Robust Rank Aggregation” package for R software. P-value in the RRA tool indicated the possibility of ranking high in the final gene list, and genes with adjusted P-value less than 0.05 were considered as significant genes in the RRA analysis. Subgroup analysis by the type of data (whole blood or PBMCs) and the types of platforms (Illumina or Affymetrix). We also further added a sensitivity analysis by normalizing gene expression values across datasets through the ComBat of SVA R package (30).

RRA Validation Study

Previous studies have uncovered that array has several disadvantages, such as cross-hybridization risk, and RNA-seq data are highly replicable and can provide a more accurate assessment of gene expression than array (31–34). Therefore, we further performed a RRA validation study by using data from RNA-seq datasets. Raw RNA-seq read count data were analyzed using DESeq2 to identify DEGs (35). Both the up-ranked and down-ranked gene lists of those 54 co-signaling molecules in each dataset were generated by their fold changes between cases and controls. RRA analysis was then performed to integrate the outcomes above. Genes with adjusted P-value less than 0.05 were considered as significant genes in the RRA analysis.

Validation Study in SLE, RA, IBD, and JIA

We compared the expressions of CD160 and CD58 between cases and controls in four RNA-seq datasets including SLE (GSE112087), RA (GSE117769), IBD (GSE112057IBD), and JIA (GSE112057JIA), respectively. The transcripts per million (TPM) values of all genes, which could provide a more accurate measure of the true abundance of mRNA (36, 37), were calculated from the raw read counts. The TPM values of CD160 and CD58 were then extracted and were compared between cases and controls.

Samples Collection

Blood samples from newly diagnosed GD patients and healthy controls were collected in our hospital. GD was diagnosed by thyroid enlargement and hyperthyroidism, accompanied by increased thyrotropin receptor-stimulating antibody (TRAb). There were 42 GD patients and 36 healthy controls in the study of quantitative real-time PCR (qRT-CPR), and 23 GD patients and 21 healthy controls in the study of flow cytometry. The characteristics of those GD patients and controls were shown in the Supplementary Tables 2, 3. PBMCs from 2 mL EDTA anticoagulated peripheral blood were isolated by density gradient centrifugation using lymphocyte separation medium, and 1 mL TRIzol reagent were then added to store total RNA. PBMCs isolated from 5 mL heparin sodium anticoagulated peripheral blood were used for the study of flow cytometry. The study was approved by the Institute Ethics Committee of Shanghai University of Medicine & Health Sciences Affiliated Zhoupu Hospital, and written informed consent was obtained from all participants in accordance with the Declaration of Helsinki.

RNA Extraction and qRT-PCR

RNA was extracted with an RNeasy Kit (Qiagen), and complementary DNA (cDNA) was then synthesized from 1 μg RNA. qRT-PCR was performed using SYBR Green PCR Master Mix with a total reaction volume of 15 μl in triplicate. The reactions were performed using Applied Biosystems 7500 Real-Time PCR System. PCR program was as follow: 30 s at 95°C for one cycle, then 40 two-step cycles of 5 s at 95°C and 34 s at 63°C. We performed qRT-PCR using three house-keeping genes as reference including B2M, ACTB (beta-actin) and GAPDH, which are abundantly expressed in the immune cells of peripheral blood (38). Primer sequences were as follows: ACTB fwd, 5′-CAT TGC CGA CAG GAT GCA G-3′; ACTB rev, 5′-CTC GTC ATA CTC CTG CTT GCT G-3′; B2M fwd, 5′-CAT CCA TCC GAC ATT GAA GTT-3′; B2M rev, 5′- ACG GCA GGC ATA CTC ATC TTT-3′; GAPDH fwd, 5′-GGA GCG AGA TCC CTC CAA AAT-3′; GAPDH rev, 5′-GGC TGT TGT CAT ACT TCT CAT GG-3′; CD160 fwd, 5′-GCC AGA AGC CAG AAG TCA GGT ATC CG-3′; CD160 rev, 5′-CCT GTG CCC TGT TGC ATT CTT C-3′; CD58 fwd, 5′-AGA GCA TTA CAA CAG CCA TCG-3′; CD58 rev, 5′-ATC TGT GTC TTG AAT GAC CGC-3′; TNFRSF14 (HVEM) fwd, 5′-GTG CAG TCC AGG TTA TCG TGT-3′; TNFRSF14 rev, 5′-CAC TTG CTT AGG CCA TTG AGG-3′; BTLA fwd, 5′-CAT CTT AGC AGG AGA TCC CTT TG-3′; BTLA rev, 5′-GAC CCA TTG TCA TTA GGA AGC A-3′; TNFSF14 (LIGHT) fwd, 5′-ATA CAA GAG CGA AGG TCT CAC G-3′; TNFSF14 rev, 5′-CTG AGT CTC CCA TAA CAG CGG-3′. Gene expression was evaluated by the comparative CT method and normalized to reference genes.

Flow Cytometry

Co-signaling molecules usually express high in certain types of immune cells, and their expression levels are intensively related to their immune functions. Because CD160 is mainly expressed in human CD8+ T cells, and is weakly expressed in CD4+ T cells, B cells and dendritic cells (Supplementary Figures 1, 2), we further studied its expression on the surface of both CD8+ T cells and CD4+ T cells in GD patients through flow cytometry. To analyze the percentages of CD8+CD160+ T cells and CD4+CD160+ T cells, PBMCs were stained with FITC-conjugated anti-CD4 mAb (Becton-Dickinson, 555346), PerCP-Cy5.5-conjugated anti-CD8 mAb (Becton-Dickinson, 565310) and PE-conjugated anti-CD160 mAb (Becton-Dickinson, 562118). PBMCs were incubated at 4°C for 30 min in the dark with the optimal dilution of each antibody according to the manufacturer's instructions. The samples were analyzed on BD FACS Melody™ instrument (Becton-Dickinson).

Statistical Analysis

Data were shown as mean with standard error (SE) or interquartile range (IQR). Based on the patterns of data distribution, the difference in the expressions of candidate genes or the percentages of immune cells between cases and controls was analyzed using either Mann-Whitney U-test or unpaired t-test. We further analyzed the diagnostic role of CD160 in autoimmune diseases through receiver operating characteristic (ROC) curve, and the area under the ROC curve (AUCs) were calculated. A two-sided P < 0.05 suggested statistically significant difference. Analyses were performed using STATA (version 12.0).

Results

Characteristics of Included Datasets

According to the eligibility criteria, a total of 19 array datasets were considered eligible in the RRA discovery study (Table 1). Moreover, 6 RNA-seq datasets were included into the RRA validation study (Table 2). Table 1 showed the main characteristics of those 19 array datasets, such as GEO accession IDs and samples information (Table 1). Table 2 showed the main characteristics of those 6 RNA-seq datasets (Table 2). Those 25 datasets included a total of 2,292 patients with autoimmune diseases and 690 controls. Among those 19 array datasets (Table 1), 7 datasets (GSE121239, GSE65391, GSE61635, GSE80060, GSE93272, GSE17590, and GSE9006) longitudinally profiled the transcriptome of some patients with different visits. To exclude the possible impact of treatment on the expressions of those co-signaling molecules, we only used data of the first visit for those patients. There were no longitudinal gene expression values in those 6 RNA-Seq datasets.

Table 1

Characteristics of 19 publicly available array datasets in the RRA discovery study.

GSE ID	Participants	Tissues	Type of analysis	Platform
GSE121239	65 SLE patients and 20 controls	PBMCs	Array	GPL13158 (Affymetrix)
GSE65391	118 SLE patients and 32 controls	Whole blood	Array	GPL10558 (Illumina)
GSE81622	30 SLE patients and 25 controls	PBMCs	Array	GPL10558 (Illumina)
GSE50772	61 SLE patients and 20 controls	PBMCs	Array	GPL570 (Affymetrix)
GSE61635	79 SLE patients and 30 controls	Whole blood	Array	GPL570 (Affymetrix)
GSE49454	157 SLE patients and 20 controls	Whole blood	Array	GPL10558 (Illumina)
GSE17755SLE	22 SLE cases and 45 controls	Whole blood	Array	GPL1291 (Hitachisoft)
GSE80060	104 JIA patients and 22 controls	Whole blood	Array	GPL570 (Affymetrix)
GSE66795	131 Sjögren's syndrome patients and 29 controls	Whole blood	Array	GPL10558 (Illumina)
GSE84844	30 Sjögren's syndrome patients and 30 controls	Whole blood	Array	GPL570 (Affymetrix)
GSE73754	52 Ankylosing spondylitis patients and 20 controls	Whole blood	Array	GPL10558 (Illumina)
GSE51092	190 Sjögren's syndrome patients and 32 controls	Whole blood	Array	GPL6884 (Illumina)
GSE17755RA	112 RA cases and 45 controls	Whole blood	Array	GPL1291 (Hitachisoft)
GSE45291	493 RA patients and 20 controls	Whole blood	Array	GPL13158 (Affymetrix)
GSE93272	65 RA patients and 35 controls	Whole blood	Array	GPL570 (Affymetrix)
GSE17755JIA	57 JIA cases and 53 controls	Whole blood	Array	GPL1291 (Hitachisoft)
GSE17590	22 JIA cases and 21 controls	Whole blood	Array	GPL6106 (Illumina)
GSE9006	43 T1DM cases and 24 controls	PBMCs	Array	GPL96 (Affymetrix)
GSE3365	85 IBD patients and 42 controls	PBMCs	Array	GPL96 (Affymetrix)

SLE, systemic lupus erythematosus; RA, rheumatoid arthritis; IBD, inflammatory bowel disease; JIA, juvenile idiopathic arthritis; T1DM, type 1 diabetes mellitus; PBMCs, peripheral blood mononuclear cells; RRA, Robust rank aggregation.

Table 2

Characteristics of 6 RNA-seq datasets in the RRA validation study.

GSE ID	Participants	Tissues	Type of analysis	Platform
GSE122459	20 SLE patients and 6 controls	PBMCs	RNA-seq	Illumina HiSeq 2500
GSE112087	31 SLE patients and 29 controls	Whole blood	RNA-seq	Illumina HiSeq 2500
GSE117769	50 RA patients and 50 controls	Whole blood	RNA-seq	Illumina HiSeq 2500
GSE112057IBD	75 IBD patients and 12 controls	Whole blood	RNA-seq	Illumina HiSeq 2000
GSE112057JIA	115 JIA patients and 12 controls	Whole blood	RNA-seq	Illumina HiSeq 2000
GSE79970	85 JIA patients and 16 controls	PBMCs	RNA-seq	Illumina HiSeq 2500

SLE, systemic lupus erythematosus; RA, rheumatoid arthritis; IBD, inflammatory bowel disease; JIA, juvenile idiopathic arthritis; PBMCs, peripheral blood mononuclear cells; RRA, Robust rank aggregation.

Outcomes in the RRA Discovery Study

Findings from the RRA discovery study suggested that CD160 was the most significant gene aberrantly expressed in autoimmune diseases (Adjusted P = 5.9E-12), followed by CD58 (Adjusted P = 5.7E-06), CD244 (Adjusted P = 9.5E-05), LGALS9 (Adjusted P = 0.001) and CD27 (Adjusted P = 0.005) (Figure 1; Table 3). As shown in Figure 1, CD160 was down-regulated in most included datasets (Figure 1). The RRA outcomes after ComBat normalization was similar with that before normalization, and CD160 was still the most aberrantly expressed co-signaling molecule in autoimmune diseases (adjusted P-value = 1.3E-11, Supplementary Figure 3, Supplementary Table 4).

Figure 1

Heatmap shows those significant genes in the RRA analysis of 19 array datasets. Data from 19 array datasets were integrated using RRA method. CD160 was the most significant down-regulated gene in autoimmune diseases (Adjusted P = 5.9E-12), while CD58 was the most significant up-regulated gene (Adjusted P = 5.7E-06). The numbers in the heatmap were for the logarithmic fold change in each dataset which was calculated by the limma package. Red indicated increased expression, and green indicated decreased expression.

Table 3

Significant aberrantly expressed genes in the RRA analysis of 19 array datasets.

Gene name	Expression change	logFC	P-value	adjPvalue
CD160	Down	−0.56	1.1E-13	5.9E-12
CD58	Up	0.36	1.1E-07	5.7E-06
CD244	Down	−0.23	1.8E-06	9.5E-05
LGALS9	Up	0.32	1.9E-05	0.001
CD27	Down	−0.28	9.1E-05	0.005

logFC, the log2 of fold change; adjPvalue, adjusted P value; RRA, Robust rank aggregation.

In the subgroup analyses, both RRA analysis of 14 array datasets using whole blood and RRA analysis of 5 array datasets using PBMCs validated CD160 as the most aberrantly expressed co-signaling molecule in autoimmune diseases (Supplementary Tables 5, 6, Supplementary Figures 4, 5), which suggested that the aberrant expression of CD160 in autoimmune diseases was not affected by the type of transcriptomes. Additionally, both RRA analysis of 9 Affymetrix datasets and RRA analysis of 7 Illumina datasets consistently identified CD160 as the most aberrantly expressed co-signaling molecule in autoimmune diseases (Supplementary Tables 7, 8, Supplementary Figures 6, 7), which suggested that the aberrant expression of CD160 in autoimmune diseases was also not significantly influenced by types of platforms (Illumina or Affymetrix).

Outcomes in the RRA Validation Study

In the RRA validation study, RRA analysis of those 6 RNA-seq datasets validated CD160 as the most significant gene aberrantly expressed in autoimmune diseases (Adjusted P = 5.9E-09) (Figure 2, Table 4). However, findings from the RRA validation study suggested that CD58 was not aberrantly expressed in autoimmune diseases (Adjusted P = 0.99).

Figure 2

Heatmap shows those significant genes in the RRA analysis of 6 RNA-seq datasets. Data from 6 RNA-seq datasets were integrated using RRA method. CD160 was the most significant down-regulated gene in autoimmune diseases (Adjusted P = 5.9E-09), while CD58 was not an obviously significant gene in the RRA validation study. The numbers in the heatmap were for the logarithmic fold change in each dataset which was calculated by the limma package. Red indicated increased expression, and green indicated decreased expression.

Table 4

Significant aberrantly expressed genes in the RRA analysis of 6 RNA-seq datasets.

Gene name	Expression change	logFC	P-value	adjPvalue
CD160	Down	−0.52	1.1E-06	5.9E-09
PDCD1LG2	Up	0.51	0.0006	0.03
CD274	Down	0.50	0.0008	0.04

logFC, the log2 of fold change; adjPvalue, adjusted P value; RRA, Robust rank aggregation.

Validation in SLE, RA, IBD, and JIA

By comparing the TPM values of CD160 and CD58 from four RNA-seq datasets including SLE (GSE112087), RA (GSE117769), IBD (GSE112057IBD), and JIA (GSE112057JIA), we found that the aberrant expression of CD160 was statistically significant in SLE (P = 0.0006), IBD (P < 0.0001) and JIA (P = 0.0005), and was marginally significant in RA (P = 0.06) (Figure 3). CD58 was aberrantly expressed between SLE patients and controls (P < 0.0001), but it not obvious in RA (P = 0.55), IBD (P = 0.11), and JIA (P = 0.09) (Figure 4).

Figure 3

Validation of the aberrant expression of CD160 in SLE, RA, IBD, and JIA patients. The transcripts per million (TPM) value of CD160 were calculated from the raw read counts of four RNA-seq datasets including SLE (GSE112087), RA (GSE117769), IBD (GSE112057IBD), and JIA (GSE112057JIA), respectively. Difference between groups was analyzed using unpaired t-test. The aberrant expression of CD160 was statistically significant in SLE (P = 0.0006), IBD (P < 0.0001), and JIA (P = 0.0005), and was marginally significant in RA (P = 0.06).

Figure 4

Validation of the expression of CD58 in SLE, RA, IBD, and JIA patients. The transcripts per million (TPM) value of CD58 were calculated from the raw read counts of four RNA-seq datasets including SLE (GSE112087), RA (GSE117769), IBD (GSE112057IBD), and JIA (GSE112057JIA), respectively. Difference between groups was analyzed using unpaired t-test. CD58 was aberrantly expressed between SLE patients and controls (P < 0.0001), but it not obvious in RA (P = 0.55), IBD (P = 0.11), and JIA (P = 0.09).

CD160 had a moderate role in diagnosing SLE, IBD and JIA, and the AUCs for CD160 to diagnose SLE, IBD, and JIA were 0.756 (95%CI 0.63–0.88), 0.823 (95%CI 0.69–0.95) and 0.778 (95%CI 0.63–0.92) (Figure 5), respectively. CD58 had a moderate role in diagnosing SLE with an AUC of 0.799 (95%CI 0.69–0.91), but it had poor value in diagnosing RA, JIA, and IBD (Figure 5).

Figure 5

ROC curve to assess the diagnostic roles of CD160 and CD58 in autoimmune diseases. The diagnostic roles of CD160 and CD58 in autoimmune diseases were assessed by receiver operating characteristic (ROC) curve. Transcripts per million (TPM) values of CD58 and CD160 in four RNA-seq datasets including SLE (GSE112087), RA (GSE117769), IBD (GSE112057IBD), and JIA (GSE112057JIA) were used. The area under the ROC curve (AUCs) were shown in the figures. (A) Shows the diagnostic roles of CD160 and CD58 in SLE, and the AUC for CD160 and CD58 in diagnosing SLE is 0.756 (95%CI 0.63–0.88) and 0.799 (95%CI 0.69–0.91), respectively. When the testing threshold for the TPM value of CD160 was defined as 4.35, the sensitivity and specificity for CD160 in diagnosing RA were 75.9 and 67.7%, respectively. (B) Shows the diagnostic roles of CD160 and CD58 in IBD, and the AUC for CD160 and CD58 in diagnosing IBD is 0.823 (95%CI 0.69–0.95) and 0.649 (95%CI 0.50–0.80), respectively. When the testing threshold for the TPM value of CD160 was defined as 8.67, the sensitivity and specificity for CD160 in diagnosing IBD were 83.3 and 58.7%, respectively. (C) Shows the diagnostic roles of CD160 and CD58 in JIA, and the AUC for CD160 and CD58 in diagnosing JIA is 0.778 (95%CI 0.63–0.92) and 0.662 (95%CI 0.53–0.80), respectively. When the testing threshold for the TPM value of CD160 was defined as 8.67, the sensitivity and specificity for CD160 in diagnosing JIA were 83.3 and 53.9%, respectively. (D) Shows the diagnostic roles of CD160 and CD58 in RA, and the AUC for CD160 and CD58 in diagnosing RA is 0.605 (95%CI 0.49–0.72) and 0.449 (95%CI 0.33–0.57), respectively.

Aberrant Expression of CD160 in GD

Outcomes of qRT-PCR suggested that CD160 was aberrantly expressed in the PBMCs between GD patients and healthy controls (P < 0.05), but CD58, BTLA, LIGHT, and HVEM were not (P > 0.05) (Figures 6, 7, Supplementary Figure 8). Moreover, CD160 had a moderate role in diagnosing GD, and the AUC was 0.725 (95%CI 0.61 to 0.84) (Figure 8). When the testing threshold was defined as 0.5 for CD160, the sensitivity and specificity for CD160 in diagnosing GD were 100.0 and 38.1%, respectively. When the testing threshold was defined as 1.0 for CD160, the sensitivity and specificity for CD160 in diagnosing GD were 50.0 and 71.4%, respectively. The outcomes from flow cytometry confirmed that CD160 was differentially expressed on the surface of CD8+ T cells between GD patients and healthy controls and the percentage of CD8+CD160+ T cells was obviously lower in GD patients than that of healthy controls (P = 0.002, Figures 9, 10), which proved the aberrant expression of CD160 in GD at the protein level.

Figure 6

qRT-PCR outcomes suggested that CD160 was obviously aberrantly expressed in GD patients. Expression level of CD160 in the PBMCs were determined using qRT-PCR. Gene expression levels were calculated using the comparative CT method and were normalized to three reference genes including B2M, ACTB, and GAPDH, respectively. There were a total of 42 GD patients and 36 healthy controls. Difference between groups was analyzed using Mann-Whitney U-test.

Figure 7

qRT-PCR outcomes suggested that CD58 was not aberrantly expressed in GD patients. Expression level of CD58 in the PBMCs were determined using qRT-PCR. Gene expression levels were calculated using the comparative CT method and were normalized to three reference genes including B2M, ACTB, and GAPDH, respectively. There were a total of 42 GD patients and 36 healthy controls. Difference between groups was analyzed using Mann-Whitney U-test.

Figure 8

ROC curve to assess the diagnostic roles of CD160, CD58, BTLA, HVEM, and LIGHT in GD. The diagnostic roles of CD160, CD58, BTLA, HVEM, and LIGHT in GD were assessed by receiver operating characteristic (ROC) curve. There were a total of 42 GD patients and 36 healthy controls. The area under the ROC curve (AUCs) were shown in the figure, and CD160 had the best diagnostic role among those 5 genes with an AUC of 0.725. When the testing threshold was defined as 0.5 for CD160, the sensitivity and specificity for CD160 in diagnosing GD were 100.0 and 38.1%, respectively. When the testing threshold was defined as 1.0 for CD160, the sensitivity and specificity for CD160 in diagnosing GD were 50.0 and 71.4%, respectively.

Figure 9

Plots of CD8+CD160+ T cells and CD4+CD160+ T cells from flow cytometry. (A) CD4+ and CD8+ T subsets were gated by flow cytometry in the PBMCs of healthy controls. (B) CD4+CD160+ T cells were gated by flow cytometry in the PBMCs of healthy controls. (C) CD8+CD160+ T cells were gated by flow cytometry in the PBMCs of healthy controls. (D) CD4+ and CD8+ T subsets were gated by flow cytometry in the PBMCs of GD patients. (E) CD4+CD160+ T cells were gated by flow cytometry in the PBMCs of GD patients. (F) CD8+CD160+ T cells were gated by flow cytometry in the PBMCs of GD patients.

Figure 10

Flow cytometry suggested decreased percentage of CD8+CD160+ T cells in GD patients than healthy controls. There were a total of 23 GD patients and 21 healthy controls. Difference between groups was analyzed using unpaired t-test.

Discussion

Co-signaling molecules play essential roles in modulating T cells-mediated immune responses, but their roles in autoimmune diseases has not yet clearly defined. We systematically assessed the expressions of 54 co-stimulatory or co-inhibitory genes in 2,292 patients with autoimmune diseases and 690 controls. The study suggests that CD160 is the most significant co-signaling gene aberrantly expressed in autoimmune diseases, and the aberrant expression of CD160 is an important characteristic of autoimmunity including GD, which indicates that the dysfunction of CD160-related pathway exerts a critical role in the development of autoimmunity.

Excessive co-stimulation and/or insufficient co-inhibitory signaling can result in the activation of autoimmune cells and thus lead to the onset of autoimmunity. Previous studies have suggested that some autoimmune diseases are caused by the disruption of the balance between co-stimulation and co-inhibition signaling pathways (21). Some studies showed that some co-stimulatory or co-inhibitory genes were aberrantly expressed in patients with autoimmune diseases, such as Programmed Death-Ligand 1 (PDL1), Programmed cell death-1 (PD-1), and CD40 (39–42). For instance, PD1 pathway was down-regulated in several autoimmune diseases such as RA, multiple sclerosis and T1DM (43–45), while inducible T cell costimulator (ICOS) is highly expressed in patients with SLE, IBD, or RA (46–48). Moreover, in vivo studies also have proven knock-out or blockade of some co-stimulatory or co-inhibitory genes can promote or delay the development of autoimmunity, such as Cytotoxic T-Lymphocye Antigen 4 (CTLA-4), CD40 ligand (CD40L), and CD28 (49–53). Finally, findings from gene-association studies also have suggested that some polymorphisms in co-stimulatory or co-inhibitory genes are intensively related to risk of autoimmune diseases, such as CLTA4, CD40, and CD28 (54–58). Therefore, dysfunction of co-stimulatory and co-inhibitory pathways exerts critical roles in the pathogeneses of autoimmunity.

The findings in the present study suggest the critical role of CD160 in the pathogenesis of autoimmunity. CD160 is a member of an emerging co-stimulatory/co-inhibitory pathway, namely HVEM/CD160/BTLA/LIGHT pathway (59–64). HVEM is known as tumor necrosis factor receptor superfamily member 14 (TNFRSF14) and has dual functional activity by either binding to co-inhibitory receptors (CD160 and BTLA), or binding to a co-stimulatory receptor LIGHT (65). The combination of LIGHT with HVEM exhibits co-stimulatory signal and can promote the immune response, while the binding of BTLA or CD160 to HVEM exhibits co-inhibitory signal and can attenuate TCR signals (66). Therefore, HVEM/ CD160/BTLA/LIGHT pathway is a bidirectional switch which can critically regulate the activation of T cells (66). Some studies have indicated HVEM/CD160/BTLA/LIGHT pathway may has an important role in the etiology of autoimmunity. For instance, several studies using animal models of autoimmune diseases proved the involvement of HVEM/CD160/BTLA/LIGHT pathway in the development of autoimmune diseases (67–71). Additionally, the aberrantly expression of HVEM, BTLA, and LIGHT has been studied in some autoimmune diseases, such as RA, T1DM, and SLE, and there was disease-specific difference in the outcomes (63, 72–75). Additionally, this study also suggest that CD160 has a moderate role in diagnosing SLE, IBD, JIA, and GD (Figures 5, 8), suggesting that CD160 may be helpful for the diagnosis of autoimmune diseases. Most autoimmune diseases are difficult to diagnose, and there is still lack of good biomarkers for many autoimmune diseases. The AUC of the ROC analysis for CD160 was generally less than 0.80, suggesting that it only had a moderate role in diagnosing autoimmune diseases and should not be used in a clinical setting as a sole diagnostic marker. The value of the combination of CD160 with previously established diagnostic methods in the diagnosis of common autoimmune diseases may be promising, and can be explored in future studies.

There are few studies focusing on the role of CD160 in autoimmune diseases. A study by Hosokawa et al. reported that CD160 was aberrantly expressed in CD8+ memory stem T cells (TSCMs) of acquired aplastic anemia (AA) patients (76). Hosokawa et al. also reported that the percentage of CD8+CD160+ TSCMs was not different between SLE patients and controls (P > 0.05) (76), but the percentage of CD8+CD160+ T cells was not analyzed. A gene-association study suggested that CD160 rs744877 was associated with RA, suggesting the involvement of CD160 in the pathogenesis of RA (77). Another study by Bouma et al. revealed that Thiopurine treatment reduced the expression of CD160 in whole blood of Crohn's disease patients, and thiopurine might produce its effect by selectively affecting effector cytotoxic CD160-positive cells (78). Collectively, the role of CD160 in the pathogenesis of autoimmunity is still largely rudimentary. Besides, few studies have investigated the expression of CD160 in those autoimmune diseases using qRT-PCR or flow cytometry, and additional studies are needed.

Some studies suggested that several co-stimulatory/co-inhibitory molecules were involved in the pathogenesis of GD, such as CTLA4, CD40, ICOS, and ICOSL (79–83). Our study firstly assessed the expression levels of HVEM, CD160, BTLA, and LIGHT in GD patients, and found that CD160 was aberrantly expressed in GD patients (Figure 6), which had not been reported in previous studies. The outcomes from flow cytometry confirmed that the percentage of CD8+CD160+ T cells was obviously lower in GD patients than that of healthy controls (Figure 10), which proved the aberrant expression of CD160 in GD at the protein level. Our findings also suggested that CD160 had a moderate role in diagnosing GD. However, the clinical significance of CD160 in GD is still not well defined, and further studies are needed.

A better understanding of the roles of co-stimulatory/co-inhibitory pathways in autoimmune diseases have important implications for the development of novel therapeutics (84, 85). Currently, some therapeutic interventions by targeting aberrantly expressed co-stimulatory/co-inhibitory genes have been studied as promising therapeutic methods for autoimmune diseases, such anti-CD40L drugs and CTLA-4Ig (86–89). The findings from our study identify CD160 as a novel therapeutic target for the treatment of autoimmune diseases, which has important implications for future studies on CD160-related pathway. Future studies are recommended to explore the feasibility of treating autoimmune diseases by targeting CD160-related pathway.

In our study, we found that CD160 was the most significant co-signaling gene aberrantly expressed in autoimmune diseases, while several widely-studied co-signaling genes, such as ICOS, PD1 and CLTA4, were not identified as significant genes in the RRA analysis. One possible explanation is the disease-specific roles of those co-signaling genes in autoimmune diseases. For instance, PD1 was differentially expressed on the CD4+ and CD8+ T cells between RA and psoriatic arthritis patients (90). Another possible explanation is the use of transcriptome data in whole blood or PBMC but not in specific types of immune cells. There is high possibility for the existence of aberrant expression of some co-signaling genes in certain types of immune cells but not in whole blood or PBMCs. Owing to the limited transcriptome datasets from specific types of immune cells in GEO, we were unable to analyze the differentially expressed genes in specific types of immune cells by RAA. Therefore, further studies are recommended to explore autoimmunity-related co-signaling genes in specific types of immune cells for the existence of enough transcriptome datasets in the future. Besides, the results for other molecules such as CD58 and CD244 were not consistent in our subgroup analyses by types of datasets, which may be resulted from the moderate difference between cases and controls or the difference in the compositions of immune cells between whole blood and PBMCs. The expression levels of these co-signaling molecules in autoimmune diseases need to be explored in more future studies.

This study suggests that CD160 is the most significant co-signaling gene aberrantly expressed in autoimmune diseases, and its dysfunction is an important characteristic of autoimmunity including GD. However, the molecular mechanism underlying the role of CD160 in autoimmunity is largely elusive, and additional studies are warranted to uncover it. Moreover, future studies are recommended to explore the feasibility of treating autoimmune diseases by targeting CD160-related pathway.

Author Contributions

WH, JZ, and SL designed the study, collected data, performed statistical analyses, and wrote the final version of the manuscript. BW and QL participated in data collection and performed statistical analyses. XJ, QY, and RS participated in data collection. All authors approved the final version of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.