Increased expression of hub gene CXCL10 in peripheral blood mononuclear cells of patients with systemic lupus erythematosus.

Systemic lupus erythematosus (SLE) is a prototypical autoimmune disease, characterized by overactive inflammation and aberrant activation of lymphocytes. Chemokine C-X-C motif ligand 10 (CXCL10) has an important role in the initiation and deterioration of SLE. However, the expression levels of CXCL10 mRNA in T-helper (Th) cells and B lymphocytes from SLE patients have remained elusive. In the present study, a Bioinformatics analysis of differentially expressed gene (DEG) profiles obtained from RNA sequencing data for three matched samples was performed to explore the hub genes, mainly through Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes and protein-protein interaction analysis. Furthermore, the expression of CXCL10 in peripheral blood mononuclear cells (PBMCs), CD4+ Th cells and CD19+ B cells of 108 subjects, including 66 SLE patients and 42 healthy controls, was confirmed by reverse transcription-quantitative PCR. In addition, 4 single-nucleotide polymorphism (SNP) loci in the 3'-untranslated region of CXCL10 were assessed using the Snapshot SNP genotyping assay. A total of 152 clustered DEGs mainly accumulated in immune-associated GO terms and interferon-associated pathways were identified. The expression of CXCL10, one of the central genes in the interaction network cluster (the degree of interaction, MCODE score=28.414), was 6.27-fold higher in SLE patients compared with control patients. Furthermore, CXCL10 mRNA was confirmed to be elevated in PBMCs and CD19+ B cells of patients with SLE (P<0.001 for the two cell types). However, no significant difference in CD4+ T lymphocytes was present (P=0.881). In addition, no polymorphism was identified in four selected loci from the samples. Taken together, the present results demonstrated that CXCL10, one of the hub genes in the pathogenesis of SLE, is upregulated in PBMCs and B lymphocytes of patients with SLE, although none of the SNPs selected for analysis in the present study were identified to have any potential associations with SLE. Copyright: © Zhang et al.

Introduction

Systemic lupus erythematosus (SLE) is a chronic inflammatory autoimmune disease (1) that may be influenced by genetic variation (2), environmental stimuli or estrogen (3,4). Pathological hallmark features include the overproduction of auto-antibodies secreted by overactive B lymphocytes, accompanied by aberrant T-cell function, accumulation of immune complexes and excessive production of pro-inflammatory cytokines (2). It has been demonstrated that SLE may be activated by a variety of environmental stimuli, including strong ultraviolet (UV) light (5), ambient fine particulate and infections (6,7). Of note, high UV irradiation is a feature associated with high altitude conditions of the Yunnan plateau (YP), China (8). However, to date, only few studies (9,10) have focused on the pathogenesis of SLE in YP, and studies into the relationship between SLE and high altitude have not been previously reported. Therefore, the mechanism of the pathological characteristics associated with SLE in populations exposed to high UV remains to be fully elucidated.

Accumulated evidence indicates that the type-I interferon (IFN) system is involved in cutaneous autoimmune diseases. IFNs induce the expression of pro-inflammatory cytokines and chemokines, which support the cellular immune response (11). Elevated IFN levels in the serum and increased IFN expression in peripheral blood mononuclear cells (PBMCs) have been indicated to serve a key role in the pathogenesis of SLE (12). In the present study, dozens of differentially expressed genes (DEGs) in SLE patients vs. normal controls were identified by using high-throughput sequencing. Further computational biological analysis revealed central and interactional targets. Among these potential key genes, chemokine C-X-C motif ligand 10 (CXCL10), as a chemokine, also known as IFN-inducible protein (IP-10) has drawn our attention; CXCL10 has a significant role in inflammation and immunity by binding to chemokine C-X-C motif receptor 3 (13).

CXCL10 is highly expressed in a wide range of human diseases. It is involved in the pathological processes of three major human disorders, namely infectious, inflammatory and autoimmune diseases, including Graves' disease, psoriasis, type I diabetes and SLE (13–15). In SLE patients, the serum levels of CXCL10 are highly elevated and correlate with the level of disease activity (16,17). However, to the best of our knowledge, the mRNA expression of CXCL10 in SLE patients particularly in crucial lymphocytes, including T-helper (Th) cells and B lymphocytes, has remained elusive and its molecular mechanisms in SLE has not been previously explored.

