Microarray profile of B cells from Graves' disease patients reveals biomarkers of proliferation.

B lymphocytes are the source of autoantibodies against the thyroid-stimulating hormone receptor (TSHR) in Graves' disease (GD). Characterization of autoimmune B-cell expression profiles might enable a better understanding of GD pathogenesis. To reveal this, the expression levels of long noncoding RNAs (lncRNAs) and mRNAs (genes) in purified B cells from patients with newly diagnosed GD and healthy individuals were compared using microarrays, which elucidated 604 differentially expressed lncRNAs (DE-lncRNAs) and 410 differentially expressed genes (DEGs). GO and pathway analyses revealed that the DEGs are mainly involved in immune response. A protein-protein interaction network presented experimentally validated interactions among the DEGs. Two independent algorithms were used to identify the DE-lncRNAs that regulate the DEGs. Functional annotation of the deregulated lncRNA-mRNA pairs identified 14 pairs with mRNAs involved in cell proliferation. The lncRNAs TCONS_00022357-XLOC_010919 and n335641 were predicted to regulate TCL1 family AKT coactivator A (TCL1A), and the lncRNA n337845 was predicted to regulate SH2 domain containing 1A (SH2D1A). TCL1A and SH2D1A are highly involved in B-cell proliferation. The differential expression of both genes was validated by qRT-PCR. In conclusion, lncRNA and mRNA expression profiles of B cells from patients with GD indicated that the lncRNA-mRNA pairs n335641-TCL1A, TCONS_00022357-XLOC_010919-TCL1A, and n337845-SH2D1A may participate in GD pathogenesis by modulating B-cell proliferation and survival. Therefore, the identified lncRNA and mRNA may represent novel biomarkers and therapeutic targets for GD.

Introduction

Graves’ disease (GD) is the most common cause of hyperthyroidism and diffuse goiter, affecting 1.0–1.5% of the population (1). GD is characterized by lymphocytic infiltration and the production of autoantibodies directed against components of the thyroid, including thyroglobulin (TG), thyroid peroxidase (TPO), and thyroid-stimulating hormone receptor (TSHR) (2, 3). TSHR autoantibodies mimic the function of TSH and lead to uncontrolled thyroid hormone production and hyperthyroidism (4, 5). B cells are the source of pathognomonic activating autoantibodies (TRAb) against TSHR in GD (6, 7). B cells also present antigens to T cells (8). Although B cells are essential for GD pathogenesis, little is known regarding the expression profiles of autoimmune B cells during GD development.

Long noncoding RNAs (lncRNAs), defined as non-protein-coding transcripts longer than 200 nucleotides, have been primarily studied for their effects on development, cancer, and cell differentiation (9, 10). There is emerging evidence that lncRNAs play crucial roles in immune regulation and autoimmunity. T lymphocytes, B lymphocytes, macrophages, dendritic cells, and NK cells all express lncRNAs, which have been shown to be involved in the differentiation and activation of immune cells (11). Accumulating evidence suggests that lncRNA dysregulation is involved in the pathogenesis of autoimmune diseases, such as systemic lupus erythematosus, rheumatoid arthritis, Crohn’s disease, multiple sclerosis, and autoimmune thyroid disease (AITD) (11, 12, 13, 14, 15). AITD, which includes GD and Hashimoto’s thyroiditis, is caused by an immune response to self-thyroid antigens (16). Single-nucleotide polymorphisms (SNPs) in the promoter of SAS-ZFAT, an antisense transcript of the ZFAT gene, were reported to determine susceptibility to AITD. SAS-ZFAT is exclusively expressed in CD19+ B cells in peripheral blood, indicating that SAS-ZFAT might control susceptibility to AITD by modulating B-cell function (12, 16). Recent studies have illustrated that lncRNAs can regulate the expression of protein-coding genes through cis-acting mechanisms on neighboring genes or through trans-acting mechanisms involving protein complexes such as transcription factors, splicing factors, or RNA-decay machinery (17, 18). In this study, we aimed to identify lncRNA–mRNA pairs that play a functional role in GD.

