Lymphocytic Thyroiditis Transcriptomic Profiles Support the Role of Checkpoint Pathways and B Cells in Pathogenesis

Abstract

Background:

Autoimmune thyroid diseases are the most common types of autoimmune diseases, but their physiopathology is still relatively unexplored. Genotype-tissue expression (GTEx) is a publicly available repository containing RNAseq data, including profiles from thyroid. Approximately 14.8% of these glands were affected by focal lymphocytic thyroiditis and 6.3% were annotated as Hashimoto. We interrogated these data to improve the characterization of infiltrating cells and to identify new molecular pathways active in autoimmune thyroiditis.

Materials and Methods:

Histological GTEx images of 336 thyroid samples were classified into three categories, that is, non-infiltrated thyroid, small focal infiltrated thyroid, and extensive lymphoid infiltrated thyroid. Differentially expressed genes among these categories were identified and subjected to in silico pathway enrichment analysis accordingly. CIBERSORTx deconvolution was used to characterize infiltrating cells.

Results:

As expected, most of the transcriptional changes were dependent on tissue infiltration. Upregulated genes in tissues include—in addition to lineage-specific B and T cell genes—a broad representation of inhibitory immune checkpoint receptors expressed by B and T lymphocytes. CIBERSORTx analysis identified 22 types of infiltrating cells showed that T cells predominate 3:1 over B cells in glands with small infiltrates, only by 1.7:1 in those with large infiltrates. Follicular helper and memory CD4 T cells were significantly more abundant in glands with large infiltrates (p < 0.0001), but the most prominent finding in these glands was an almost sixfold increase in the number of naive B cells (p < 0.0001). A predominance of M2 macrophages over M1 and M0 macrophages was observed in the three gland categories (p < 0.001).

Conclusions:

Analysis of transcriptomic RNA-seq profiles constitutes a rich source of information for the analysis of autoimmune tissues. High-resolution transcriptomic data analysis of thyroid glands indicates the following: (a) in all infiltrated glands, active autoimmune response coexists with suppressor counteracting mechanisms involving several inhibitory checkpoint receptor pairs, (b) glands with small infiltrates contain an unexpected relatively high proportion of B lymphocytes, and (c) in highly infiltrated glands, there is a distinct transcriptomic signature of active tertiary lymphoid organs. These results support the concept that the autoimmune response is amplified in the thyroid tissue.

Introduction

Autoimmune thyroid diseases (AITDs) are the most common organ-specific autoimmune diseases, with prevalence ranging between 2% and 5% of the population (1,2). In addition to well-established forms of thyroid autoimmunity, focal thyroiditis and subclinical thyroiditis affect up to 15–20% of the population older than 60 years (3,4). Despite extensive analysis of the immune response to thyroid autoantigens (5), the fundamental aspects of the etiology and physiopathology of AITD remain unexplained, including whether intrinsic alterations in thyroid follicular cells (TFCs) play a role in disrupting tolerance (6), and which factors determine the transition from latent focal thyroiditis to fully developed AITDs.

In the past 20 years, several new levels of regulation of the immune response have been elucidated, one of them being immune checkpoint receptors (7 –10). Comparison of the transcriptomic profiles of Graves' disease (GD) and control thyroid glands led us to propose the coexistence of mechanisms responsible for activating and suppressing immune responses in autoimmune glands (11). Therefore, it was not surprising to find a partially active programmed cell death protein 1/programmed death-ligand 1 (PD-1/PD-L1) pathway in GD, suggesting that TFCs play an active role in suppressing thyroid autoimmunity once it has been established (12). However, PD-1/PD-L1 is only one of the several regulatory checkpoints that may contribute to the maintenance of peripheral tolerance, which is why a more comprehensive investigation is required.

Genotype-tissue expression (GTEx) is an extraordinary resource in which RNA-seq data from multiple human tissues derived from over 500 cadaveric organ donors have been made available (13), including >300 thyroid gland samples. Given the high prevalence of autoimmune thyroiditis in the general population (3,4), 70 glands with histological thyroiditis were identified in this data repository. This constituted an opportunity to explore the activity of additional regulatory circuits in autoimmune thyroid glands. We interrogated GTEx data to expand the knowledge on the features of infiltrating cells, to identify new molecular pathways involved in autoimmune thyroiditis pathogenesis, and to validate our previous results in GD glands (11). In addition, the analysis of the composition of leukocyte infiltration using CIBERSORTx (14,15) provides a better understanding of the similarities and differences between incipient and well-established lymphocytic thyroiditis.

Materials and Methods

Samples, histological classification, and demographics

The GTEx database v.6 included histological data corresponding to 336 thyroid glands obtained from cadaveric organ donors. After exclusion of 4 glands containing thyroid carcinomas, the degree of lymphocytic infiltration was assessed in 332 samples using GTEx Histology Viewer utility (Supplementary Table S1); of them, matching RNA-seq raw data from 292 thyroids were downloaded from the GTEx database (Fig. 1A). This study was approved by GTEx and Hospital Vall d'Hebron ethics committees.

FIG. 1.