In the present study, DEG profiles in PBMCs from SLE patients and healthy controls (Ctrl) were first analyzed. Subsequently, the expression levels of CXCL10 in PBMCs, CD4+ Th cells and CD19+ B lymphocytes were examined using reverse transcription-quantitative PCR (RT-qPCR). Furthermore, a tentative experiment was performed to screen for polymorphisms in the 3′-untranslated region (UTR) of CXCL10 to explore any possible association with susceptibility to SLE.

Materials and methods

Participants

A total of 3 patients meeting the American College of Rheumatology (ACR) classification criteria (18) for SLE [SLE disease activity index (SLEDAI) score ≥4] and three gender-matched healthy individuals were enrolled in the high-throughput screening experiment (all female subjects, with an average age of 39.08±7.23 and 35.71±6.83 years, respectively). A total of 66 patients meeting the ACR criteria (SLEDAI score ≥4) and 42 unrelated, unaffected controls were then recruited for mRNA detection. The characteristics and laboratory data of the cohort are listed in Table I. All cases of inpatients were recruited from the Department of Dermatology in the Second Affiliated Hospital of Kunming Medical University (Kunming, China), and the health controls were recruited from the check-up department of the same hospital at random during the same period (January to October 2018). Informed consent for the anonymized analysis of their blood samples, in addition to personal data regarding all the parameters listed in Table I were obtained from all subjects. The present study was approved by the Ethics Committee of the Second Affiliated Hospital of Kunming Medical University (Kunming, China).

Table I.

Demographic data of the two groups and SNPs in the 3′-untranslated region of chemokine C-X-C motif ligand 10.

