Long non‑coding RNA expression profiles identify lncRNA‑XLOC_I2_006631 as a potential novel blood biomarker for Hashimoto's thyroiditis.

Long non‑coding RNAs (lncRNAs) have been increasingly recognized as important immune checkpoints involved in the pathogenesis of autoimmune diseases. However, the exact role of lncRNAs in Hashimoto's thyroiditis (HT) has been rarely studied. The aim of the present study was to investigate the role of lncRNAs and the potential biomarkers in HT, a total of 33 patients with HT and 32 healthy volunteers were enrolled in the present study, and five patients and five healthy controls were investigated using next generation sequencing. A total of 218 dysregulated lncRNAs, including 94 upregulated and 124 downregulated lncRNAs, were identified and examined in the peripheral blood mononuclear cells (PBMCs) from patients with HT. The majority of the lncRNAs were intergenic and exonic (66.06%). Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analysis demonstrated that abnormally expressed lncRNAs were enriched in the 'NF‑kB expression', in the 'TGF‑β signaling pathway' and in the 'JAK‑STAT signaling pathway', which are associated with the immunopathogenic mechanisms of HT. In total, three lncRNAs (LOC729737, XLOC_I2_006631 and BC041964) were validated and had a trend identical to that detected by the sequencing results. The expression of lncRNA‑XLOC_I2_006631 was upregulated and was positively correlated with the serum concentrations of anti‑thyroperoxidase antibody in patients with HT. Methyl‑CpG‑binding protein 2 (MECP2) was identified as the potential regulatory gene of lncRNA‑XLOC_I2_006631 using a prediction program. The expression of MECP2 was increased and was positively correlated with the elevated expression levels of lncRNA‑XLOC_I2_006631 and anti‑thyroperoxidase antibody in patients with HT. Furthermore, lncRNA‑XLOC_I2_006631 was able to regulate MECP2 expression in vitro. Receiver operating characteristic curve analysis suggested that lncRNA‑XLOC_I2_006631 has a potential diagnostic value. Collectively, the present results indicated the important role of dysregulated lncRNAs in HT and demonstrated that lncRNA‑XLOC_I2_006631 functioned as a positive regulator of MECP2 expression, suggesting a potential mechanism. Thus, lncRNA‑XLOC_I2_006631 may be used as a biomarker of HT.

Introduction

Hashimoto's thyroiditis (HT) is a common organ-specific autoimmune thyroid disease characterized by autoantibodies against thyroid tissue and lymphocytic infiltration in the thyroid parenchyma, resulting in the progressive destruction and atrophy of thyrocytes and a diffusely enlarged thyroid gland, which eventually evolves into hypothyroidism (1,2). HT was initially described by Dr Hakaru Hashimoto in 1912 as a new type of lymphomatous thyroid tumor in Japan (3), and was recognized as an autoimmune disease of the thyroid several years later (4). HT is currently the most frequently occurring autoimmune disease, with an incidence of ~1% per year (5). According to National Health and Nutrition Examination Survey of the United States, the prevalence of HT is estimated to be 5% (6). HT is more common in women compared with men, and in countries with iodine deficiency (7). The serological features of patients with HT are characterized by elevated production of anti-thyroglobulin-Ab (TgAb) and anti-thyroperoxidase-Ab (TPOAb) (8,9). Additionally, some patients have high concentrations of thyrotropin (10). HT is currently considered a complex and continuously expanding disease caused by genetic susceptibility, environmental factors, age, sex and immune disorders (11). However, the pathogenesis of HT remains unknown.

Long non-coding RNAs (lncRNAs) have been recently recognized as a new class of non-coding RNAs >200 nucleotides in size (12,13). lncRNAs are shorter than mRNA, have low evolutionary sequence conservation and are highly prevalent in the eukaryotic transcriptome (14). Based on various types of transcriptional regulation and the position in the genome relative to the protein-coding genes and enhancer regulatory elements, lncRNAs have been classified as intronic lncRNAs, antisense lncRNAs, long intergenic ncRNAs, enhancer RNAs and transcribed pseudogene lncRNAs (15). Multiple regulatory mechanisms of lncRNAs have been reported to influence transcription and post-transcriptional events, as well as the translation of diverse individual genes and the whole gene networks (16). The regulatory mechanisms of lncRNA vary depending on their location in the nucleus or cytoplasm (17). Numerous studies have reported that lncRNAs serve key roles in pathogenesis of various diseases, including cancer, cardiovascular diseases, nervous system disease and immune diseases (18-21). Certain lncRNAs are involved in autoimmune diseases, such as multiple sclerosis, rheumatoid arthritis (RA) and systemic lupus erythematosus (SLE) (22-24). However, the expression profiles and function of lncRNAs in patients with HT are yet to be elucidated.

The present study aimed to identify lncRNAs in the peripheral blood mononuclear cells (PBMCs) from patients with HT and to investigate the potential role of lncRNAs in the pathogenesis of HT using next generation sequencing (NGS). Using this approach, the present study aimed to elucidate the mechanism of lncRNAs and to identify the potential diagnostic value of lncRNAs in patients with HT.

Materials and methods

Subjects and samples