We found that patients with newly diagnosed GD had significantly higher numbers of B cells among their peripheral blood mononuclear cells (PBMCs) than healthy individuals. To reveal the mechanism of the increase in peripheral B cells, we studied the lncRNA and mRNA profiles of peripheral B cells from patients with GD by microarray analysis. We focused on lncRNA–mRNA pairs that were abnormally regulated in GD, especially those involving protein-coding genes that were related to cell proliferation. Among the identified genes, TCL1 family AKT coactivator A (TCL1A), and SH2 domain containing 1A (SH2D1A) have been reported to be associated with B-cell proliferation and survival. Using computational algorithms, we predicted that the lncRNAs TCONS_00022357-XLOC_010919 and n335641 regulate TCL1A and that the lncRNA n337845 regulates SH2D1A. The expression of both target genes was verified by qRT-PCR. The results suggest that the lncRNA–mRNA pairs n335641–TCL1A, TCONS_00022357-XLOC_010919–TCL1A, and n337845–SH2D1A might be useful biomarkers related to the proliferation of B lymphocytes. Furthermore, drugs targeting these genes might represent a novel therapy for GD.

Materials and methods

GD patients and healthy individuals

Peripheral blood was obtained from patients with newly diagnosed GD without thyroid associated ophthalmopathy (n = 34) and age-matched healthy individuals (n = 34). The patients with GD were clinically diagnosed with recent-onset GD in the Fifth People’s Hospital of Shanghai medical ward from 2014 to 2016 on the basis of the presence of overt hyperthyroidism, diffuse goiter or normal thyroid volume on ultrasonography, and a high rate thyroid iodine uptake. The characteristics of the patients are summarized in Table 1. Healthy controls normalized for sex, age, smoking status, thyroid function, and thyroid autoantibodies levels were recruited from the physical checkup center of the Fifth People’s Hospital of Shanghai, Shanghai, China. The controls were all in good health and had no family or personal history of thyroid disease. The exclusion criteria were a history of other autoimmune diseases, pregnancy during the study period, or a history of immunosuppressive medication. All procedures were in accordance with the ethical standards of the institutional and national research committees and the Declaration of Helsinki. The study protocols were approved by the Medical Ethics Committee of the Fifth People’s Hospital of Shanghai, Fudan University, Shanghai, China (NO. 002). All participants provided written informed consent.

Table 1

Clinical data of GD patients.

Subject	Age range (years)	Duration of clinical manifestation (days)	Hormone			Autoantibody
			FT3 (pmol/L)	FT4 (pmol/L)	TSH (mIU/L)	TGAb (IU/mL)	TPOAb (IU/mL)	TRAb (IU/L)
1	60–65	30	37.4	>100.0	<0.005	1364	23	8.11
2	35–40	60	13.0	>100.0	<0.005	18	97	13.87
3	25–30	150	4.9	1.5	0.050	0	67	12.50
4	60–65	5	17.9	52.5	<0.005	13	5	2.97
5	45–50	240	32.4	>100.0	0.005	195	>600	>40.00
6	40–45	2	22.4	71.2	<0.005	229	274	8.11
7	25–30	120	17.8	57.6	<0.005	446	539	6.20
8	50–55	60	21.7	48.9	<0.005	582	105	1.52
9	35–40	150	41.1	>100.0	<0.005	473	>600	6.38
10	40–45	30	46.7	>100.0	<0.005	229	274	17.84
11	45–50	420	22.3	47.2	<0.005	337	406	>40.00
12	20–25	13	30.3	83.3	<0.005	39	23	12.20
13	45–50	2	>50.0	>100.0	<0.005	610	>600	8.10
14	40–45	15	10.4	26.3	0.010	488	180	4.52
15	30–35	1	8.1	21.5	<0.005	>4000	343	29.45
16	50–55	90	22.0	58.0	<0.005	100	216	7.56
17	55–60	2	32.6	76.9	<0.005	68	342	20.04
18	60–65	10	19.4	43.5	<0.005	705	102	13.78
19	45–50	1	16.7	44.6	0.008	408	>600	>40.00
20	60–65	15	13.9	63.5	<0.005	524	56	8.11
21	35–40	30	>50.0	>100.0	<0.005	263	>600	35.76
22	40–45	30	22.5	71.2	<0.005	446	325	8.62
23	40–45	1	22.5	84.9	<0.005	229	364	18.90
24	10–15	60	32.3	92.8	0.010	29	157	11.07
25	45–50	75	30.6	71.2	<0.005	229	>600	8.11
26	45–50	30	34.5	90.0	<0.005	10	152	14.34
27	40–45	30	22.8	63.2	<0.005	113	595	4.68
28	30–35	8	6.9	15.6	0.010	23	121	4.00
29	20–25	25	30.3	83.3	<0.005	29	23	7.94
30	45–50	210	11.2	36.5	<0.005	778	>600	1.77
31	55–60	600	22.5	54.8	<0.005	497	>600	2.77
32	30–35	100	40.0	87.6	<0.005	118	410	10.70
33	40–45	30	41.3	>100.0	<0.005	564	8	12.72
34	30–35	2	19.8	79.5	<0.005	132	333	4.26