Classification of samples based on the degree of thyroid infiltration. (A) From left to right, representative thyroid sections for samples corresponding to the NIT, SFIT, and ELIT categories as obtained from the GTEx repository. SFIT samples revealed the presence of small infiltrates and ELIT glands had large areas of lymphocytic infiltration that in some cases contained organized germinal centers. (B) Sex distribution of donors for each of the three categories. Extensive lymphocyte infiltration is significantly more common in women. ELIT, extensive lymphoid infiltrated thyroid; GTEx, genotype-tissue expression; NIT, not infiltrated thyroid; SFIT, small focal infiltrated thyroid.

Differential gene expression analysis

Differential gene expression (DGE) was analyzed using DESeq2 library (v.1.18.1) for R (16) and expressed as log₂-fold change. Variance stabilization was achieved using the variance-stabilizing transformation method. The median of ratios was used to correct for different sequencing depths among the samples (17). The false discovery rate was corrected using the Benjamin-Hochberg procedure.

Pathway enrichment analysis

Biological significance of differentially expressed genes (DEGs) was analyzed using Reactome Pathway Knowledge (18), which provides a list of the most represented signaling pathways for each comparison among the sample categories based on the DGE between them.

Immune infiltration profiling

Composition of thyroid leukocyte infiltrates was inferred using CIBERSORTx deconvolution software (14,15,19). Briefly, for each histologically defined category, raw gene counts were normalized to transcripts per million (TPM) and processed using CIBERSORTx algorithm for RNA-seq data. The default LM22 file was selected as the reference file to identify 22 different leukocyte subsets, which were grouped into 6 major populations (Supplementary Table S2). Differences between groups were analyzed using the Kruskal–Wallis test for nonparametric data with Dunn's test correction using GraphPad Prism^® 9.3. (San Diego, CA, USA). The correlation between the abundance of the 22 leukocyte subpopulations was analyzed using R.

Results

Thyroid infiltration is more frequent in female and older patients

A trained histopathologist examined thyroid images and classified the samples into three categories: not infiltrated thyroid (NIT), small focal infiltrated thyroid (SFIT), and extensive lymphoid infiltrated thyroid (ELIT). Lymphocytic infiltration was not observed in 262 samples (79%), while different grades of infiltration were observed in 70 samples (21%), 49 (15%) of which were classified as SFIT and 21 (6%) as ELIT, most of which contained organized lymphoid follicles with germinal centers that could be considered tertiary lymphoid tissue and be suggestive of Hashimoto thyroiditis (HT). As expected, focal thyroid infiltration was more frequent in females and in donors older than 60 years (Table 1 and Fig. 1B). As the proportion of glands with large lymphocytic infiltrates was close to the estimated prevalence of HT in the United States (20), we can consider that this sample is representative of the U.S. adult population.

Table 1.

Demographic Data of Donor Individuals Included in the Study Classified by the Extension of Lymphocytic Infiltration Visible in Genotype-Tissue Expression Bank of Images Into Not Infiltrated Thyroid, Small Focal Infiltrated Thyroid, and Extensive Lymphoid Infiltrated Thyroid

Sample information	NIT	SFIT	ELIT	p ^*
Sample no., n (%)	262 (78.9)	49 (14.8)	21 (6.3)
Sex, n (%)				< 0.01
Female (n = 122)	89 (71.8)	21 (16.9)	14 (11.3)
Male (n = 212)	173 (83.2)	28 (13.5)	7 (3.4)
Age group, n (%)				< 0.01
<40 (n = 41)	32 (78)	7 (17)	2 (5)
40–60 (n = 179)	147 (82)	10 (12)	4 (7)
>60 (n = 112)	83 (74)	21 (19)	8 (7)

p Values calculated using chi-square test.

ELIT, extensive lymphoid infiltrated thyroid; GTEx, genome-tissue expression; NIT, not infiltrated thyroid; SFIT, small focal infiltrated thyroid.

DGE is highly dependent on tissue infiltration

RNA-seq data were first analyzed by principal component analysis to compare the infiltration categories. Samples corresponding to the NIT and ELIT categories were mapped to the opposite sides of the plot. SFIT samples were distributed across and overlapped with the ELIT and NIT area clusters (Fig. 2A), suggesting that the SFIT category might be a transitional state between the NIT and ELIT categories.

FIG. 2.

DGE in SFIT and ELIT glands. (A) PCA distribution of the samples based on gene expression. Analysis showed two differentiation clusters corresponding to NIT (light grey) and ELIT (black) categories, with SFIT (dark grey) samples overlapping both. (B) Volcano plots showing DEGs in SFIT (top) and ELIT (bottom) glands compared to those in NIT glands. Highlighted are genes with an absolute fold change expression >3. (C) Expression changes in genes with structural and functional roles in the normal thyroid tissue in SFIT and ELIT categories compared to NIT glands.

As expected, the transcriptomic profiles of the ELIT glands were highly enriched in lymphocyte-specific transcripts with a moderate reduction in TFCs' restricted genes. This reduction was more pronounced for genes encoding proteins that are involved in thyroid-specific functions (TG, TPO, DIO1, SLC5A5, and DUOX1) than for genes maintaining the thyroid epithelial cell structure (KRT7, KRT18, CLDN3, OCLN, and GJA19) (Table 2). All these changes were very evident in the comparison of the ELIT glands with the NIT glands, but much less evident when comparing SFIT to NIT (Fig. 2B, C and Table 3). Besides the reduction in genes specific to thyroid epithelial cells, there was an even more marked reduction in levels of calcitonin, gastrin-releasing peptide , somatostatin, and chromogranin as if C cells were more negatively affected by the lymphocytic infiltration than thyroid epithelial cells; however, there is no evidence of calcitonin deficiency in patients with thyroid autoimmune disease (21).