Characteristics	SLE (n=66)	Ctrl (n=42)	P-value
Age (years; mean ± SD, M (Q3-Q1))	39±7.23 39 (45.25–34)	36±6.83 36 (42–28.75)	0.826
BMI (kg/m2; mean ± SD)	22.44±2.35	21.10±2.75	0.293
Gender [n (%)]			0.958
Female	63 (95.45)	40 (95.24)
Male	3 (4.55)	2 (4.76)
Ethnicity [n (%)]			0.793
Han	57 (86.36)	37 (88.10)
Others	9 (13.64)	5 (11.90)
Disease duration [months; M (Q3-Q1)]	78 (156.0–21.5)	NA	NA
SLEDAI [score; M(Q3-Q1)]	9 (15–7)	NA	NA
24-h urinary protein [g; M (Q3-Q1]	0.29 (1.39–0.06)	ND	NA
Low C3 [<0.79 g/l, n (%)]	54 (81.82)	ND	NA
Low C4 [<0.10 g/l, n (%)]	47 (71.21)	ND	NA
Anti-dsDNA [+, n (%)]	14 (21.21)	ND	NA
Anti-ENA [+, n (%)]	19 (28.79)	ND	NA
Anti-Sm [+, n (%)]	16 (24.24)	ND	NA
SSA (anti-Ro) [+, n (%)]	40 (60.71)	ND	NA
SSB (anti-La) [+, n (%)]	2 (3.03)	ND	NA
SNPs
rs35795399 (C/G/T)	0/0/66	0/0/42	NA
rs58658570 (A/T)	66/0	42/0	NA
rs34836828 (−/G)	66/0	42/0	NA
rs148141229 (A/G)	66/0	42/0	NA

BMI, Body mass index, where <18.5 was defined as underweight, 18.5–24.0 was defined as normal and >24.0 was defined as overweight; SD, standard deviation; M (Q3-Q1), median (quartile3-quartile1); NA, not available; ND, not detected; Ctrl, healthy controls; SLE, systemic lupus erythematosus patients; SNP, single nucleotide polymorphism; SLEDAI, SLE disease activity index; C3, complement 3; dsDNA, double-stranded DNA; ENA, extractable nuclear antigen; Sm, Smith antibody.

High-throughput sequencing and analysis

PBMCs were purified from heparinized venous blood by density-gradient centrifugation over Lymphoprep™ (Axis-Shield), performed by fractionation on Lymphoprep™ and separation by density-gradient centrifugation (400 × g, 20 min, 18–20°C). To isolate mononuclear cells, the Lymphoprep™ interface was aspirated without disturbing the erythrocyte/granulocyte pellet and then washed once with Phosphate-buffered saline (PBS). The cells were then counted and re-suspended in PBS at a concentration of 2×106/ml. Furthermore, the leukocytes deposited at the bottom were stored at −80°C.

RNA-sequencing (RNA-seq) was performed by Vazyme Biotech Co., Ltd. After mRNA was extracted, concentrated and sheared into fragments, complementary DNA (cDNA) was synthesized, followed by purification of fragments, terminal repair, polyA-tailing and ligation of Illumina sequencing adapters. Finally, raw data from a chain-specific library were obtained by amplification using qPCR. The Q-value was strictly >30 according to the quality-control assessment criteria. As a procedure for the final report, genes were annotated according to the National Center for Biotechnology Information (NCBI, http://www.ncbi.nlm.nih.gov/refseq/) (19), University of California Santa Cruz Genome Browser database (https://genome.ucsc.edu) (20) and Ensembl genome databases (http://www.ensembl.org/Homo_sapiens/Info/Index) (21) after data were cleaned according to alignment, assembly and qualification. mRNA abundance was estimated by the expected number of fragments per kilobase of transcript sequence per million base pairs sequenced. DEGs, identified by comparison of SLE patients and controls with a cutoff of |log2 fold change| ≥1 and P≤0.05, were included in the final report.

Cell isolation

Isolation of PBMCs was performed as previously described within 2 h after the blood specimen was obtained. CD4+ cells were enriched by positive selection with CD4 magnetic microbeads (Miltenyi Biotec) according to the manufacturer's protocol. To obtain highly purified CD4 cells, PBMCs were incubated with the microbeads (20 µl/107 cells) at 4°C for 15 min prior to being passed through LS columns using a quadroMACS separator (Miltenyi Biotec). During the process, degas buffer containing PBS (pH 7.2), bovine serum albumin (0.5% BSA, Sigma-Aldrich; Merck KGaA) and EDTA (2 mM) was used for washing the columns. A purity of >90% was confirmed by fluorescence staining with CD4-FITC (OKT-4; cat. no. 11-0048-42; eBioscience; Thermo Fisher Scientific, Inc.) and analysis using an LSR Fortessa fluorescence-activated cell sorting instrument (BD Biosciences). The same procedure was used to isolate and identify CD19+ lymphocytes by using B Cell Isolation Kit II (cat. no. 130-091-151; Miltenyi Biotec GmbH) from the remaining cells according to the manufacturer's instructions.

RT-qPCR

RNA was extracted from PBMCs, Th cells and B lymphocytes using the RNAclean Kit (Tiangen Biotech Co., Ltd.). cDNA was synthesized from total RNA with a cDNA synthesis kit (Takara Biotechnology Co., Ltd.). qPCR for cDNA sample amplification was performed using the Roche 480 real-time PCR system with SYBR Green Master Mix (KAPA Biosystems; Roche Diagnostics). The sequences of the primers are listed in Table II. The relative expression levels of genes were normalized to those of human GAPDH, calculated using the 2−∆∆Cq method and were log-transformed (22).

Table II.

Primer sequences.

Gene name	Primer sequence (5′-3′)
IFI27	F: CGTCCTCCATAGCAGCCAAGAT
	R: ACCCAATGGAGCCCAGGATGAA
OLFM4	F: GACCAAGCTGAAAGAGTGTGAGG
	R: CCTCTCCAGTTGAGCTGAACCA
IFI44	F: GTGAGGTCTGTTTTCCAAGGGC
	R: CGGCAGGTATTTGCCATCTTTCC
SLC12A1	F: AGGCTCTTTCCTACGTGAGTGC
	R: GCCACTGTTCTTGGTAAAGGCG
KRT1	F: CAGCATCATTGCTGAGGTCAAGG
	R: CATGTCTGCCAGCAGTGATCTG
TUBB2A	F: TTGGGAGGTCATCAGCGATGAG
	R: AGGCTCCAGATCCACCAGGATG

F, forward; R, reverse; IFI27, interferon α inducible protein 27; OLFM4, olfactomedin 4; IFI44, interferon-induced protein 44; SLC12A1, solute carrier family 12 member 1; KRT1, keratin 1; TUBB2A, tubulin β 2A class IIa.

Selection and detection of single-nucleotide polymorphisms (SNPs)

A total of four candidate SNPs were investigated in the 3′UTR of CXCL10: rs35795399, rs58658570, rs34836828 and rs148141229. These had been identified from Genome-Wide Association Studies from the NCBI (https://www.ncbi.nlm.nih.gov/snp/?term=CXCL10). Genomic DNA was extracted from the remaining leucocytes using the TIANamp Blood DNA Kit (Tiangen Biotech Co., Ltd.) according to the manufacturer's protocol. The gene polymorphisms were analyzed using the SNaPshot assay (Sangon Biotech Co., Ltd.). To confirm the genotyping results, PCR-amplified DNA samples were examined by DNA sequencing and the results were 100% concordant.

Bioinformatics analysis

Networks of DEGs were algorithmically generated based on the potential connectivity of their products using the Database for Annotation Visualization and Integrated Discovery online database (DAVID, version 6.8; http://david.ncifcrf.gov) (23), the Search Tool for the Retrieval of Interacting Genes and proteins (STRING) database (version 10.5; http://string-db.org) (24) and Cytoscape (version 3.6.1) (25), which is an open-source Bioinformatics software platform for visualizing molecular interaction networks. Enrichment analysis of DEGs, including Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, was performed using DAVID. STRING was used to analyze the functional interactions between proteins and construct protein-protein interaction (PPI) networks. In the present study, the threshold of protein interaction was set to medium confidence, which is the default setting of the website.

Furthermore, the plug-in MCODE (version 1.4.2) of Cytoscape that is able to visualize PPI networks was used to cluster a given network based on topology to identify densely connected regions. Finally, hub genes were selected with degrees of connectivity >10.

Statistical analysis. Statistically significant differences between two groups were determined using an unpaired Student's t-tests for data conforming to a normal distribution, or a Mann-Whitney U-test for those with an abnormal distribution. Kolmogorov-Smirnov test was used to determine whether the data conformed to normal distribution. Pearson's χ2 test was used to investigate any association with gender and ethnicity. A two-tailed P<0.05 was considered to indicate statistical significance. All statistical tests were performed using GraphPad Prism 7.0 (GraphPad, Inc.).

Results

Bioinformatics analysis of DEGs

To identify the most significant genes, DEGs were first sorted by their ascending P-value and then listed by the absolute values of logarithmic conversion (|log2SLE/Ctrl|) in descending order. A positive value indicated higher expression in SLE patients relative to controls; correspondingly, a negative value indicated a lower relative expression. The top 20 significantly up- or downregulated DEGs are listed in Table III. Of note, the expression of IFN-α-inducible protein 27 (IFI27) in the SLE group was 102.62 times that of the value in the control group. Of all of the DEGs, 76.97% (117/152) were upregulated and 23.03% were downregulated.

Table III.

Details of the top 20 up-/downregulated differentially expressed genes.

	Upregulation			Downregulation
Rank	Gene name	log2(SLE/Ctrl)	P-value	Gene name	log2(SLE/Ctrl)	P-value
1	IFI27	6.68123	5.00×10−5	SLC12A1	−4.01723	5.00×10−5
2	FBRSL1	4.38689	5.00×10−5	INTS1	−3.73539	5.00×10−5
3	OLFM4	3.91972	5.00×10−5	PRR12	−3.20664	5.00×10−5
4	SPTSSB	3.90728	5.00×10−5	TUBB2A	−2.46642	5.00×10−5
5	MMP8	3.65456	5.00×10−5	FAM153B	−2.38627	5.00×10−5
6	IFI44L	3.57125	5.00×10−5	ETS1	−2.145136	5.00×10−5
7	SIGLEC1	3.258	5.00×10−5	SLC38A7	−2.14331	5.00×10−5
8	RSAD2	3.24366	5.00×10−5	HECW2	−2.09099	5.00×10−5
9	USP18	3.2127	5.00×10−5	PFDN4	−1.8775	5.00×10−5
10	DLGAP1	3.15073	5.00×10−5	KRT1	−1.73932	5.00×10−5
11	EPSTI1	3.03301	5.00×10−5	DDIT4	−1.56707	5.00×10−5
12	MPO	3.00225	5.00×10−5	FOS	−1.4891	5.00×10−5
13	G0S2	2.99775	5.00×10−5	RNF182	−1.34662	5.00×10−5
14	IFI44	2.78607	5.00×10−5	MPZL1	−1.32971	5.00×10−5
15	LTF	2.77532	5.00×10−5	FCRL2	−1.26173	5.00×10−5
16	AGBL2	2.77292	5.00×10−5	FCER1A	−1.2243	5.00×10−5
17	LY6E	2.70831	5.00×10−5	BCAP29	−1.21236	5.00×10−5
18	ISG15	2.67862	5.00×10−5	KLRF1	−1.07264	5.00×10−5
19	CXCL10	2.64825	5.00×10−5	GNLY	−1.06844	5.00×10−5
20	SPATS2L	2.62536	5.00×10−5	ADGRG1	−1.06523	5.00×10−5

Ctrl, healthy controls; SLE, systemic lupus erythematosus patients; FBRSL1, fibrosin like 1; SPTSSB, serine palmitoyltransferase small subunit B; MMP8, matrix metallopeptidase 8; IFI44L, interferon unduced protein 44 like; SIGLEC1, sialic acid binding Ig like lectin 1; RSAD2, Radical S-Adenosyl Methionine Domain Containing 2; USP18, ubiquitin specific peptidase 18; DLGAP1, DLG associated protein 1; EPSTI1, epithelial stromal interaction 1; MPO, myeloperoxidase; G0S2, G0/G1 switch 2; LTF, lactotransferrin; AGBL2, ATP/GTP binding protein like 2; LY6E, lymphocyte antigen 6 family member E; ISG15, ISG15 ubiquitin like modifier; CXCL10, chemokine C-X-C motif ligand 10; SPATS2L, spermatogenesis associated serine rich 2 like; INTS1, integrator complex subunit 1; PRR12, proline rich 12; TUBB2A, tubulin β 2A class Iia; FAM153B, family with sequence similarity 153 member B; ETS1, ETS proto-oncogene 1, transcription factor; SLC38A7, solute carrier family 38 member 7; HECW2, HECT, C2 and WW domain containing E3 ubiquitin protein ligase 2; PFDN4, prefoldin subunit 4; DDIT4, DNA damage inducible transcript 4; FOS, Fos proto-oncogene, AP-1 transcription factor subunit; RNF182, ring finger protein 182; MPZL1, myelin protein zero like 1; FCRL2, Fc receptor like 2; FCER1A, Fc fragment of IgE receptor Ia; BCAP29, B cell receptor associated protein 29; KLRF1, killer cell lectin like receptor F1; GNLY, granulysin; ADGRG1, adhesion G protein-coupled receptor G1.

To cluster DEGs into ontology subsets, they were further analyzed using DAVID. GO analysis indicated that in the category Biological Process (BP), the DEGs were significantly enriched in the following terms: ‘Defense response to virus’ (GO:0051607), ‘type I interferon signaling pathway’ (GO:0060337), ‘response to virus’ (GO:0009615), ‘negative regulation of viral genome replication’ (GO:0045071), ‘innate immune response’ (GO:0045087), ‘complement activation’ (GO:0006956, GO:0006958), ‘phagocytosis’ (GO:0006911, GO:0006910) and ‘positive regulation of B-cell activation’ (GO:0050871) (Fig. 1A). In the category Cellular Component (CC; Fig. 1B), the significantly enriched genes were concentrated in the following terms: ‘Blood microparticle’ (GO:0072562), ‘circulating immunoglobulin complex’ (GO:0042571), ‘external side of plasma membrane’ (GO:0009897) and ‘nucleosome’ (GO:0000786). Furthermore, in the category Molecular Function (MF), the DEGs were mainly enriched in ‘immunoglobulin receptor binding’ (GO:0034987), ‘antigen binding’ (GO:0003823), ‘serine-type endopeptidase activity’ (GO:0004252), ‘double-stranded RNA binding’ (GO:0003725) and ‘2′-5′-oligoadenylate synthetase activity’ (GO:0001730) (Fig. 1C). All of the aforementioned terms were significantly enriched by the DEGs (P<0.05). The most significantly different and the most abundant DEGs identified from the sequencing results were clustered in the term ‘innate immune response’, followed by enrichment of virus- and type 1 IFN-associated BP terms, and immune-associated GO terms in the categories MF, CC and BP.

Figure 1.

Enrichment analysis of DEGs in high-throughput RNA sequencing. Mapping of Gene Ontology annotation and enrichment terms in the categories (A) Biological Process, (B) Cellular Component and (C) Molecular Function. (D) KEGG pathway enrichment analysis of significant DEGs. The upper x-axis/red graph displays the negative logarithm of the P-value, whereas the lower x-axis and blue bars represent the number of genes involved with each term. (E) Protein-protein interaction network generated by the Search Tool for the Retrieval of Interacting Genes and proteins and (F) interaction network of 30 hub genes obtained with the plug-in MCODE of Cytoscape. The color represents upregulation or downregulation, as defined in the chart. DEG, differentially expressed gene; KEGG, Kyoto Encyclopedia of Genes and Genomes.

KEGG pathway analysis revealed that the DEGs were significantly enriched in ‘systemic lupus erythematosus’ (hsa05322; 12 DEGs), followed by ‘influenza A’ (hsa05164) and ‘herpes simplex infection’ (hsa05168), as presented in Fig. 1D.

STRING analysis was employed to determine the roles of various DEGs in the pathogenesis of SLE. The PPI network was obtained and is displayed in Fig. 1E. The network statistics were as follows: Number of edges, 678; number of nodes, 139; average node degree, 9.76; average local clustering coefficient, 0.486; PPI enrichment P-value <1.0×10−16. This level of enrichment indicated that the proteins of DEGs at least partially biologically interact as a group during the pathogenesis of SLE. In addition, the most significant module was obtained from the aforementioned PPI network (Fig. 1F). A total of 30 upregulated molecules were identified as hub genes with degrees of ≥10 and a cluster score of 28.414, as determined by MCODE of Cytoscape.

RT-qPCR validation

To verify the reliability of the high-throughput sequencing data, several candidate mRNAs of interest were initially identified for further analysis, including IFI27, olfactomedin 4 (OLFM4), interferon-induced protein 44 (IFI44), solute carrier family 12 member 1 (SLC12A1), keratin 1 (KRT1) and tubulin β 2A class IIa (TUBB2A). The expression levels of these mRNAs were confirmed by RT-qPCR. The results demonstrated that IFI27, OLFM4 and IFI44 were upregulated and that SLC12A1, KRT and TUBB2A were downregulated in SLE patients relative to those in the healthy group (all P<0.05; Fig. 2). These results demonstrated that the expression patterns of IFI27, OLFM4, IFI44, SLC12A1, KRT1 and TUBB2A were consistent with those obtained by high-throughput sequencing analysis.

Figure 2.

Relative expression levels of genes in SLE patients and healthy controls determined by reverse transcription-quantitative PCR. (A) IFI27, (B) OLFM4 and (C) IFI44 were significantly upregulated in SLE patients, whilst (D) SLC12A1, (E) KRT1 and (F) TUBB2A were significantly downregulated in SLE patients compared with those in the healthy group. The expression levels of mRNA were normalized against GAPDH, which served as an internal control. *P<0.05, ***P<0.001. Ctrl, healthy controls; SLE, systemic lupus erythematosus patients; IFI27, interferon-α-inducible protein 27; OLFM4, olfactomedin 4; IFI44, interferon-induced protein 44; SLC12A1, solute carrier family 12 member 1; KRT1, keratin 1; TUBB2A, tubulin β 2A class IIa.

Expression of CXCL10 mRNA in monocytes

Through Bioinformatics analysis of the high-throughput RNA-seq profiling data, it was determined that CXCL10 was upregulated by 6.27-fold in PBMCs from SLE patients compared with those in controls (P<0.01). In the mRNA PPI network, CXCL10 was one of the hub genes and significantly interacted with interferon induced with helicase C domain 1, 2′-5′-oligoadenylate synthetase (OAS)-like, OAS3, IFI27, receptor transporter protein 4, HECT and RLD domain containing E3 ubiquitin protein ligase family member 6 (HERC6), ubiquitin specific peptidase 18 (USP18) and HERC5. Since CXCL10 appeared to contribute to the pathogenesis of numerous autoimmune diseases, the expression of CXCL10 was further assessed in a total of 66 patients with SLE and 42 controls, with demographic data and SLEDAI scores listed in Table I.

It is well known that in the PBMC population, Th cells and B lymphocytes have pivotal roles in the initiation and progression of SLE (26). Therefore, CXCL10 mRNA was measured not only in PBMCs but also in Th cells and B lymphocytes. The results, presented in Fig. 3A, demonstrated that the expression of CXCL10 was significantly increased in PBMCs (P<0.001) and was similar in B lymphocytes (P<0.001), yet unchanged (P=0.881) in Th cells when comparing SLE patients with normal controls. In the present study, CXCL10 was identified as a DEG at the transcriptional level in the screen as well as in the verification experiments using RT-qPCR analysis. The trend of increased expression present in PBMCs from SLE patients was consistent with the results of the RNA-seq. Of note, the expression of CXCL10 in the Th-cell subpopulation was unchanged.

Figure 3.

Expression of CXCL10 and SNPs in CXCL10. (A) mRNA expression of CXCL10 in PBMCs, Th cells and B lymphocytes of SLE patients and healthy control. (B) SNP in the 3′-untranslated region of CXCL10 genotyped by PCR sequencing. The long horizontal line in the middle represents the mean and the bars indicate the standard deviation. n.s., no significance; ***P<0.001. Ctrl, healthy controls; SLE, systemic lupus erythematosus patients; SNP, single nucleotide polymorphism; CXCL10, chemokine C-X-C motif ligand 10; Th cell, T-helper cell; PBMCs, peripheral blood mononuclear cells.

Allele-specific expression of CXCL10

SNPs are considered to have a vital role in the function of CXCL10, which contributes to the pathogenesis of several diseases, including Graves' disease (27). Therefore, all of the DNA specimens of the present study were used to detect polymorphisms at selected SNP loci. The results indicated that when assessing the samples for rs35795399 (C/G/T), rs58658570 (A>T), rs34836828 (−>G) or rs148141229 (A>G), no mutation was identified in any of them, as demonstrated in Table I and Fig. 3B.

Discussion

Epidemiological surveys have indicated a five-year survival rate for SLE of 93.8% in China (28,30), and the prevalence in Asian patients is 2–3-fold higher than that in Caucasians; these rates in Asian patients are accompanied by serious clinical manifestations and higher fatality rates (1,28,30). Strong UV light unique to the plateau in Yunnan province may have an effect on the incidence and mortality associated with SLE (6,29,30).

Therefore, the high-throughput mRNA-seq profiles of PBMC from SLE patients from Yunnan province were first compared with those of health controls. As expected, the results indicated that the DEGs were enriched in autoimmune-associated pathways, including SLE, Rheumatoid arthritis and RIG-I-like receptor signaling pathway, as shown in Fig. 1D. In addition, GO and KEGG analysis indicated that the DEGs were enriched in viral infection, including defense response to virus and influenza A, as well as herpes simplex infection. Previous studies reported the roles of viruses in the pathogenesis of autoimmune diseases as epitope spreading, molecular mimicry, cryptic antigens and bystander activation (30). Viral infection has been frequently associated with a high level of IFN in most lupus patients (31). This was consistent with differentially expressed IFN-associated genes (IFI27, IFI44, CXCL10, IFI44L and sialic acid-binding Ig-like lectin 1) that were identified in the present study.

In past decades, RNA-seq, an innovative technology for comprehensive transcriptome profiling on a genomic scale, performed in a high-throughput and quantitative manner, has provided clues for identifying aberrantly expressed RNAs. First, the present analysis regarding the pathogenesis of SLE indicated that the enriched DEGs had an important role in the IFN signaling which is consistent with the results of Becker et al (32). This result has also been confirmed by a large number of studies for numerous years (33). Furthermore, IFN-α is an essential factor that acts against viral infection in innate immunity (34,35). Consequently, it is to be expected that the enriched molecules and pathways are linked to viral infection and initiation of an inherent anti-viral immune response. Finally, a variety of autoimmune cellular components and immune-associated pathways have also been identified in present study, which are identical to previous findings from patients with SLE (13,31,32). In summary, the conclusions of the present Bioinformatics analysis were highly consistent with the results of previous studies, including relatively high expression of IFI44, IFI44L, USP18, IFN-induced protein with tetratricopeptide repeats 3 and OAS3 (32,36). However, discrepancies with previous results were also present, including the upregulation of CXCL10 and downregulation of proline-rich 12 in the present study (32,36).

Elevated CXCL10 levels in serum has been considered to be an important potential biomarker in predicting SLE flares and making overall clinical decisions (30,37–39), with a particular role in the diagnosis of arthritis (30). However, to the best of our knowledge, few studies have investigated its involvement in the molecular mechanisms of SLE. Based on the dysregulation identified from the RNA-seq results of the present study, RT-qPCR was used and confirmed increased CXCL10 in PBMCs from SLE patients compared to controls. Recently, Lee et al (40) reported that IP-10/CXCL10 expression was 1.5-fold higher in the plasma of SLE patients relative to that in a healthy group. The fold change was even greater in the present study, demonstrating an up to 6.28-fold upregulation in SLE vs. Ctrl. It may be speculated that this is associated with the strong UV light in Yunnan. UV light, as an important initiation factor triggering the mechanisms of SLE, may activate an innate and subsequent adaptive immune response through Toll-like receptor or Th lymphocytes (6,29,30). In addition, CXCL10 in serum and urine increases in patients with SLE, depending on the level of disease activity and of anti-double stranded DNA antibody (41). In the present study, a correlation analysis of CXCL10 mRNA and disease activity was performed. The results indicated no difference in CXCL10 mRNA levels between patients with active and inactive SLE (data not shown). Furthermore, CXCL10 mRNA levels were not associated with any clinical indicators (data not shown). Therefore, it is likely that CXCL10 may be involved in the pathogenesis of SLE through other mechanisms at the transcriptional level, including cAMP-dependent protein kinase A signaling pathways, mitogen-activated protein kinase/PI3K and JNK/STAT (42). It has been reported that CXCL10 was able to polarize effector type I Th (Th1) and Th17 cells, whose imbalance contributes to the symptoms of SLE (30,43).

In the present study, CXCL10 mRNA expression in CD4+ Th cells and CD19+ B cells was compared between SLE patients and health controls for the first time, to the best of our knowledge. The absence of a significant difference of CXCL10 in Th cells between the two groups was not expected, but it may be explained by multiple antagonistic CD4+ lymphocyte subpopulations. It is well-known that the balance of Th1/Th2 and Th17/T-regulatory (Treg) cells as a large regulatory network is important in biological processes, and opposite alternations of CXCL10 in these cells may lead to unchanged total expression. Further research by our group will investigate possible causes. Furthermore, CXCL10 expression in the B lymphocytes of SLE patients was increased compared with that in controls. This result was similar to previous results in immature dendritic cells (DCs) (30,44). The expected result may be due to the UV light, as previously indicated (5,45), particularly when directly considering DCs (46). As a result, enhanced expression of CXCL10 may contribute to the interaction of immune and inflammatory responses and may have a vital role in the activity of SLE.

Finally, potential factors that may facilitate the expression of CXCL10 were explored in the present study. A previous study from Germany reported that polymorphisms in the CXCL10 gene may change the chemotaxis of activated T cells to targets and directly influence the severity of Grave's disease (27). Furthermore, polymorphisms in CXCL10 determine the secretion of the chemokine by monocyte-derived DCs upon exposure to Aspergillus fumigatus, and are closely associated with the risk of invasive aspergillosis after allogeneic stem-cell transplantation (44). Numerous studies on the function of genes have focused on the binding of micro (mi)RNA with SNPs of their target gene (47,48). In addition, it is well known that miRNAs inhibit the translation of their corresponding mRNA by directly binding to their 3′UTR, so it may be hypothesized that SNPs at the 3′UTR of CXCL10 are associated with its expression. Thus, four potential SNPs within the 3′UTR of CXCL10 we examined in the present study, but no mutation was identified in these polymorphic sites. However, correlations of CXCL10 gene promoter polymorphisms with the protein in plasma and disease susceptibility have been reported in infectious diseases, including malaria (49) and hepatitis B virus (50). Similar research on CXCL10 will be performed by our group in the future.

In conclusion, the present study indicated that CXCL10, a pivotal node in transcriptional gene expression networks, synergistically functions with other genes in the pathogenesis of SLE. It was revealed that CXCL10 mRNA is highly and significantly expressed in the PBMCs and B lymphocytes of SLE patients compared with that in controls, while its expression is unaltered in Th cells. However, no mutation was identified within four SNPs in the 3′UTR of CXCL10. Future studies on the expression of CXCL10 in other subgroups of Th cells, including Th1, Th2, Th17 and Treg cells, are required in order to better understand the role of CXCL10 in the mechanisms of SLE, or of SNPs in the promoter region of CXCL10.