Reference range: FT3, 3.1–6.8 pmol/L; FT4, 12.0–22.0 pmol/L; TSH, 0.270–4.20 mIU/L; TGAb <115 IU/mL; TPOAb <34 IU/mL; TRAb <1.75 IU/L.

FT3, free triiodothyronine; FT4, free thyroxine; TGAb, thyroglobulin antibody; TPOAb, thyroid peroxidase antibody; TRAb, TSH receptor antibody; TSH, thyroid stimulating hormone.

Thyroid function and antithyroid antibodies detection

The levels of free triiodothyronine (FT3) (reference range, 3.1–6.8 pmol/L), free thyroxine (FT4) (reference range, 12.0–22.0 pmol/L), thyroid stimulating hormone (TSH) (reference range, 0.270–4.20 mIU/L), thyroglobulin antibody (TGAb) (cut-off level, 115 IU/mL), thyroid peroxidase antibody (TPOAb) (cut-off level, 34 IU/mL), and TSH receptor antibody (TRAb) (cut-off level, 1.75 IU/L) were measured automatically by Cobas c8000 machine via utilizing a commercial electro-chemiluminescence immunoassay kit (Roche Diagnostics, Cobas 8000 analyzer e602 module). The intra-assay CVs of FT3, FT4, TSH, TGAb, TPOAb, and TRAb were 1.4–6.5%, 2.6–4.3%, 1.1–3.0%, 1.3–4.9%, 2.7–6.3%, and 0.9–7.6%, respectively, and the inter-assay CVs of coefficients of these indexes were 1.9–7.2%, 5.1–8.4%, 3.3–7.2%, 3.4–6.3%, 4.2–9.5%, and 1.9–11.4%, respectively.

Isolation of peripheral blood mononuclear cells (PBMCs)

Discontinuous Lymphoprep (Axis-Shield PoC As, Oslo, Norway) gradient centrifugation was used to collect PBMCs, which were resuspended in PBS containing 2% fetal bovine serum (FBS).

Flow cytometry analysis

Isolated PBMCs were blocked by Fc receptors before staining with a specific antibody. The PE-conjugated anti-CD19 antibody (Biolegend, San Diego, CA, USA) was used for surface staining following the manufacturer’s recommended protocol. Flow cytometry was conducted using a CyAn™ ADP Analyzer (Beckman, Coulter, Inc., Carlsbad, CA, USA) and data were analyzed using the FlowJo software (Tree Star, Ashland, Or, USA).

B-cell purification

Human CD19+ B cells were isolated by positive selection using PE-conjugated anti-CD19 (Biolegend) and anti-PE microbeads (Miltenyi Biotec, Bergisch Gladbach, Germany) according to the manufacturer’s instructions, and the purity was confirmed to be over 90.0% as judged through flow cytometry (Supplementary Fig. 1, see section on supplementary materials given at the end of this article).

RNA extraction and purification

RNA was extracted from the purified CD19+ B cells from 21 healthy individuals and 24 GD patients (patients 11–34, Table 1) using the TRIzol Reagent (Life technologies) as per the manufacturer’s protocol. The purity and concentrations of RNA samples were measured through NanoDrop ND-2000 spectrophotometry (Thermo Scientific). For Microarray experiments, total RNA was further purified using the RNeasy micro kit (QIAGEN) and RNase-Free DNase Set (QIAGEN).

lncRNA and mRNA expression profiling

CD19+ B cells from the peripheral blood of three GD patients and three healthy individuals were analyzed using microarrays. The RNA from the three GD patients (patients 11–13, Table 1) and the three healthy individuals were pooled separately. lncRNA and mRNA expression profiling was performed using Affymetrix GeneChip Human Transcriptome Array 2.0 microarrays (Affymetrix). Total RNA was amplified, labeled, and purified using the Affymetrix WT PLUS Reagent Kit (Affymetrix) to obtain biotin-labeled complimentary DNA (cDNA) according to the manufacturer’s instructions. The labeled cDNA was hybridized and washed using a GeneChip Hybridization, Wash, and Stain Kit in a Hybridization Oven 645 and a Fluidics Station (all from Affymetrix) according to the manufacturer’s instructions.

Data acquisition

Arrays were scanned in an Affymetrix GeneChip Scanner 3000 (Affymetrix). The Command Console Software (Affymetrix) was used to control the scanner and summarize probe cell intensity data (CEL file generation) with default settings. Then, raw data were normalized by the Expression Console.