Table 2.

Top 10 Upregulated and Downregulated Genes in Extensive Lymphoid Infiltrated Thyroid Glands Compared to Not Infiltrated Thyroid Control Samples

Gene	Protein coded	Tissue or cell expression	Log₂ fold change	p-value adjusted
Upregulated
SERPINA9	Serpin family A member 9	Germinal center associated B cells	8.838	4.99E-11
TCL1A	TCL1 family AKT co-activator A	B cells and	7.880	2.47E-21
FCAMR	Fc fragment of IgA and IgM receptor	Lymphoid tissue, renal epithelium	7.769	8.64E-14
FCRL4	Fc receptor-like 4	Epithelial memory B cells	7.301	2.36E-16
CD19	CD19 molecule	B cells	7.243	1.43E-49
PAX5	Paired box 5	B cells	7.240	1.61E-28
FCRL2	Fc receptor-like 2	B cells	7.119	5.25E-27
FCRL1	Fc receptor-like 1	B cells	7.101	2.08E-29
CR2	Complement C3d receptor 2	B cells	7.038	9.69E-45
FDCSP	Follicular dendritic cell-secreted protein	Dendritic cells	7.028	2.08E-07
Downregulated
CALCA	Calcitonin-related polypeptide alpha	Thyroid C cells	−6496	1.09E-04
GRP	Gastrin-releasing peptide	Ubiquitously expressed	−6417	6.39E-01
CHGA	Chromogranin A	Parathyroid and adrenal gland, other neuroendocrine cells	−4337	6.15E-03
SST	Somatostatin	Digestive and nervous systems, adrenal glands	−3771	9.72E+00
SCGN	Secretagogin, EF-hand calcium binding protein	Digestive and nervous systems, hypophysis	−3527	2.88E-04
VGF	VGF nerve growth factor inducible	Nervous system, hypophysis	−3418	8.49E-01
VSTM2A	V-set and transmembrane domaincontaining 2A	Nervous system	−3274	2.32E+00
LOC100128317	Uncharacterized LOC100128317	—	−2889	1.77E-03
PTPRN	Protein tyrosine phosphatase receptor type N	Nervous system, hypophysis, adrenal glands	−2844	2.57E+00
CSMD3	CUB and Sushi multiple domains 3	Nervous system, hypophysis, testis	−2827	6.59E-03

Table 3.

Top 10 Upregulated and Downregulated Genes in Small Focal Infiltrated Thyroid Glands Compared to Not Infiltrated Thyroid Control Samples

Gene	Protein coded	Tissue or cell expression	Log₂ fold change	p-Value adjusted
Upregulated
IGLL5	Immunoglobulin lambda-like polypeptide 5	B cells	4.154	4.33E-10
CR2	Complement C3d receptor 2 (CD21)	B cells	3.903	2.88E-34
FCAMR	Fc fragment of IgA and IgM receptor	Lymphoid tissue. renal epithelium	3.848	1.04E-06
CXCL13	C-X-C motif chemokine ligand 13	Secondary lymphoid organs	3.773	2.30E-11
FCRL4	Fc receptor-like 4	Epithelial memory B cells	3.767	1.90E-08
CD79A	CD79a molecule	B cells	3.707	2.20E-21
LINC02422	Long intergenic nonprotein coding RNA 2422	—	3.654	8.76E-08
SERPINA9	Serpin family A member 9	Germinal center associated B cells	3.648	9.01E-01
TNFRSF13C	TNF receptor superfamily member 13C (BAFF-R)	B cells	3.555	6.81E-22
FCRLA	Fc receptor-like A	B cells	3.546	3.65E-12
Downregulated
THEG	Testicular haploid expressed gene protein	Spermatocytes	−2.043	3.88E-03
LINC00051	Long intergenic nonprotein coding RNA 51	—	−1.797	3.87E-03
LESD1	Proteoglycan 3. proeosinophil major basic protein 2 pseudogene	Peripheral blood	−1.071	6.61E-03
USP6	Ubiquitin specific peptidase 6	Ubiquitously expressed	−0.989	7.40E-04
LINC02232	Long intergenic non-protein coding RNA 2232	—	−0.812	9.82E-03
ART5	ADP-ribosyltransferase 5	Testis	−0.742	3.10E-04
HSD11B2	Hydroxysteroid 11-beta dehydrogenase 2	Kidney, digestive system, other tissues	−0.576	5.77E-03
RGS17	Regulator of G protein signaling 17	Ubiquitously expressed	−0.441	2.68E-03
MDH1B	Malate dehydrogenase 1B	Testis	−0.395	9.07E-03
KLHL10	Kelch-like family member 10	Testis	−0.376	3.61E-03

The most DEGs specific to immune cells in SFIT and ELIT glands as opposed to those in NIT glands were mostly B lymphocyte-coding genes. Some genes were highly lineage specific (e.g., CD19, CD79A, BLK, and BTK) and other genes were highly expressed in B cells, but shared with other cell types (CR2 or BAFFR). Although not among the top DEGs, there were many upregulated T cell-specific genes in the infiltrated samples compared to those in the NIT group, many of which exhibited high-fold change values in most ELIT glands.

In silico characterization of thyroid immune infiltration

