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)(4)(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)(10)(11). Recent epigenetic studies have further proven the essential role of epigenetics in the pathogeneses of autoimmune diseases (12)(13)(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)(18)(19). Co-signaling molecules include costimulatory 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 cellsmediated 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 cellsmediated 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 Tcell 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.

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 log 2transformed 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 RNAseq data are highly replicable and can provide a more accurate assessment of gene expression than array (31)(32)(33)(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,

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 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.

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).

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 (
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).

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.

DISCUSSION
Co-signaling molecules play essential roles in modulating T cellsmediated 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 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.
the disruption of the balance between co-stimulation and coinhibition 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)(40)(41)(42). For instance, PD1 pathway was down-regulated in several autoimmune diseases such as RA, multiple sclerosis and T1DM (43)(44)(45), while inducible T cell costimulator (ICOS) is highly expressed in patients with SLE, IBD, or RA (46)(47)(48). Moreover, in vivo studies also have proven 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.
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)(60)(61)(62)(63)(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)(68)(69)(70)(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)(73)(74)(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/coinhibitory molecules were involved in the pathogenesis of GD, such as CTLA4, CD40, ICOS, and ICOSL (79)(80)(81)(82)(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/coinhibitory 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 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.
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 CD160related 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 cosignaling 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 cosignaling 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.

FUNDING
The present work was supported by grants from the National Natural Science Foundation of China (Grant Nos. 81670722 and 81873636).