Enrichment analysis of DEGs

To identify differential gene expression between healthy individuals and GD patients, a |Fold Change (linear)| ≥2 was considered as differently expressed. Gene Ontology (GO, http://www.gneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg) analyses were conducted on these differentially expressed genes (DEGs) to categorize the biological functions and pathways that they are associated with.

lncRNA target prediction

To identify the cis- or trans-regulatory targets of lncRNAs, two independent algorithms were used. The first algorithm searched for cis-regulated target genes with the help of gene annotations from UCSC (http://genome.ucsc.edu/). Genes located within 10 kb upstream or downstream of the lncRNAs were considered cis-regulated target genes. The second algorithm was based on mRNA sequence complementarity and prediction of RNA duplex energy. After a first round of screening using BLAST software, the RNAplex software was used to identify trans-regulated target genes (19). The RNAplex parameters were set as in previous studies (19).

Construction of a protein–protein interaction (PPI) network

To construct a PPI immune network involving differentially expressed genes (DEGs), all possible combinations of gene pairs were defined using InnateDB (http://www.innatedb.ca/), a publicly available database of the genes, proteins, experimentally verified interactions, and signaling pathways involved in the innate immune response (20). A final interaction network was constructed using the force-directed layout in the Cytoscape software. Differently expressed lncRNAs and their associated protein-coding genes were identified on the basis of chromosomal localization and BlAST sequence alignment.

Quantitative real-time polymerase chain reaction (qRT-PCR) validation

Total RNA from each sample was reverse transcribed into cDNA (Takara Biotechnology Co.). The cDNA was used in qRT-PCR assays to validate the presence of differentially expressed mRNAs and lncRNAs in B cells from GD patients. Real-time PCR amplifications were performed on an ABI 7500 thermocycler in 20 μL reaction volumes containing cDNA, primers, and SYBR Green Real-time PCR Master Mix (Toyobo Co., Ltd., Osaka, Japan). The primer sequences are listed in Table 2. Thermocycling was performed as follows: initial denaturation at 95°C for 30 s; 45 cycles of denaturation at 95°C for 5 s, and elongation at 60°C for 15 s. The specificity of the PCR reaction was controlled by assessing the melting curves. The ΔΔCt method was used to compute the target gene expression relative to that of ACTB. All primers were synthesized and purified by Shanghai Sangon Biotechnology Co. Ltd. (Shanghai, China).

Table 2

Primer information.

Primer	Sequence (5′–3′)
BIRC3-F	TTGCTGTGATGGTGGACTCAGGTGTTG
BIRC3-R	GGTAACTGGCTTGAACTTGACGGATGAACT
EAF2-F	GCCTTCCACACTGTGCGCTATGACTTC
EAF2-R	GTCACCTGTTCACCTTCACCAACCTCAAG
IL7R-F	AAGGCTTCTGGAGTGAATGGAGTCCAAGTT
IL7R-R	GCCAAGATGACCAACAGAGCGACAGAGA
FYN-F	CGTCTTTGGAGGTGTGAACTCTTCGTCTCA
FYN-R	AACTCAGGTCATCTTCTGTCCGTGCTTCAT
TCL1A-F	ACATCAAGATTGACGGCGTGGAGGACAT
TCL1A-R	AAGGAGACAGGTGCTGCCAAGACATCA
SH2D1A-F	AAATCAGCAGGGAAACCGGCGAGAAG
SH2D1A-R	ACCGTGATACAGCACACATAGGCAGTACA
n335641-F	AAGCCTCTGGGTCAGTGGT
n335641-R	CAACCTCTATCCCTTTTCCAC
TCONS_00022357-XLOC_010919-F	GGAGGCAGAAGCAGGAGAAT
TCONS_00022357-XLOC_010919-R	ACTGACCCAGCCTCATCAAT
n337845-F	AAAGGAAATGTGTTTGGCAGT
n337845-R	TCTTACTGATGCTTTTGCCAGA
ACTB-F	CTGGAACGGTGAAGGTGACA
ACTB-R	CGGCCACATTGTGAACTTTG

F, forward; R, reverse.

Statistical analysis

All statistical analyses were performed using SPSS 17.0 (SPSS) and PRISM 6.0 (GraphPad Software Inc.). Results were presented as the mean ± s.e.m. For comparison between two groups, Student’s unpaired t-test or unpaired t-test with Welch’s correction was used. P values <0.05 were considered statistically significant.

Results

Increased peripheral B lymphocyte levels in GD patients

Compared with the normal reference values, the levels of free triiodothyronine (FT3), free thyroxine (FT4), TRAb, TPOAb, and TGAb in the blood of GD patients were remarkably high, whereas that of TSH was strikingly low. The basic information about the patients is shown in Table 1. The average proportion of CD19+ B cells in the PBMCs of GD patients was 10.7%, which was significantly higher than that in the PBMCs of the healthy controls (Fig. 1).

Figure 1

The ratio of CD19+ B cells to PBMCs increased in GD patients. (A) Representative percentage of CD19+ B cells in healthy individuals and new-onset GD patients as determined by flow cytometry. (B) CD19+ B-cell frequencies in total PBMCs were measured in 34 healthy individuals and 34 new onset GD patients, with each point representing a single person. ***P < 0.001. The data shown are the mean ± s.e.m.

Differential mRNA and lncRNA expression profiles in B lymphocytes from GD patients

For revealing the mechanism of the abnormality of B cells in GD, the expression profiles of lncRNA and mRNA were investigated. In total, 604 lncRNAs were differentially expressed (DE-lncRNAs) between the purified B cells of healthy individuals and those of patients with GD. Of those, 303 were upregulated in GD patients relative to the healthy controls and the other 301 were downregulated in GD patients relative to the healthy controls (Fig. 2A). Furthermore, a total of 410 differentially expressed genes (DEGs) were also identified, 331 (80.7%) of which were downregulated in GD patients compared with the healthy controls and 79 (19.3%) of which were upregulated in GD patients compared with the healthy individuals (Fig. 2B). Our data showed that molecular events in peripheral blood B cells, such as lncRNAs and mRNAs, have changed during GD development.

Figure 2

Differentially expressed lncRNA and mRNA between healthy individuals and GD patients. (A) Scatterplot showing that the differentially expressed lncRNA, n335641, TCONS_00022357-XLOC_010919, and n337845 were differentially expressed in B cells of GD patients relative to the healthy individuals. (B) Scatterplot showing that the differentially expressed mRNA, TCL1A, and SH2D1A were differentially expressed in B cells of GD patients relative to the healthy individuals. Red dots, upregulated genes; Blue dots, downregulated genes.

Functional enrichment analysis

Gene ontology (GO) and KEGG pathway enrichment analyses were performed to identify the functions of the DEGs. GO project used to classify the significantly regulated genes according to biological process (BP), cellular component (CC), and molecular function (MF). The DEGs were significantly enriched with GO terms related to the immune response. The most common biological processes associated with the DEGs were immune response, T cell receptor signaling pathway, and regulation of immune response (Fig. 3A). Consistent with the enriched GO categories, the KEGG pathway analysis showed that the major signal-transduction pathways associated with the DEGs were related to immune signaling, such as T cell receptor signaling pathway, natural killer cell mediated cytotoxicity, and antigen processing and presentation (Fig. 3B).

Figure 3

GO and pathway enrichment analyses. (A) Significantly enriched GO terms in biological process clusters (red), cellular component clusters (green), and molecular function clusters (blue) were analyzed. (B) KEGG pathway analysis.

Confirmation of the microarray results by qRT-PCR

To confirm the results of the microarray analysis, two mRNAs that were upregulated (BIRC3 and EAF2) and two mRNAs that were downregulated (IL7R and FYN) in the B lymphocytes of the patients with GD were randomly selected and analyzed by qRT-PCR in 21 healthy individuals and 24 GD patients. As expected, the expression levels of BIRC3 and EAF2 were increased in GD patients relative to those in the healthy individuals, whereas the expression levels of IL7R and FYN were decreased in GD patients relative to those in the healthy individuals, consistent with the microarray results (Fig. 4).

Figure 4

Validation of microarray results by qRT-PCR. (A, B, C and D) The RNA of purified CD19+ B cells from another 21 healthy individuals and 24 GD patients were used to verify the randomly selected mRNA. BIRC3 and EAF2 were upregulated, whereas IL7R and FYN were downregulated significantly in purified B cells from GD patients as compared with healthy controls. The data shown are the means ± s.e.m. Three replicates showed similar results. ***P < 0.001.

Analysis of differentially expressed lncRNA and mRNA associations

Because the roles of lncRNAs are mainly executed through effects on protein-coding target genes, the cis-regulatory targets of DE-lncRNAs were identified on the basis of the distance between the lncRNAs and their neighboring protein‐coding genes (adjacent upstream or downstream) and the trans-regulatory targets of DE-lncRNAs were analyzed by RNAplex software. The results showed that 84 DE-lncRNAs were significantly associated with 81 differentially expressed protein-coding genes to generate 85 lncRNA–mRNA pairs. Of the 84 lncRNA candidates, 61 were identified as sense lncRNAs, 5 as antisense lncRNAs, 3 as bidirectional lncRNAs, 2 as intergenic lncRNAs, and 13 as other lncRNAs.

To identify lncRNA–mRNA associations and evaluate the impact of the DE-lncRNAs on gene expression, the predicted target genes of the DE-lncRNAs were identified and combined with the mRNA transcriptomes of the study participants. Because the patients with GD had higher percentages of peripheral B lymphocytes than the healthy controls, we hypothesized that the regulation of B-cell proliferation was abnormal in the patients with GD. Fourteen pairs of DE-lncRNAs and adjacent protein‐coding genes associated with cell proliferation were identified (Table 3). Among the identified genes, those encoding TCL1A (21, 22) and SH2D1A (23, 24, 25) have been reported to be involved in B-cell proliferation and survival. Thus, n335641–TCL1A, TCONS_00022357-XLOC_010919–TCL1A, and n337845–SH2D1A were identified as candidate lncRNA–mRNA pairs that were differentially regulated between GD patients and healthy individuals.

Table 3

Deregulated lncRNAs and their corresponding protein-coding genes (associated with cell proliferation).

lncRNA	Log2 (FC)-lncRNA	Type	mRNA	Log2 (FC)-mRNA
n326505	−1.21	antisense	CDC14A	−1.39
n332938	−1.23	sense	EIF2AK2	+1.18
n407800	−2.38	sense	HPGD	−1.06
n409270	−2.76	sense	TCF7	−1.07
TCONS_00010118-XLOC_004584	−1.09	other	NDFIP1	−1.90
n384644	−2.58	sense	IL6ST	−3.09
n335470	−3.84	sense	TNFAIP3	−3.72
n378743	−2.72	bidirectional	PRKCQ	−2.13
n335567	−1.58	other	IFITM1	−1.16
TCONS_00022357-XLOC_010919	+2.31	bidirectional	TCL1A	+1.05
n335641	+1.97	sense	TCL1A	+1.05
n383211	−1.21	other	ERN1	−1.20
n409231	−2.30	sense	TYROBP	−2.66
n337845	−1.50	intergenic	SH2D1A	−1.39

FC, fold change.

To validate the candidate lncRNA–mRNA pairs in independent patients with GD, the expression of the pairs was measured in an additional 21 healthy individuals and 24 patients with GD by qRT-PCR. As expected, the lncRNAs n335641 and TCONS_00022357-XLOC_010919 were upregulated, whereas n337845 was downregulated, in the purified peripheral B cells of the patients with GD relative to the healthy controls (Fig. 5A, B and C). Likewise, TCL1A, the predicted target of n335641 and TCONS_00022357-XLOC_010919, was upregulated, whereas SH2D1A, the predicted target of n337845, was downregulated in the B cells of the patients with GD relative to the healthy controls (Fig. 5D and E), which were consistent with the results obtained from the bioinformatics analysis. These potential candidate lncRNAs (Fig. 2A) and mRNAs (Fig. 2B) are also shown in scatter plots. However, the roles and mechanisms of these lncRNA–mRNA pairs in dysfunction of B cells in patients with GD were not clear and need to be studied in depth.

Figure 5

Validation of mRNA and lncRNA microarray results by qRT-PCR. (A, B, C, D, and E) The RNA of purified CD19+ B cells from another 21 healthy individuals and 24 GD patients were used to verify the significantly regulated lncRNA and their related target mRNA. (A and B) lncRNA n335641 and TCONS_00022357-XLOC_010919, which potentially regulate TCL1A were upregulated. (C) lncRNA n337845, which potentially regulates SH2D1A was downregulated. (D and E) TCL1A mRNA was upregulated, whereas SH2D1A was downregulated significantly in purified B cells from GD patients as compared with healthy controls. The data shown are the means ± s.e.m. Three replicates showed similar results. ***P < 0.001.

Correlation analysis for DEGs and DE-lncRNAs

To reveal the PPI networks of the DEGs, immunologically focused network analyses of 126 DEGs (13 upregulated and 113 downregulated in the patients with GD relative to the controls) and 193 potential interactions were performed using the InnateDB platform (Fig. 6). The downregulated and upregulated DEGs with high degrees are illustrated in Table 4. Accordingly, the downregulated DEG SH2D1A had a connective degree of 7, potentially interacting with FYN (26, 27, 28, 29, 30), LCK (28), CD247 (31), CD5 (32) and so on, and the upregulated DEG TCL1A had a connective degree of 2, potentially interacting with JUNB (33) and FOS (33). To further uncover the correlations between the DEGs and the DE-lncRNAs, the potential functions of the DE-lncRNA–protein pairs were identified. In total, there were 28 regulatory relationships between DEGs and DE-lncRNAs, including 27 DEGs and 28 DE-lncRNAs (3 upregulated in GD and 25 downregulated in GD). For instance, n335641 and XLOC_010919 were predicted to regulate TCL1A, whereas n337845 was predicted to regulate SH2D1A (Fig. 6).

Figure 6

The integrated network formed by differentially expressed genes annotated in InnateDB and lncRNAs. Rounded nodes represent genes and triangular nodes represent lncRNAs. Red nodes represent the upregulated ones, and green nodes represent downregulated ones. Black lines represent protein–protein interaction (PPI). Dashed lines represent lncRNA–mRNA interaction.

Table 4

The downregulated genes with the connectivity degree ≥6 and upregulated genes with the connectivity degree ≥2 in the protein–protein interaction network.

Downregulated DEGs	Connective degree	Downregulated DEGs	Connective degree	Upregulated DEGs	Connective degree
FYN	22	LGALS1	8	BIRC3	4
ZAP70	18	GAPDH	7	NSF	3
IL7R	18	HSPA5	7	IFIT3	3
LCK	17	SH2D1A	7	IFIT1	3
ARRB1	13	CD5	7	STX7	3
LCP2	12	TNFAIP3	6	TCL1A	2
CD247	10	CD2	6	XAF1	2
HSPA8	10	CD3E	6
CD4	9	TUBB	6
ATXN1	8	ITK	6

Discussion

In GD, B cells are the source of TRAb against TSHR (2, 6, 34). It is widely accepted that the quantity and function of B cells are anomalous in GD patients (7, 8, 35). We previously reported that B cells infiltrate the thyroid tissues of patients with GD (36). In the present study, we found that the percentage of peripheral B cells was higher in GD patients than in healthy individuals; however, the molecular determinants of that difference remain unknown. Microarray analysis identified 604 lncRNAs and 410 mRNAs with differential expression between the B cells of patients with GD and those of healthy individuals. The differential expression of a subset of those lncRNAs and mRNAs was confirmed by qRT-PCR. GO and KEGG pathway analyses of the DEGs showed enrichment of biological functions and signaling pathways involved in the immune response, such as T-cell receptor signaling pathway, regulation of immune response, and antigen processing and presentation. The PPI network constructed by InnateDB (experimentally validated interactors) showed that the mRNA expression of most of the DEGs was decreased in GD compared with that in healthy individuals. To our knowledge, no studies have examined lncRNA and mRNA expression profiles of peripheral B cells in GD.

The functions of the DE-lncRNAs might provide hints about the mechanisms of B cells dysfunction and pathogenesis in GD. Previous research reported that lncRNAs can regulate the expression of neighboring protein-coding genes via cis-acting mechanisms (37), and we predicted that some of the DE-lncRNAs in GD would function based on a cis-acting mechanism. Because the functions of most of the DE-lncRNAs were not known, we predicted them on the basis of the positional relationships between the lncRNAs and nearby genes. lncRNAs can also regulate the expression of protein-coding genes via trans-acting mechanisms including interactions with chromatin modifiers, transcription factors, splicing factors, or RNA-decay machinery (38, 39, 40). Owing to this, 85 lncRNA–mRNA pairs were identified between patients with GD and healthy controls.

In our study, we found the percentage of peripheral B lymphocytes was increased in GD patients. Previous research also reported that the percentage and absolute count of peripheral B lymphocytes were markedly increased in GD patients than in normal controls, but the proportions and absolute counts of B lymphocytes were reduced to normal in euthyroid patients with GD after antithyroid drug therapy (41). Thus, we hypothesized that abnormal regulation of proliferation was the cause of abnormal B-cell function in GD. Therefore, we focused on protein-coding genes linked to cell proliferation. After combining the lncRNA and mRNA expression data, we identified 14 pairs of DE-lncRNAs and corresponding protein-coding genes associated with cell proliferation (Table 3). Among those protein-coding genes, TCL1A acts as co-activator of AKT kinase and thus promotes AKT-induced cell survival and proliferation (21), which are part of the normal development of early B cells and T cells (42). TCL1A has also been reported to regulate cell survival and cell death pathways in human naive and IgM memory B cells (22). The PPI network shows that TCL1A is related to the pathogenesis of B-cell chronic lymphocytic leukemia (B-CLL) via physical interactions with JUNB and FOS (33), which play important roles in the regulation of thyroid cell function (43, 44, 45). That suggests that TCL1A might be involved in thyroid autoimmune disorders through interactions with JUNB and FOS. In addition, X-linked lymphoproliferative disease 1 (XLP1), manifested as lymphomas resulting from uncontrolled B-cell proliferation upon infection with Epstein-Barr virus, is a primary immunodeficiency caused by mutations or deletions of SH2D1A gene, which encodes the SLAM-associated protein (SAP) (46, 47). B-lymphoblastoid cell lines (B-LCLs) from XLP1 patients have an intrinsic defect that affects cell activation, apoptosis, and proliferation (23), and SAP-deficient mice also have a marked defect in B-cell proliferation (24, 48). SAP could form a ternary complex with the Lyn kinase and the inhibitory IgG Fc receptor FCGR2B to regulate B-cell proliferation and survival, leading to a decrease in B-cell proliferation and an increase in activation-induced cell death (25). Among the proteins that interacted with SH2D1A in the PPI network, LCK and CD247 were previously shown to be hypermethylated in CD4+ T cells of GD patients (49). Furthermore, CD247 is a susceptibility gene in AITD (50). CD5 peripheral B lymphocytes were increased in GD and were associated with disease activity (51, 52, 53), indicating that SH2D1A may also play roles in GD pathogenesis via its interactors. The lncRNAs TCONS_00022357-XLOC_010919 and n335641 were predicted to target TCL1A, and the lncRNA n337845 was predicted to target SH2D1A. The differential expression of both target genes in GD was confirmed by qRT-PCR. The results indicate that the n335641–TCL1A, TCONS_00022357-XLOC_010919–TCL1A, and n337845–SH2D1A lncRNA–mRNA pairs might be novel biomarkers associated with the proliferation of B lymphocytes in GD. Further study should be conducted to verify whether those are biomarkers for the diagnosis, treatment, and prognosis of GD.

In summary, we analyzed lncRNA and mRNA expression profiles in B cells to identify molecules that influence GD pathogenesis by regulating B-cell function, proliferation, and survival. Our study is preliminary and therefore has some limitations. First, we only isolated B lymphocytes from peripheral blood and not from the thyroid or local neck lymph nodes surrounding the thyroid. B lymphocytes of different origins might contain different genetic information. Second, the lncRNA database is updated continuously, so the datasets used in our study might be limited. Therefore, the lncRNA signatures might need further identification. Furthermore, the relationships between lncRNAs and mRNAs predicted by the two algorithms were not verified by experiments in vivo and in vitro. Third, the mechanisms by which the lncRNA–mRNA pairs influence B-cell proliferation have not been uncovered. Further study should be carried out to verify whether the lncRNA–mRNA pairs can be regarded as biomarkers for the diagnosis, treatment, and prognosis of GD. Except for three pairs of lncRNA–mRNA pairs mentioned previously, other differently expressed lncRNA–mRNA pairs may also be related to B-cell proliferation or other abnormal B-cell functions such as B-cell activity or suppression ability. We believe, however, that the lncRNA–mRNA expression pattern provides a basis for understanding the pathogenesis of GD and will lead to diagnostic biomarkers or therapeutic targets in future.

Conclusions

Our study provides an integrated analysis of lncRNA and mRNA expression profiles in peripheral B cells of GD patients. The results suggest that the lncRNA–mRNA pairs n335641–TCL1A, TCONS_00022357-XLOC_010919–TCL1A, and n337845–SH2D1A participate in GD pathogenesis by modulating B-cell proliferation and survival. The identified lncRNA–mRNA pairs should be investigated further as potential biomarkers and therapeutic targets in GD, which might provide a better understanding of the molecular basis of GD.

Declaration of interest

The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

Funding

This work was supported by the Natural Science Foundation of Shanghai (Grant Number: 18ZR1429900) and the National Natural Science Foundation of China (Grant number: 81302570, 81701488/H0422).

Data availability

All microarray raw data files are available from the Gene Expression Omnibus (GEO) database (accession number GSE136709), https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE136709.

Author contribution statement

J X C and Z B B designed the research. J X C, Z B B, and W Y H performed experiments. J X C and Y Q analyzed data. W W and L X Y revised the manuscript. J X C, W Y H, and H L Q wrote the manuscript, and Z B B and L J edited the manuscript.