For each subpopulation identified using CIBERSORTx deconvolution software (Supplementary Table S2), its relative fraction in relation to total leukocyte infiltration as well as an absolute score that measures the overall abundance within each sample were calculated. Unexpectedly, the relative presence of each major population in the three gland categories showed that glands histologically classified as NIT contained a small, but not negligible number of leukocyte transcripts, which did not reflect the distribution of leukocyte populations in the blood, but was similar to those infiltrating the SFIT and ELIT glands. This suggests that leukocytes present in NIT glands were not cells trapped in the capillaries, but a mix of myeloid and lymphoid cells that are probably normal constituents of the thyroid tissue (Fig. 3A).

FIG. 3.

Relative presence of leukocyte subpopulations identified by CIBERSORTx. (A) Mean relative distribution of major leukocyte subpopulations within NIT, SFIT, and ELIT categories. (B) Comparison of the relative presence of major leukocyte populations between the NIT (dark grey), SFIT (black), and ELIT (light grey) categories. Infiltrated glands are characterized by a relative increase in B cells and a decrease in NK cells, monocytes, and macrophages. (C) Correlograms showing the relationship between all leukocyte subpopulations identified by CIBERSORTx for the NIT, SFIT, and ELIT categories. (D) From top to bottom, correlation in ELIT samples between M2 macrophages, memory B cells, and follicular CD4 cells and M1 macrophages, naive B cells, and activated NK cells. *p < 0.05, ***p < 0.001, ****p < 0.0001.

In all categories, lymphocytes were the main infiltrating population. While the proportion of T cells in the total infiltrating cells remained around 50% in the three categories, B cells showed a pronounced increase in SFIT and ELIT compared to those in NIT samples (15.0% ± 9.4% vs. 6.9% ± 4.5%, p < 0.0001 and 31.2% ± 13.0% vs. 6.9% ± 4.5%, p < 0.0001, respectively). By contrast, SFIT and ELIT samples showed a relative decrease in NK cells (6.4% ± 3.6% vs. 9.5% ± 4.4%, p < 0.0001 and 4.1% ± 3.3% vs. 9.5% ± 4.4%, p < 0.0001, respectively) and in monocytes and macrophages (26.8% ± 7.9% vs. 30.5% ± 8.3%, p < 0.05 and 16.9% ± 5.3% vs. 30.5% ± 8.3%, p < 0.0001, respectively; Fig. 3B).

To determine whether NIT, SFIT, and ELIT were quantitatively or qualitatively different, the correlations between the proportions of infiltrating leukocyte subpopulations were analyzed separately for each group of samples (Fig. 3C). In NIT samples, the relative infiltration of different leukocyte subpopulations was only weakly correlated with each other; however, multiple mutual correlations were observed in the ELIT category glands. In this category, a remarkable negative correlation was observed between the relative presence of M2 macrophages and memory B cells (r = −0.72, p < 0.01) and follicular CD4 T cells (r = −0.67, p < 0.01). A negative correlation between M1 macrophages and naive B cells (r = −0.58, p < 0.05) and a positive correlation between M1 macrophages and activated NK cells (r = 0.57, p < 0.05) were also observed (Fig. 3D). From the above, it can be deduced that the differences in the infiltrates in NIT, SFIT, and ELIT are both quantitative and qualitative.

The absolute infiltration score obtained for thyroid samples showed a significant increase of leukocytes in glands of SFIT and ELIT categories compared to that in NIT glands (1.25 ± 0.56 vs. 0.91 ± 0.20, p < 0.0001 and 2.57 ± 0.92 vs. 0.91 ± 0.20, p < 0.0001, respectively), confirming the classification generated by histopathological examination of the GTEx images (Supplementary Fig. 1A, B). Naive B cells, memory B cells, and plasma cells were significantly increased in infiltrated tissues (SFIT and ELIT) compared to those in NIT samples. While the relative proportion of total T lymphocytes showed no difference among glands of different infiltration categories, there was an increase in the absolute number of infiltrating resting CD4 memory cells, follicular helper CD4 cells, and CD8 cells in SFIT and ELIT groups compared with NIT glands. Interestingly, infiltrated glands also showed higher activation of NK cells and increased numbers of M1 and M2 macrophages, but not M0 macrophages. All these changes were especially prominent in most infiltrated glands belonging to the ELIT category (Fig. 4A and Supplementary Table S3).

FIG. 4.

(A) Comparison of the absolute presence of different leukocyte subpopulations identified by CIBERSORTx between the NIT (dark grey), SFIT (black), and ELIT (light grey) categories. There is a significant increase in the counts of most B and T cell subpopulations and activated NK cells in infiltrated glands, especially in ELIT samples. In addition, the counts of M1 and M2, but not those of M0 macrophages are increased in SFIT and ELIT samples. (B) DEGs associated with the Th1, Th2, and Th17 responses in SFIT and ELIT categories compared to those in NIT glands. Th1 is the predominant CD4 response in infiltrated tissues. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Th1-associated genes are increased in infiltrated thyroid tissues

Resting CD4 memory cells were the major subpopulation present in the infiltrates of NIT, SFIT, and ELIT glands. This version of CIBERSORTx gives no information on T cell polarization, such as Th1, Th2, and Th17. To overcome this shortcoming, we investigated the DEGs associated with Th polarization in the GTEx data (22). This analysis indicated that Th1 was the predominant CD4 response associated with lymphocytic thyroiditis, which increased further in the more infiltrated tissues belonging to the ELIT category (Fig. 4B).