Patients with HT usually have no obvious clinical symptoms at the early stage, and most only identify the elevated serum levels of TgAb and TPOAb via physical examination. To obtain a definitive diagnosis, strict inclusion criteria are required. The samples were collected between December 2017 and May 2018 at The Affiliated People's Hospital of Jiangsu University. In total, 33 patients with HT were enrolled in this study. Among them, five female patients with HT were used for NGS. Moreover, 28 patients with HT, including 21 women and seven men, were used to verify the sequencing results. Their ages ranged from 27-75. The main clinical characteristics of these patients are listed in Table I.

Table I

Clinical features of patients with HT and healthy controls.

Variable	Patients with HT	Healthy controls	Range	P-values
Number	28	27
Sex (M/F)
Male	7	14
Female	21	13
Age, years	41±13	35±9		0.06
FT3, pmol/l	5.36±0.65	5.41±0.53	3.28-6.47	0.77
FT4, pmol/l	9.72±1.81	10.41±1.29	7.64-16.03	0.11
TSH, uIU/ml	4.95±4.54	1.73±0.54	0.56-5.91	<0.01
TgAb, IU/ml	195.1±247.6	0.4±0.6	0-4	<0.01
TPOAb, IU/ml	709.9±536.2.5	0.9±1.3	0-9	<0.01

Data are presented as the mean ± SD, and were compared using unpaired t-tests. FT3, free triiodothyronine; FT4, free thyroxine; TSH, thyroid stimulating hormone; TgAb, anti-thyroglobulin-Ab; TPOAb, anti-thyroperoxidase-Ab; HT, Hashimoto's thyroiditis.

All patients were diagnosed based on clinical manifestations and auxiliary examinations, including B-flow ultrasonic imaging and laboratory criteria. The serum concentrations of free triiodothyronine (3.28-6.47 pmol/l), free thyroxine (FT4; 7.64-16.03 pmol/l), thyroid stimulating hormone (TSH; 0.56-5.91 uIU/ml), TgAb (0-4 IU/ml) and TPOAb (0-9 IU/ml) were measured using an LDX-800 system (Beckman Coulter, Inc.) according to the manufacturer's instructions. All patients had positive tests for TgAb or TPOAb. A total of nine patients with HT and hypothyroidism had high levels of TSH with or without low levels of FT4, and other patients with euthyroid had physiological levels of TSH. In total, 32 age- and sex-matched healthy subjects were included as controls, including five healthy female volunteers for NGS and 27 healthy volunteers for validation of the sequencing results. All healthy subjects were euthyroid and free of thyroid-specific autoantibodies, as well as had no history of thyroid disease, chronic inflammatory diseases, tumors and other autoimmune diseases (25). All subjects showed no obvious clinical symptoms of infection and were not taking immunosuppressive drugs. The numbers of peripheral leukocytes in all individuals were within the normal range (3.50-9.50×109/l). After rigorous screening, five female patients and five female volunteers were selected at random and were involved in sequencing due to the significantly higher incidence of HT in women (10).

The HT group criteria included: i) All patients with HT were clinically diagnosed based on clinical manifestations and auxiliary examinations; ii) there was no other autoimmune diseases except HT in patients with HT; iii) all the individuals had not recently taken immunosuppressive drugs for other diseases; iv) they had not been infected recently; and v) the age ranged between 20-50. Healthy control group criteria were as follows: i) Thyroid function or thyroid B-ultrasound were normal in healthy subjects; ii) there was no occurrence of a recent infection; iii) all had no autoimmune diseases and had not recently taken immunosuppressive drugs; and iv) the age ranged between 20-50. Peripheral blood samples (2 ml) were obtained from all patients and healthy controls.

All protocols were approved by the Ethics Committee of The Affiliated People's Hospital of Jiangsu University, and the samples were collected at The Affiliated People's Hospital of Jiangsu University after obtaining the written informed consent of all subjects. All procedures described in this study were compliant with the standard biosecurity and institutional safety procedures (26).

Cell isolation and purification

Human PBMCs were isolated via density-gradient centrifugation over Ficoll-Hypaque solution (Tianjin Haoyang Biological Technology Co., Ltd.) according to the manufacturer's instructions. The samples were centrifuged at 600 × g for 20 min and the centrifugal temperature was maintained at 18-20°C. PBMCs were stored at −80°C until use in lncRNA sequencing and reverse transcription-quantitative (RT-q)PCR. Human PBMCs were cultured in RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS (Gibco; Thermo Fisher Scientific, Inc.) at 37°C in 5% CO2.

Library preparation for RNA and lncRNA sequencing