Immune regulatory pathways are upregulated in SFIT and ELIT samples

Six of the top 10 pathways were common between the SFIT and ELIT gland categories. These pathways included immunoregulatory interactions, co-stimulation by CD28 family molecules, cellular migration, and B cell activation (Fig. 5A and Supplementary Tables S4 and S5). Since PD-1 signaling was found to be within the 10 most enriched biological pathways in the ELIT category, but not in the SFIT category, we investigated whether the most infiltrated glands would also show increased expression of other inhibitory immune checkpoint receptors.

FIG. 5.

(A) Reactome-enriched biological pathways in SFIT (left) and ELIT (right) samples compared to those in NIT glands. Most represented pathways showed a high overlap between the SFIT and ELIT samples and corresponded to processes involved in B and T cell stimulation and regulation of their response. The size of each central node is proportional to the number of upregulated genes associated with them. (B) Differential expression of immune checkpoint molecule-encoding genes in SFIT and ELIT categories compared to those in NIT samples. Infiltrated glands are characterized by the upregulation of different immune checkpoint molecules associated with B and T cell regulation. Changes observed were of a higher magnitude and involved an increased number of genes in ELIT samples. (C) Differential expression of genes involved in the formation and regulation of ectopic lymphoid tissue in SFIT and ELIT glands.

SFIT samples overexpressed only three inhibitory immune checkpoint coding genes with an absolute fold change greater than two: FCRL3 and FCRL4, which were expressed by B cells and TIGIT, a negative receptor mainly expressed by T cells (Fig. 5B). ELIT samples showed increased expression of a broader variety of immune checkpoint receptor genes, including PD-1, LAG3, CTLA-4, and TIM3, which are inhibitory checkpoint receptors of T lymphocytes. Among genes coding for immune checkpoints, only CD112, CD113 and CD115 showed moderate downregulation in ELIT glands when compared to that in the NIT category (Supplementary Table S6).

B cells and germinal center activity revealed

As the formation of tertiary lymphoid organs (TLOs) has been postulated to play a pathogenic role in autoimmune diseases, including thyroid autoimmune diseases (23), we specifically examined the DEGs involved in the formation and activation of lymphoid follicles. Lymphotoxin-alpha, CXCR4, CXCR5, CXCL13, and CCL19 participate in the formation of lymphoid follicles and were highly upregulated in ELIT glands, but only CXCR5 and CXCL13 were significantly upregulated in SFIT glands. The enzyme AICDA, which is essential for B cell somatic hypermutation and class switch recombination processes specific for activated germinal centers, was markedly increased in ELIT (fivefold) and SFIT glands, although by a lesser extent (1.7-fold) (Fig. 5C). In addition, SERPIN9 and CR2, genes typically expressed by centroblasts, were among the top 10 upregulated genes in SFIT and ELIT samples, as was follicular dendritic cell secreted protein (FDCSP), a gene that is specifically expressed by follicular dendritic cells.

CIBERSORTx analysis revealed a marked expansion of Tfh and memory B cells in ELIT glands, but only a very limited expansion in the SFIT category (Fig. 4B). Notably, FCRL4, which codes for an inhibitory receptor originally described in activated B cells from mucosal lymphoid tissues, is one of the DEGs in SFIT and ELIT categories.

Discussion

We report an in-depth analysis focusing on the changes associated with focal and extensive lymphocytic thyroiditis aided by public RNA-seq repositories such as GTEx. A previous analysis of the same dataset by Cho et al. (24) described transcriptomic changes in thyroid glands related to aging, including the upregulation of immune system genes. However, this was not employed in analyzing autoimmune pathogenesis. Both studies are limited by data available at the GTEx, and therefore, clinical and laboratory details relevant to thyroid pathology could not be analyzed.

We have demonstrated that majority of the most upregulated genes in samples with small infiltrates (SFIT) were specifically related to the immune system. Although this pattern was also observed in samples with extensive lymphoid infiltrates (ELIT), this category was characterized by a higher number of DEGs. Most downregulated genes in SFIT and ELIT tissues belonged to ubiquitous signaling pathways, although the expression of functional and structural genes specific to TFCs was also reduced, especially in ELIT samples. This reduction may be an indication of parenchymal destruction.

Pathway enrichment analysis confirmed that most overrepresented biological pathways in SFIT and ELIT categories were related to immune system activation and regulation. In addition, genes coding for immune checkpoint molecules expressed on B and T lymphocytes were significantly increased, including BTLA, CTLA-4, TIGIT, FCRL3, FCRL4, and PD-1 (25 –27). The expression of immune checkpoints was further enhanced in the ELIT category samples, suggesting that their upregulation may be related to the degree of tissue infiltration. Antigen presentation and interferon signatures genes, which were the focus of our previous analysis of GD by array hybridization (11), are also upregulated in the ELIT category of samples. However, due to the different methodologies in the two studies, it is difficult to make a productive comparison.

Characterization of thyroid infiltration using the CIBERSORTx algorithm led to the identification of changes presumably associated with the progression of autoimmune responses. In agreement with previous reports, thyroid infiltrates were mainly composed of T and B lymphocytes (28,29). Other changes observed include an increase of M2 macrophages and the polarization of T cells toward the Th1 response, as reported previously (11,30). In addition, ELIT samples showed a significant increase in infiltrating naive and memory B cell subpopulations and follicular helper T cells, which reflect the known process of TLO formation in thyroid autoimmunity (31 –33). This is further supported by the upregulation of genes involved in germinal center formation in the more infiltrated glands, including LT-α, CXCR5, CXCL13, and CCL19 (34 –37).

Interestingly, among the top 10 DEGs in SFIT and ELIT, we detected that FCRL4, a low-affinity immunoglobulin A receptor with immunoregulatory activity, expressed by a subpopulation of memory B cells initially found to be restricted to mucosa-associated lymphoid tissue (38 –40). However, recent studies have identified an increase of FCRL4+ B cells in TLOs formed in patients with primary Sjögren syndrome or rheumatoid arthritis, and in peripheral mononuclear cells in GD, but not in HT patients. These data suggest a possible role of this cell subset in the pathogenesis of organ-specific autoimmune diseases (41 –43).

Overall, these results provide a better understanding of transcriptomic changes associated with the initiation and progression of thyroid-specific autoimmune responses. Despite the overlap in DGE observed between SFIT and ELIT samples, CIBERSORTx analysis revealed different compositions of infiltrating leukocyte subpopulations without the bias inherent to procedures that include tissue digestion. These results suggest that immune responses in autoimmunity evolve during the course of the disease toward a more chronic state, which is characterized by a better coordinated interplay between infiltrating leukocytes, the presence of TLOs, and the upregulation of immune checkpoint molecules, indicating a shift toward exhaustion of infiltrating cells (44,45). Along with previous works, these results help to understand better the mechanisms underlying thyroid-specific immune-related adverse events observed in cancer immunotherapy after treatment with anti-immune checkpoint receptors antibodies such as Nivolumab (46,47).

Changes that lead to the transition from subclinical forms of thyroiditis to Hashimoto's thyroiditis could not be clearly identified in this study. However, the data point to the recruitment of Th1 cells and the formation of TLO as critical steps that may give rise to an expansion loop that ultimately leads to clinical autoimmunity (23), being the significance of last process further supported by the therapeutic success of anti-CD20 in a large variety of autoimmune diseases, including Graves' ophthalmopathy (48,49).

As the transcriptional data obtained by RNA-seq are dominated by the infiltration itself, it accounts for most expression changes observed and possibly masks those occurring in TFCs. Additional transcriptional studies using isolated TFCs are needed to assess the role of parenchymal cells in local inflammation, and their contribution to peripheral tolerance mechanisms regulating autoimmune responses.

Footnotes

Authors' Contributions

D.Á.-S. contributed to the design of the study, analyzed the data, and wrote part of the article. A.M.-S. contributed to data analysis. A.G.-B., I.B., O.G., E.C., P.M.-L., A.P., C.Z., and C.I. contributed to the review of the article and the interpretation of results. R.P.-B. contributed to the design of the study, to the interpretation of the results, wrote part of the article, and was responsible for the approval of the final draft.

Acknowledgments

We thank the GTEx Project and database of Genotypes and Phenotypes (dbGaP) for granting access to the transcriptomic and clinical information of thyroid tissue donors used for this work. The GTEx Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by the National Cancer Institute (NCI), the National Human Genome Research Institute (NHGRI), the National Heart, Lung, and Blood Institute (NHLBI), the National Institute on Drug Abuse (NIDA), the National Institute of Mental Health (NIMH), and the National Institute of Neurological Disorders and Stroke (NINDS). The data used for the analyses described in this article were obtained from dbGaP accession number phs000424.vN.pN. We also thank Ricardo Gonzalo and the Bioinformatics and Statistics Unit (UEB) from the Vall d'Hebron Institute of Research (VHIR) for their contribution to data analysis.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

This work was supported by grants PI17/00324 from the Instituto de Salud Carlos III and was co-financed by the European Regional Development Fund (ERDF). The GTEx Project was supported by the Commond Fund of the Office of the Director of the National Institutes of Health (NIH), and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS ().

Supplementary Material

Supplementary Figure S1

Supplementary Table S1

Supplementary Table S2

Supplementary Table S3

Supplementary Table S4

Supplementary Table S5

Supplementary Table S6

References

McLeod

DSA

, Cooper

. 2012. The incidence and prevalence of thyroid autoimmunity. Endocrine, 42:252–265.

Madariaga

, Santos Palacios

, Guillén-Grima

, Galofré

. 2014. The incidence and prevalence of thyroid dysfunction in Europe: a meta-analysis. J Clin Endocrinol Metab, 99:923–931.

Kurashima

, Hirokawa

. 1985. Focal lymphocytic infiltration in thyroids of elderly people. Histopathological and immunohistochemical studies. Surv Synth Pathol Res, 4:457–466.

Williams

, Doniach

. 1962. The post-mortem incidence of focal thyroiditis. J Pathol Bacteriol, 83:255–264.

Fröhlich

, Wahl

. 2017. Thyroid autoimmunity: role of anti-thyroid antibodies in thyroid and extra-thyroidal diseases. Front Immunol, 8:51.

Bottazzo

, Hanafusa

, Pujol-Borrell

, Feldmann

. 1983. Role of aberrant Hla-Dr expression and antigen presentation in induction of endocrine autoimmunity. Lancet, 322:1115–1119.

Chen