The lncRNAs were derived from various genome annotation projects, specialized lncRNA databases and lncRNA-related literatures, including 'RefSeq' (https://www.ncbi.nlm.nih.gov/refseq/), 'UCSC_KnownGene' (http://genome.ucsc.edu/cgi-bin/hgGateway), 'Gencode' (https://www.genco-degenes.org/), 'Ensembl' (http://asia.ensembl.org/index.html), 'lncRNAdb' (https://lncipedia.org/), 'lncRNAdisease' (http://www.cuilab.cn/lncrnadisease), 'Genbank' (https://www.uniprot.org/database/DB-0028), 'NRED' (http://jsm-research.imb.uq.edu.au/nred/) and 'RNAdb' (https://rnacentral.org/). Library preparation of the RNA and high-throughput sequencing were performed by Cloud-Seq Pte Ltd. A TruSeq stranded total RNA library prep kit (Illumina, Inc.) was used for pretreatment of the RNA prior, extracted with TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.), to the construction of the sequencing library according to the manufacturer's instructions. Quality control and quantitative analysis of the libraries were performed using a Bio-Analyzer 2100 system (Agilent Technologies, Inc.). The 10 pM libraries were desaturated into single-stranded DNA molecules, captured on Illumina Flowcells (Illumina, Inc.) and amplified into clusters in situ, which were sequenced for 150 cycles on an Illumina Novaseq 6,000 sequencing instrument (Illumina, Inc.) in the two-terminal mode (PE mode) according to the Illumina's instructions.

lncRNA profiling analysis

Paired-end reads were obtained from an Illumina Novaseq 6000 sequencing instrument, and quality control was performed using Q30 (27). High quality reads were obtained using Cut-adapt software (v1.9.3) to remove low quality reads for the analysis of lncRNAs (28). High quality reads were subsequently compared with the human reference genome (UCSC HG19) using Hisat2 software (v2.0.4) (http://ccb.jhu.edu/software/hisat2/index.shtml). The fragments per kilobase of exon per million fragments mapped (FPKM) values of lncRNAs, as the lncRNA expression profiles in PBMCs of patients with HT and healthy volunteers, were calculated using Cuffdiff software (v2.2.1) (29), as described in the GTF gene annotation file (https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.13/). Multiple changes and P-values of statistical indicators between the experimental and control groups were calculated to screen for lncRNAs with significantly different expression levels (fold change ≥2.0 and P-value <0.05).

Gene Ontology (GO) (http://www.geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (https://www.genome.jp/kegg/) analyses of differentially expressed lncRNA-associated genes were performed to predict the potential function of lncRNAs (30,31).

RNA extraction and RT-qPCR

Total RNA was isolated from PBMCs using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) at 4°C for 15 min according to the manufacturer's instructions. RNA concentrations were measured using a Nano-Drop ND-1000 instrument (Thermo Fisher Scientific, Inc.). The quality of RNA was determined via the optical density 260/280 ratio. RNA integrity was measured via 1% agarose gel electrophoresis. RT was performed to synthesize the cDNA with random primers in a ReverTraAca® RT-qPCR kit (Toyobo Life Science), which contained oligo-dT and random primers, according to the manufacturer's instructions. RT-qPCR was performed in triplicate using TB GreenTM Premix Ex Taq II (Takara Bio, Inc.). The thermocycling conditions were as follows: Initial denaturation at 95°C for 3 min, followed by 40 cycles at 95°C for 12 sec, 62°C for 40 sec, and 72°C for 32 sec. The primer sequences are presented in Table SI. The expression level of each gene was normalized to the endogenous expression level of the β-actin transcript and was calculated using 2−∆∆Cq method (32). The data were analyzed using Applied Biosystems 7500 Manager software v2.3 (Thermo Fisher Scientific, Inc.).

Small interfering (si)RNA transfection

siRNA (Guangzhou RiboBio Co., Ltd.) was designed to target the sequence of lncRNA-XLOC_I2_006631. A non-specific scrambled siRNA was used as a negative control (NC). The siRNA sequences were as follows: siRNA1, 5′-TTCTCCATAGAGATTTCGG-3′; siRNA2, 5′-CTTGCATAGCTGAGCAACC-3′; and siRNA3, 5'-GCCTCTCAACAACCTCTTT-3′. The isolated human PBMCs were transfected with the lncRNA-XLOC_ I2_006631 siRNA or NC at 100 nM concentration using the Entranster-R reagent (Engreen Biosystem, Co., Ltd.), according to the manufacturer's instructions in the presence of 0.5 μg/ml functional anti-human CD3 Ab (cat. no. 300310; 20 μg/ml) and 2 μg/ml functional anti-human CD28 Ab (cat. no. E-AB-F 11950; 0.5 mg/ml; Miltenyi Biotec GmbH). After incubation at 37°C in 5% CO2 for 24 h, the PBMCs were used for RT-qPCR.

Statistical analysis

GraphPad Prism version 5 software (GraphPad Software, Inc.) was used to manage and analyze the data. A Student's unpaired t-test was used for comparisons of two groups when variables passed the normal distribution test. Mann-Whitney U test was performed to analyze the differences between the two groups of non-normally distributed data. One-way ANOVA test was used for comparisons of multiple groups. Tukey's test was applied for pair-to-pair comparisons of multiple groups after ANOVA. Correlation between the variables was determined using the Pearson's correlation coefficient. Receiver operating characteristic (ROC) curve and the area under ROC curve (AUC) were generated based on the logistic regression model. Sensitivity was used as the Y-axis to represent true positive rate, while 100%-specificity% was used as the X-axis to represent false positive rate. All experiments were performed three times. P<0.05 was considered to indicate a statistically significant difference.

Results

lncRNAs expression profiles in patients with HT

NGS was performed to determine lncRNAs (GEO ID: GSE156468) in the samples. Hierarchical clustering analysis identified that lncRNA expression patterns in PBMCs were different between the HT and control groups, with red and green representing upregulated and downregulated lncRNAs, respectively (Fig. 1A). Volcano plots (Fig. 1B) and scatter plot (Fig. 1C) were used to identified significantly dysregulated lncRNAs. A 2-fold expression difference was used as a cutoff, it was found that the lncRNA expression profiles in the HT group were distinctly different from those in the control group. A total of 97,286 known lncRNAs were identified in PBMCs of patients with HT and healthy controls, including 69,544 upregulated and 27,742 downregulated lncRNAs. A total of 218 novel lncRNAs, including 94 upregulated and 124 downregulated, were verified as having significantly aberrant expression (fold change ≥2; P<0.05; Fig. 1D). Additionally, the 218 identified lncRNAs were divided into six different categories according to the position of lncRNA in the genome. Intergenic lncRNAs accounted for the largest proportion of 38.99% (85/218), exonic sense-overlapping lncRNAs accounted for 27.06% (59/218), intronic sense-overlapping lncRNAs accounted for 6.88% (15/218), natural antisense lncRNAs accounted for 7.8% (17/218), intronic antisense lncRNAs accounted for 9.17% (20/218) and bidirectional lncRNAs accounted for 10.09% (22/218) (Fig. 1E). These results were used in the initial analysis of the sequencing data and were subsequently validated by expanding the sample size.

Figure 1

lncRNA expression profiles in patients with HT. PBMCs of five patients with HT and five healthy volunteers (controls) were enrolled into next-generation high-throughput sequencing. (A) Hierarchical clustering analysis of differentially expressed lncRNAs. Red and green represented upregulated and downregulated lncRNAs, respectively. (B) Volcano plots and (C) scatter plots were used to distinguish differentially expressed lncRNAs. The red squares of volcano plots represented statistically significant dysregulated lncRNAs. The red dots above the diagonal line in the middle of the scatter plot indicated significantly upregulated lncRNAs, and the green dots below represented significantly downregulated lncRNAs. (D) A total of 97,286 known differentially expressed lncRNAs, including 69,544 upregulated and 27,742 downregulated lncRNAs, were identified. Moreover, 218 novel lncRNAs were identified, 94 lncRNAs were significantly upregulated and 124 lncRNAs were significantly downregulated. (E) Number of upregulated (red) and downregulated (green) lncRNAs according to the categories of formation mode. HT, Hashimoto's thyroiditis; lncRNA, long non-coding RNA.

GO and KEGG pathway analysis of differentially expressed lncRNAs in patients with HT

GO and KEGG analyses were performed to determine the potential biological functions of dysregulated lncRNAs. GO enrichment analysis is composed of biological process (BP), cellular component (CC) and molecular function (MF). A total of 257 GO enrichment items, including 128 upregulated and 129 downregulated lncRNAs, were counted (Fig. S1A). The top 10 GO enrichment terms are presented in Fig. 2. The most prominent GO enrichment terms of upregulated lncRNAs were 'membrane assembly' in BP, 'proton-transporting two-sector ATPase complex, catalytic domain' in CC and 'proton-transporting ATP synthase activity, rotational mechanism' in MF (Fig. 2A). Moreover, the most prominent GO enrichment terms of downregulated lncRNAs were 'embryonic viscerocranium morphogenesis' in BP, 'transcription export complex' in CC (A total of eight terms) and 'acetylcholine binding' in MF (Fig. 2B).

Figure 2

GO analysis of differentially expressed lncRNAs in patients with Hashimoto's thyroiditis. GO enrichment analysis is composed of biological process, cellular component and molecular function. The top 10 GO terms of source genes regulated by (A) upregulated and (B) downregulated lncRNAs, respectively. GO, Gene Ontology; lncRNA, long non-coding RNA.

According to the KEGG classification, 17 signaling pathways were associated with potential target genes of upregulated lncRNAs (Fig. S1B), including 'thyroid hormone signaling pathway' (Fig. S2A), and 19 signaling pathways were associated with downregulated lncRNAs (Fig. S1B), including 'calcium-regulated signaling pathway' (Fig. S2B). The top 10 KEGG pathways of dysregulated lncRNAs were shown in Fig. 3A (upregulated lncRNAs) and Fig. 3B (downregulated lncRNAs). Among the 36 signaling pathways, 'longevity regulating pathway-mammal', 'non-alcoholic fatty liver disease' and 'adipocytokine signaling pathway' (Fig. S3) were involved in the regulation of the NF-κB, PI3K-Akt, TGF-β, MAPK and JAK-STAT signaling pathways, which have been reported to participate in the pathogenesis of HT (33-37).

Figure 3

KEGG pathway analysis of differentially expressed lncRNAs in patients with Hashimoto's thyroiditis. The top 10 KEGG pathways for (A) upregu-lated and (B) downregulated lncRNAs, respectively. lncRNA, long non-coding RNA; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Validation of selected lncRNAs via RT-qPCR

To confirm the NGS data, six differentially expressed lncRNAs, including three upregulated [XLOC_I2_006631 (P=0.0051), AL137655_2 (P=0.0055) and LOC729737 (P=0.0053)] and three down-regulated [LOC100288778 (P=0.0378), EPT1 (P=0.0181) and BC041964 (P=0.0010)], were selected based on the differential expression multiplier (fold change ≥2; P<0.05), the FPKM in ≥1 group (FPKM ≥0.5) and the abnormal uniform expression in each HT sample.

These lncRNAs may be involved in the regulation of potential target genes associated with autoimmune thyroid disease. RT-qPCR was performed to verify the differences in selected lncRNAs in PBMCs from 28 patients with HT and 27 healthy volunteers. The results of the assays of five lncRNAs (LOC729737, XLOC_I2_006631, AL137655_2, BC041964 and LOC100288778) were consistent with the NGS data (Fig. 4A-E). However, only lncRNA LOC729737, XLOC_I2_006631 and BC041964 expression levels had statistically significant differences. The results of EPT1 were inconsistent with the trend present in the sequencing data (Fig. 4F). The fold changes of in the expression of six selected lncRNAs in RT-qPCR and sequencing are presented in Fig. 4G. The results indicated that validation via RT-qPCR and sequencing results were inconsistent with regard to the magnitude of changes and significance.

Figure 4

Validation of selected lncRNAs via RT-qPCR. PBMCs were obtained from 28 patients with HT and 27 healthy controls. The transcript levels of (A) lncRNA-LOC729737, (B) lncRNA-XLOC_I2_006631, (C) lncRNA-AL137655_2, (D) lncRNA-BC041964, (E) lncRNA-LOC100288778 and (F) lncRNA-EPT1 in PBMCs from patients with HT and healthy controls were determined via RT-qPCR. (G) Fold changes of six selected lncRNAs expression between RT-qPCR and sequencing were determined. Horizontal lines show the mean. *P<0.05, ***P<0.001. NS, no significance; lncRNA, long non-coding RNA; RT-qPCR, reverse transcription-quantitative PCR; HT, Hashimoto's thyroiditis.

ROC curve analysis of confirmed lncRNAs

The double-blind method was used in the experiments. ROC curve analysis was used to evaluate the potential diagnostic value of lncRNAs in the peripheral blood of patients with HT. The AUC of lncRNA-XLOC_I2_006631 was 0.8479 (95% CI=0.743-0.953; P<0.0001), and the sensitivity, specificity, likelihood ratio and Jorden index were 88.89%, 75.00%, 3.56 and 0.64, respectively. The AUC of lncRNA-LOC729737 was 0.8294 (95% CI=0.718-0.941; P<0.0001), and the sensitivity, specificity, likelihood ratio and Jorden index were 74.07%, 89.29%, 6.91 and 0.63, respectively. The AUC of lncRNA-BC041964 was 0.6032 (95% CI=0.451-0.756; P=0.1892), and the sensitivity, specificity, likelihood ratio and Jorden index were 37.04%, 89.29%, 3.46 and 0.26, respectively (Fig. 5). These data suggested that lncRNA-XLOC_I2_006631 and lncRNA-LOC729737 could potentially be used to differentiate patients with HT from healthy controls and, thus, may be suitable biomarkers of HT.

Figure 5

ROC curve analysis of confirmed lncRNAs. ROC curve analysis of lncRNA-LOC729737, lncRNA-XLOC_I2_006631 and lncRNA-BC041964 was used to distinguish patients with HT from healthy volunteers. ROC, receiver operating characteristic curve; lncRNA, long non-coding RNA; HT, Hashimoto's thyroiditis; AUC, area under the curve.

Potential prediction of lncRNAs target genes

Prediction of the target genes was performed to evaluate the potential functions of validated lncRNAs in HT and to determine whether they can participate in the regulation of gene expression. lncRNA-LOC729737, lncRNA-XLOC_I2_006631 and lncRNA-BC041964 were selected for these studies due to the consistency of the results between the NGS and validation data. Based on the list of the potential target genes associated with the pathogenesis of HT, the potential target gene of lncRNA-LOC729737 was predicted to be STAT3; the potential target gene of lncRNA-XLOC_I2_006631 was predicted to be MECP2 and the potential target gene of lncRNA-BC041964 was predicted to be IL-21R. To confirm these findings, the transcript expression levels of MECP2, STAT3 and IL-21R were determined via RT-qPCR. The mRNA expression of MECP2 was increased in PMBCs of patients with HT compared with that in the healthy controls (Fig. 6A).

Figure 6

Potential prediction of lncRNAs target genes and the effect of lncRNA-XLOC_I2_006631 on the transcription of MECP2. (A) Relative expression of MECP2 mRNA in PBMCs from patients with HT and healthy controls was detected via RT-qPCR. (B) Correlation between the transcript levels of MECP2 and lncRNA-XLOC_I2_006631 in patients with HT. The correlations between the transcript levels of MECP2 and the serum concentrations of (C) TgAb and (D) TPOAb. (E) Purified PBMCs was transfected with lncRNA-XLOC_I2_006631-specific siRNAs and NC (100 nM) in the presence of functional anti-human CD3 mAb and anti-human CD28 mAb. (F) Transcript levels of lncRNA-XLOC_I2_006631 were detected via RT-qPCR after transfection with lncRNA-XLOC_I2_006631-siRNA1-3 and NC. (G) Transcript level of MECP2 was detected after transfection. Each data point represents an individual subject, horizontal lines represent the mean. *P<0.05, **P<0.01, ***P<0.001. lncRNA, long non-coding RNA; HT, Hashimoto's thyroiditis; RT-qPCR, reverse transcription-quantitative PCR; siRNA, small interfering RNA; NC, negative control; MECP2, Methyl-CpG-binding protein 2; TgAb, anti-thyroglobulin-Ab; TPOAb, anti-thyroperoxidase-Ab.

A moderate positive correlation was observed between the transcript levels of MECP2 and the transcript levels of lncRNA-XLOC_I2_006631 (r=0.5747; P=0.0014; Fig. 6B). Additionally, the transcript levels of STAT3 were decreased in patients with HT (Fig. S4A), and there was no correlation between the transcript levels of STAT3 and the transcript levels of lncRNA-LOC729737 (r=0.3318; P=0.0846; Fig. S4B). The transcript levels of IL-21R were increased in patients with HT compared with controls (Fig. S4C). An inconsistent relationship was detected between downregulated lncRNA-BC041964 expression and upregulated IL-21R expression (Fig. S4D); the expression levels of these two RNA species followed an opposing trend, but demonstrated a positive correlation (r=0.5681; P=0.0016; Fig. S4D). These data suggested that lncRNA-XLOC_I2_006631 was associated with MECP2 expression in patients with HT.

Effect of lncRNA-XLOC_I2_006631 on the transcription of MECP2

Utilizing lncRNA database and genome annotation configuration, a gene, XLOC_I2_006631 (transcript ID: TCONS_I2_00012362), was identified to transcribe a spliced, non-coding RNA transcript of 2,768 in length and was located in chromosome 19:18133183-18144152. According to the position compared with the protein coding gene in the genome, lncRNA-XLOC_I2_006631 belonged to the intergenic lncRNA.

To determine whether lncRNA-XLOC_I2_006631 influences MECP2 transcription in patients with HT, the relationship between MECP2 expression and the levels of autoantibodies was initially analyzed. The results indicated that the elevated levels of MECP2 were moderately positively correlated with the serum levels of TPOAb (r=0.4998; P=0.0068; Fig. 6D). However, there was no correlation between the levels of MECP2 and the levels of TgAb in patients with HT (r=0.2251; P=0.2495; Fig. 6C).

Subsequently, purified human PBMCs were transfected with lncRNA-XLOC_I2_006631-specific siRNA and NC. Manipulation of lncRNA-XLOC_I2_006631-specific siRNA resulted in a decrease in the transcript levels of lncRNA-XLOC_I2_006631 compared with that in NC group (Fig. 6E and F). The expression of MECP2 was significantly declined in PBMCs transfected with the lncRNA-XLOC_ I2_006631-specific siRNA compared with that in NC (Fig. 6G). These results suggested that lncRNA-XLOC_I2_006631 regulated the expression of MECP2 in vitro.

Increased expression of lncRNA-XLOC_I2_006631 corre- lates with clinical severity of the disease

TgAb and TPOAb are important antibodies for clinical diagnosis of HT, and their levels can reflect the severity of the disease (10). Correlations between the expression of lncRNA-XLOC_I2_006631 and the serum levels of TgAb or TPOAb were analyzed. A significant moderate positive correlation between the expression of lncRNA-XLOC_I2_006631 and the serum levels of TPOAb were detected (r=0.4717; P=0.0113; Fig. 7B); however, there was no correlation between the expression of lncRNA-XLOC_I2_006631 and the serum levels of TgAb (r=0.3080; P=0.1108; Fig. 7A). These data indicated that dysregulated lncRNA-XLOC_I2_006631 expression was, to a certain extent, associated with HT.

Figure 7

Increased expression of lncRNA-XLOC_I2_006631 correlates with clinical severity of the disease. The correlation between the transcript levels of lncRNA-XLOC_I2_006631 and the serum levels of (A) TgAb and (B) TPOAb in patients with Hashimoto's thyroiditis. Each data point represents an individual subject. TgAb, anti-thyroglobulin-Ab; TPOAb, anti-thyroperoxidase-Ab; lncRNA, long non-coding RNA.

Discussion

Early symptoms of HT can be non-specific and are easily confused with a number of other thyroid disorders, such as simple goiter (38). Furthermore, the diagnosis of HT requires a combination of clinical symptoms and a variety of laboratory tests, including B-flow ultrasonic imaging and indicators of thyroid function and autoantibodies, but there is no single specific laboratory test that can be currently used to diagnose HT (39). Hence, understanding of the mechanisms of pathogenesis may lead to identification of biomarkers that can predict the course of HT. The present study focused on a class of intracellular regulatory molecules, lncRNAs.

lncRNAs are increasingly recognized as important versatile molecules involved in the regulation of DNA, proteins or RNA (40). The mechanisms of action of lncRNAs include modulation of adjacent protein-coding genes in a cis or trans manner by binding to the promoter or enhancer region of the target DNA based on recognition of specific chromatin features, or by functioning as a bridge to connect the DNA and regulatory elements (41). Additionally, lncRNAs can act as precursor molecules of small RNAs or as 'sponges' for microRNAs (miRs) (42). Numerous studies have examined the roles of lncRNAs in the immune system. For instance, HOTAIR acts as a scaffold to induce polycomb repressive complex 2 to the target gene loci to suppress the transcription of the HOXD gene by catalyzing histone H3 lysine 27 methylation (43,44). Maternally expressed 3 has been identified as a competing endogenous RNA of miR-17 that regulates the level and function of Th17 cells in asthma (45). It has also been shown that NONHSAT079547.2 could promote cell proliferation, as well as control IL-17 expression and secretion via the miR-4716-5p/IL-17 axis in CD4+ T cells (25). lnc-DC regulates dendritic cells differentiation by directly binging to STAT3 and promoting STAT3 phosphorylation (46). However, the functional mechanism of lncRNAs in HT has not been fully elucidated.

Previous studies have reported that the inflammatory cells, including Th1 cells, Th17 cells and follicular helper T cells, are significantly increased in the PBMCs from patients with HT, and the proportion of regulatory T cells (Treg) is decreased in the PBMCs from patients with HT (2,47,48). The present study aimed to investigate the underlying causes of these unbalanced inflammatory cells in the peripheral blood of patients with HT. Hence, human PBMCs from patients with HT and healthy volunteers were selected for lncRNAs sequencing. The present results demonstrated differential expression profiles of lncRNAs in patients with HT and healthy volunteers. A total of 218 significantly differentially expressed lncRNAs, including 94 upregulated and 124 down-regulated lncRNAs, were identified in PBMCs of patients with HT. To verify the accuracy of the NGS data performed in a small set of samples, six lncRNAs (XLOC_I2_006631, AL137655_2, LOC729737, LOC100288778, EPT1 and BC041964) were selected for RT-qPCR. The data indicated that five lncRNAs had a trend similar to the NGS results, suggesting a certain accuracy of the NGS data. The differences of lncRNA-XLOC_I2_006631, lncRNA-LOC729737 and lncRNA-BC041964 between the HT and healthy groups were statistically significant. To further investigate potential biological functions of these lncRNAs, GO and KEGG analyses were subsequently performed to examine the NGS data. The results indicated that 'activated T cell proliferation', 'regulation of activated T cell proliferation', 'negative regulation of T cell proliferation' and 'immunoglobulin binding' were the GO terms that may be involved in the occurrence of HT. Additionally, the 'NF-kB', 'JAK-STAT', 'PI3K-Akt', 'TGF-β' and 'MAPK signaling pathways' of the 36 relevant KEGG pathways are involved in the mechanisms of HT (36,37,49-51). These data suggest the potential biological functions of lncRNAs associated with HT.

Non-coding RNAs usually exert regulatory roles in the organism; hence, the present study investigated the role of validated lncRNAs in patients with HT. Prediction software was used to identify the potential regulatory mRNAs of lncRNAs (LOC729737, XLOC_I2_006631 and BC041964). The predicted mRNAs had to be associated with HT. MECP2 was identified as the potential regulatory gene of lncRNA-XLOC_I2_006631, and STAT3 was identified as the potential target gene of lncRNA-LOC729737. It was also demonstrated that lncRNA-BC041964 may regulate IL-21R.

MECP2 is a member of the MBD family and is a pleio-tropic DNA binding protein that preferentially binds to methylated cytosine-phosphate-guanine (CpG) and regulates the expression of multiple methylation-sensitive genes, including T lymphocytes overexpress CD70 (52), Foxp3 (53), secreted frizzled-related protein 4 (54) and patched 1 (55,56). MECP2 can selectively bind to methylated DNA or interact with the histone deacetylase-containing complexes, leading to two distinct epigenetic repression mechanisms, including DNA methylation and histone deacetylation (57). Moreover, MECP2 can facilitate the transcription of genes by binding to unmethylated CpG DNA or combining with cAMP-response element binding protein 1 in promoters (58). Previous studies have revealed that MECP2 is involved in the pathogenesis of various autoimmune diseases, such as RA, SLE and systemic sclerosis (52,54,59). In RA, MECP2 exerts its role by inhibiting secreted frizzled-related protein 4 expression, a negative regulator of the canonical Wnt pathway (54). Research regarding SLE identified the role of MECP2 in aberrant histone modifi-cations within the T lymphocytes overexpress CD70 promoter in CD4+T cells (52). Other studies have focused on single nucleotide polymorphisms of the MECP2 gene in Sjögren's syndrome and autoimmune thyroid diseases (AITD) (60,61). Therefore, it was hypothesized that MECP2 may be associated with the pathogenesis of HT.

In the present study, MECP2 mRNA expression was increased in PBMCs from patients with HT. Moreover, the transcript levels of MECP2 were positively correlated with the serum levels of TPOAb. There was also a significant positive correlation between the elevated levels of lncRNA-XLOC_ I2_006631 and the transcript levels of MECP2. Finally, transfection with lncRNA-XLOC_I2_006631-specific siRNA resulted in a decline in MECP2 mRNA expression. Thus, the present results suggested that lncRNA-XLOC_I2_006631 regulated the expression of MECP2 in patients with HT. However, the underlying mechanism is yet to be elucidated, and further investigations are required to expand on the current findings.

Possible mechanisms of lncRNA-XLOC_I2_006631 involvement in HT were subsequently investigated in the present study. The involvement of Treg in HT development may be one of the possible mechanisms. Foxp3 is the master transcription factor of Treg cells (62). MECP2 is a central element of the upstream CpG-rich enhancer in the Foxp3 gene that is methylated in naive CD4+T cells, activated CD4+T cells and peripheral Treg cells (63). Importantly, induction of demethylation of this CpG site by a DNA methyltransferase inhibitor activates the upstream Foxp3 enhancer to induce Foxp3 expression (64). Treg cells serve a key role in autoimmune pathogenesis by maintaining immune homeostasis and controlling activation of autoreactive CD4+T effector cells (65,66). Our previous study reported that the proportion of Treg cells and Foxp3 mRNA expression were decreased in PBMCs from patients with HT (48). These data suggest that lncRNA-XLOC_I2_006631 may be involved in the pathogenesis of HT by regulating MECP2 expression. However, additional studies are required.

STAT3 is the key transcription factor of Th17 cells, which was identified as a new subset of CD4+T cells involved in the pathogenesis of HT (64). A previous study revealed that MECP2 was indispensable for the differentiation of Th17-mediated pathologies by reinforcing STAT3 signaling in mice (67). The present study hypothesized that the transcript levels of STAT3 should be consistent in Th17 cells, whose proportion is increased in patients with HT. However, the transcript levels of STAT3 were decreased in patients with HT. A possible explanation for this phenomenon is that STAT3 functions via phosphorylation regardless of its high or low levels in HT.

The IL-21/IL-21R signaling pathway is involved in the development of AITD (68). Interestingly, there was a positive correlation between the low expression of lncRNA-BC041964 and high expression of IL-21R in PBMCs from patients with HT in the current study. A positive correlation suggests that two variables change in the same direction, such as when the expression of IL-21R is downregulated as the expression of lncRNA-BC041964 is downregulated in patients with HT. However, the present data suggested that the expression of IL-21R was upregulated in patients with HT. The correlation trend and expression levels of lncRNA-BC041964 and IL-21R in patients with HT were contradictory. Therefore, it was suggested that lnc-BC041964 may participate in HT process via other target genes, but this requires further investigation.

HT is an autoimmune disease characterized by combined effects of multiple autoantibodies, and it is generally accepted that elevated serum concentrations of TgAb and TPOAb are the most impactful manifestations of HT (69). The present data demonstrated that there was a positive correlation between lncRNA-XLOC_I2_006631 expression and serum concentrations of TPOAb. However, lncRNA-XLOC_I2_006631 expression did not correlate with serum concentrations of TgAb. Detection of TgAb and TPOAb is important in HT; however, the levels of these antibodies are not specific in the diagnosis of thyroid disease types (70). High concentrations of TgAb and TPOAb can also be detected in other thyroid diseases, such as Graves' disease and primary hypothyroidism (70,71). Moreover, thyroid autoantibodies are present in a disease-free population (72). Hence, identification of novel biomarkers for HT is important. Thus, the diagnostic value of lncRNA-XLOC_I2_006631 in HT was investigated in the present study. The findings indicated that the AUC of lncRNA-XLOC_I2_006631 in PBMCs from patients with HT was 0.8479, and the sensitivity, specificity, likelihood ratio, and Jorden index were 88.89%, 75.00%, 3.56, and 0.64, respectively. These data suggested that lncRNA-XLOC_I2_006631 was in part associated with the disease process of HT and may be as a potential biomarker, which in combination with thyroid autoantibodies can improve the diagnosis of HT.

The present study reports an intriguing phenomenon, which provides a foundation for subsequent studies on the mechanism and diagnostic value of lncRNAs in HT. However, there are certain limitations to the current study. The validation cohort of lncRNA-XLOC_I2_006631 expression was small. This study focused on lncRNAs profiles in female patients with HT and did not examine differences in lncRNAs profiles between the w and female participants, potentially missing some lncRNAs signatures of male vs. female participants. Additionally, the differences in lncRNAs expression between various regions and ethnicities were not examined. Larger cohorts of patients with HT are required in the future studies to assess the possibility that lncRNA-XLOC_I2_006631 is a novel biomarker. In addition, the underlying mechanism and specific functions of lncRNA-XLOC_I2_006631 in HT should be examined using cellular and animal model experiments.

In conclusion, the present results demonstrate the differential expression patterns of lncRNAs in patients with HT and identified a significant increase in the expression of lncRNA-XLOC_I2_006631. The relationship between lncRNA-XLOC_I2_006631 and disease severity suggests that lncRNA-XLOC_I2_006631 serves an important role in the development of HT. Moreover, lncRNA-XLOC_I2_006631 may positively regulate MECP2 expression. Overall, the present results indicated that lncRNA-XLOC_I2_006631 may be involved in pathogenesis of HT via regulation of MECP2. However, additional investigations are required to confirm and expand these findings.