, Flies

. 2013. Molecular mechanisms of T cell co-stimulation and co-inhibition. Nat Rev Immunol, 13:227–242.

Nüssing

, Trapani

, Parish

. 2020. Revisiting T cell tolerance as a checkpoint target for cancer immunotherapy. Front Immunol, 11:1–16.

Iwai

, Ishida

, Tanaka

, Okazaki

, Honjo

, Minato

. 2002. Involvement of PD-L1 on tumor cells in the escape from host immune system and tumor immunotherapy by PD-L1 blockade. Proc Natl Acad Sci U S A, 99:12293–12297.

10.

Ahmadzadeh

, Johnson

, Heemskerk

, Wunderlich

, Dudley

, White

, Rosenberg

. 2009. Tumor antigen-specific CD8 T cells infiltrating the tumor express high levels of PD-1 and are functionally impaired. Blood, 114:1537–1544.

11.

Ruiz-Riol

, Armengol Barnils

MDP

, Colobran Oriol

, Sánchez Pla

, Borràs Serres

, Lucas-Martin

, Martínez Cáceres

, Pujol-Borrell

. 2011. Analysis of the cumulative changes in Graves' disease thyroid glands points to IFN signature, plasmacytoid DCs and alternatively activated macrophages as chronicity determining factors. J Autoimmun, 36:189–200.

12.

Álvarez-Sierra

, Marín-Sánchez

, Ruiz-Blázquez

, de Jesús-Gil

, Iglesias-Felip

, Gonzalez

, Casteras

, Ferrer-Costa

, Nuciforo

, Colobran

, Pujol-Borrell

. 2019. Analysis of the PD-1/PD-L1 axis in human autoimmune thyroid disease: insights into pathogenesis and clues to immunotherapy associated thyroid autoimmunity. J Autoimmun, 103:102285.

13.

GTEx Consortium. 2013. The Genotype-Tissue Expression (GTEx) project. Nat Genet, 45:580–585.

14.

Newman

, Liu

, Green

, Gentles

, Feng

, Xu

, Hoang

, Diehn

, Alizadeh

. 2015. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods, 12:453–457.

15.

Newman

, Steen

, Liu

, Gentles

, Chaudhuri

, Scherer

, Khodadoust

, Esfahni

, Luca

, Steiner

, Diehn

, Alizadeh

. 2019. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol, 37:773–782.

16.

Love

, Huber

, Anders

. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15:1–21.

17.

Anders

, Huber

. 2010. Differential expression analysis for sequence count data. Genome Biol, 11:R106.

18.

Fabregat

, Sidiropoulos

, Garapati

, Gillespie

, Hausmann

, Haw

, Jassal

, Jupe

, Korninger

, McKay

, Matthews

, May

, Milacic

, Rothfels

, Shamovsky

, Webber

, Weiser

, Williams

, Wu

, Stein

, Hermjakob

, D'Eustachio

. 2016. The reactome pathway knowledgebase. Nucleic Acids Res, 44:D481–D487.

19.

Chen

, Khodadoust

, Liu

, Newman

, Alizadeh

. 2018. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol, 1711:243–259.

20.

Staii

, Mirocha

, Todorova-Koteva

, Glinberg

, Jaume

. 2010. Hashimoto thyroiditis is more frequent than expected when diagnosed by cytology which uncovers a pre-clinical state. Thyroid Res, 3:11.

21.

Grani

, Nesca

, del Sordo

, Calvanese

, Carbotta

, Bianchini

, Fumarola

. 2012. Interpretation of serum calcitonin in patients with chronic autoimmune thyroiditis. Endocr Relat Cancer, 19:345–349.

22.

Radens

, Blake

, Jewell

, Barash

, Lynch

. 2020. Meta-analysis of transcriptomic variation in T-cell populations reveals both variable and consistent signatures of gene expression and splicing. RNA, 26:1320–1333.

23.

Aloisi

, Pujol-Borrell

. 2006. Lymphoid neogenesis in chronic inflammatory diseases. Nat Rev Immunol, 6:205–217.

24.

Cho

, Yoo

, Song

, Kim

, Lee

, Shong

, Park

, Seo

. 2018. Transcriptome network analysis reveals aging-related mitochondrial and proteasomal dysfunction and immune activation in human thyroid. Thyroid, 28:656–666.

25.

Wherry

. 2011. T cell exhaustion. Nat Immunol, 12:492–499.

26.

Murphy

, Nelson

, Šedý

. 2006. Balancing co-stimulation and inhibition with BTLA and HVEM. Nat Rev Immunol, 6:671–681.

27.

, Ehrhardt

, Davis

, Hsu

, Leu

, Ehrhardt

, Cooper

. 2003. The inhibitory potential of Fc receptor homolog 4 on memory B cells. Proc Natl Acad Sci U S A, 100:13489–13494

28.

Pearce

, Farwell

, Braverman

. 2003. Thyroiditis. N Engl J Med, 348:2646–2655.

29.

Margolick

, Hsu

, Volkman

, Burman

, Fauci

. 1984. Immunohistochemical characterization of intrathyroid lymphocytes in Graves' disease. Interstitial and intraepithelial populations. Am J Med, 76:815–821.

30.

Nanba

, Watanabe

, Inoue

, Iwatani

. 2009. Increases of the Th1/Th2 cell ratio in severe Hashimoto's disease and in the proportion of Th17 cells in intractable Graves' disease. Thyroid, 19:495–501.

31.

Weetman

. 2004. Cellular immune responses in autoimmune thyroid disease. Clin Endocrinology, 61:405–413.

32.

Armengol

, Juan

, Lucas-Martín

, Fernández-Figueras

, Jaraquemada

, Gallart

, Pujol-Borrell

. 2001. Thyroid autoimmune disease: demonstration of thyroid antigen-specific B cells and recombination-activating gene expression in chemokine-containing active intrathyroidal germinal centers. Am J Pathol, 159:861–873.

33.

Mohr

, Trésallet

, Monot

, Bauvois

, Abiven

, Atil

, Claër

, Malhotra

, Mayer

, Balderas

, Vaarala

, Deniziaut

, Brocheriou

, Buffet

, Leenhardt

, Gorochov

, Miyara

. 2020. Tissue infiltrating LTi—like Group 3 innate lymphoid cells and T follicular helper cells in Graves' and Hashimoto's thyroiditis. Front Immunol, 11:1–12.

34.

Tang

, Zhu

, Qiao

, Fu

. 2017. Lymphotoxin signalling in tertiary lymphoid structures and immunotherapy. Cell Mol Immunol, 14:809–818.

35.

Neyt

, Perros

, GeurtsvanKessel

, Hammad

, Lambrecht

. 2012. Tertiary lymphoid organs in infection and autoimmunity. Trends Immunol, 33:297–305.

36.

Furtado

, Marinkovic

, Martin

, Garin

, Hoch

, Hubner

, Chen

, Genden

, Skobe

, Lira

. 2007. Lymphotoxin β receptor signaling is required for inflammatory lymphangiogenesis in the thyroid. Proc Natl Acad Sci U S A, 104:5026–5031.

37.

B xuan

, Zhao

L dan

, Zhang × 2019. CXCR4 and CXCR5 orchestrate dynamic germinal center reactions and may contribute to the pathogenesis of systemic lupus erythematosus. Cell Mol Immunol, 16:724–726.

38.

Davis

, Wang

, Kubagawa

, Cooper

. 2001. Identification of a family of Fc receptor homologs with preferential B cell expression. Proc Natl Acad Sci U S A, 98:9772–9777.

39.

Liu

, Goroshko

, Leung

LYT

, Dong

, Khan

, Campisi

, Propst

, Wolter

, Grunebaum

, Ehrhardt

GRA

. 2020. FCRL4 Is an Fc receptor for systemic IgA, but not mucosal secretory IgA. J Immunol, 205:533–538.

40.

Falini

, Tiacci

, Pucciarini

, Bigerna

, Kurth

, Hatzivassiliou

, Droetto

, Verducci-Galleti

, Gambacorta

, Orazi

, Pasqualucci

, Miller

, Kuppers

, Dalla-Favera

, Cattoretti

. 2003. Expression of the IRTA1 receptor identifies intraepithelial and subepithelial marginal zone B cells of the mucosa-associated lymphoid tissue (MALT). Blood, 102:3684–3692.

41.

Verstappen

, Ice

, Bootsma

, Pringle

, Haacke

, de Lange

, van der Vries

, Hickey

, Vissink

, Spijkervet

FKL

, Lessard

, Kroese

FGM

. 2020. Gene expression profiling of epithelium-associated FcRL4+ B cells in primary Sjögren's syndrome reveals a pathogenic signature. J Autoimmun, 109:102439.

42.

Spengler

, Amara

, Clay

, Yeo

, Ramsköld

, Spengler

, Sippl

, Cameron

, Israelsson, Titcombe

, Grönwall

, Sahbudin

, Filer

, Raza

, Malmström

, Scheel-Toellner

. 2017. B cells expressing the IgA receptor FcRL4 participate in the autoimmune response in patients with rheumatoid arthritis. J Autoimmun, 81:34–43.

43.

Rostamzadeh

, Dabbaghmanesh

, Shabani

, Hosseini

, Amirghofran

. 2015. Expression profile of human Fc receptor-like 1, 2, and 4 molecules in peripheral blood mononuclear cells of patients with Hashimoto's thyroiditis and Graves' disease. Horm Metab Res, 47:693–698.

44.

McKinney

, Smith KGC 2016 T cell exhaustion: understanding the interface of chronic viral and autoinflammatory diseases. Immunol Cell Biol. 1–22.

45.

Sharpe

, Pauken

. 2018. The diverse functions of the PD1 inhibitory pathway. Nat Rev Immunol, 18:153–167.

46.

Byun

, Wolchok

, Rosenberg

, Girotra

. 2017. Cancer immunotherapy—immune checkpoint blockade and associated endocrinopathies. Nat Rev Endocrinol, 13:195–207.

47.

Chang

, Barroso-Sousa

, Tolaney

, Hodi

, Kaiser

, Min

. 2018. Endocrine toxicity of cancer immunotherapy targeting immune checkpoints. Endocr Rev, 40:17–65.

48.

GRAID investigators. 2011. Safety and clinical outcomes of rituximab therapy in patients with different autoimmune diseases: experience from a national registry (GRAID). Arthritis Res Ther, 13:R75.

49.

Chen

, Chen

, Sun

. 2000. Intravenous rituximab therapy for active Graves' ophthalmopathy: a meta-analysis. Hormones, 20:279–286.