Quality control and first analysisWe profiled and sequenced spatial gene expression from eight thyroid tissue samples (three HT: HT1, HT2, HT3; three GD: GD1, GD2, GD3, and two controls: C1 and C2). We sequenced each sample with a mean of 84602 ± 38060.11 reads under spots, with a median of 1197 (Q1:740.750-Q3:1792.625) genes per spot. Number of counts, number of features, and percentage of ribosomal and mitochondrial genes were similar between samples (Supplementary Fig. 2a). 407 spots were removed from the HT3 sample due to artificial scratches that correlated with low quality measures (Supplementary Fig. 2b).Single cells from HT samples were analyzed to identify the cell populations reported by ref. 27. After the removal of doublets, we analyzed 27,546 cells. Then, we performed a sample integration with harmony to remove batch effects. At 0.3 resolution, we found 11 major clusters annotated based on their top markers (Supplementary Fig. 3a, b). Correlation with published profiles is represented in Supplementary Fig. 3c.Analysis of histology-based annotation profiling of ST dataTo determine purer cell-type signatures within the tissue, we first annotated each spot based on its histological pattern in the hematoxylin and eosin image (Supplementary Fig. 1). The histologic annotation identified eight regions according to the predominant cell tissue: TFCs, thyrocytes from nodule (Nodule TFCs), TILs, germinal centers (GC), connective tissue (CT), vessels (V), colloid-macrophages (Col) and Borderline (BL). Borderline margins were assigned in those areas between TFCs and TILs (Fig. 1b). The results, as expected, showed a diminished representation of thyrocyte spots and a high presence of immune infiltration and germinal centers in HT samples (Fig. 1b).When we evaluated the performance of histological classification and its reliability using the molecular data, we first integrated the eight samples using Harmony28. UMAP plot showed a correct merge of the spots and a divergence based on the predominant cell population (Fig. 1c). Results from testing the distribution of the samples and conditions in the UMAP suggested an optimal integration, in which TILs and GC were exclusively present in AITD spots and then separated from TFCs and stroma, as expected (Fig. 1d, e). Next, each group was profiled by different canonical markers such as: TG, TPO, and TSHR for TFCs and Nodule (TFCs); CD3D, T Cell Receptor Alpha Constant (TRAC) for TILs; Follicular Dendritic Cell Secreted Protein (FDCSP) and C-X-C Motif Chemokine Ligand 13 (CXCL13) for GC; CD34 and vimentin (VIM) for CT; Von Willebrand Factor (VWF) and Myosin Heavy Chain 11 (MYH11) for V and finally lactotransferrin (LTF) and CD68 for Col. These molecular markers confirmed the correct classification of the histological pattern (Fig. 1f). Supplementary Data 1 lists the top marker genes for each of the seven groups already described.To confirm the unbiased groups for the histology-based approach, we also performed a clustering analysis of the integrated object using the eight samples at 0.1, 0.3, and 0.5 resolution (Supplementary Fig. 4). Regarding the similarities between both approaches, we opted for histological annotation instead of unsupervised analysis to prevent mixed spots and borderline areas.To further increase our knowledge of the stromal microenvironment in AITD pathogenesis, we specifically analyzed four of the major compartments comprising TFCs, connective tissue, vessels, and TILs separately, to depict molecular differences between conditions.Molecular profiling of TFCs revealed altered subpopulations of TFCs in AITDFirst, we focused on TFC state using the Thyrocyte Differentiation Score (TDS). This signature quantifies their differentiation grade by averaging the expression of TG, TPO, Solute Carrier Family 26 Member 4 (SLC26A4), Iodothyronine Deiodinase 2 (DIO2), TSHR, Paired Box (PAX8), Dual Oxidase 1 (DUOX1), Dual Oxidase 2 (DUOX2), NK2 homeobox 1 (NKX2-1), GLIS Family Zinc Finger 3 (GLIS3), Forkhead Box E1 (FOXE1), Trefoil Factor (TFF3) and Four-And-A-Half LIM Domains 1 (FHL1) genes29,30 through AddModuleScore function (Supplementary Data 2). When we statistically analyzed differences in TDS between the three conditions, we found that HT TFCs (median = −0.009 [−0.154, 0.142]) were poorly differentiated compared to control and GD samples (median = 0.785 [0.576, 0.975] and 0.755 [0.539, 0.958], respectively, p adjusted < 0.001). Interestingly, controls had slightly higher TDS than GD samples (p adjusted = 0.006) (Fig. 2a). Spatial visualization of the TDS revealed a differential distribution, with lower values observed in TFCs spots closer to the immune infiltrates (Fig. 2b). Further validation in single cell data from HT thyroids corroborated the existence of two thyrocyte populations (Supplementary Fig. 5a) with high and low TDS (Supplementary Fig. 5b).Fig. 2: Molecular profiling of TFCs identifies altered subpopulations in HT and GD samples.a Violin plots showing the comparison of Thyrocyte Differentiation Score (TDS) of TFCs spots between conditions. HT (n = 3): spots = 1247, median = −0.009 [−0.154, 0.142]; GD (n = 3): spots = 4133, 0.755 [0.539, 0.958]; controls (n = 2): spots = 1884, median = controls: 0.785 [0.576, 0.975]; all adjusted two-sided p values. b Spatial distribution of the TDS in TFCs spots. c UMAP at 0.3 resolution of TFCs spots. d Spatial disposition of TFCs clusters at 0.3 resolution. e UMAP of TFCs spots colored according to conditions: controls, GD, and HT. f Volcano Plot of differential expression analysis between cluster T1 (red) and the rest of the clusters. Blue spots represent statistically significant genes with logarithmic fold changes of more than +0.25 and less than −0.25. g Spatial distribution of the CD74-MIF score in TFCs spots. h Box plot of the CD74 and MIF score in all conditions. HT (n = 3): spots = 1247, median = 2.202 [2.016, 2.394]; GD (n = 3): spots = 4133, median = 1.620 [1.092, 2.036]; controls (n = 2): spots = 1884, median = 0.115 [−0.485, 0.605]; all adjusted two-sided p values. i Representative images of CD74 and MIF immunostaining in serial sections of thyroid tissue from controls, HT, and GD patients. Scale bar: 200 μm. Stainings were confirmed in at least seven biological replicates. j MIF signaling networks: Heatmaps display the condition of each spatial cluster as sender or receiver of MIF ligand and their contribution. Circle plot represents the strengthening of the TFCs clusters as senders of MIF and their receivers. Spatial plot of MIF and CD74-CD44 in controls, HT, and GD samples. k Re-clustering of the AITD TFCs cluster (T1). Top left: UMAP of the re-clustering of the AITD TFCs cluster. Top right: Re-clustering colored by conditions. Below: Barplot of the top KEGG pathways from the differential gene expression between HT (dark red) and GD (dark yellow) in the re-clustering. Dunn test for statistical analysis and adjusted two-sided p values with Holm (a, h) and with Bonferroni (f). TDS thyroid differentiation score, HT Hashimoto´s thyroiditis, GD Graves´ disease, UMAP Uniform Manifold Approximation and Projection for Dimension Reduction, TFCs thyroid follicular cells, CT connective tissue, V Vessels, TILs thyroid infiltrating lymphocytes, SMCs smooth muscle cells, IGs immunoglobulins.To study the molecular differences between conditions, we performed an unsupervised clustering from the integration of the ST spots corresponding to TFCs regions. At a resolution of 0.5 we found six different clusters (Fig. 2c). Cluster T1 (Fig. 2c, d cluster in red) was almost exclusively represented in HT and GD spots (AITD TFCs) (Fig. 2e). Concerning the other clusters, TFC markers such as TG, TFF3, PAX8, and SLC26A7 were upregulated in T0 and T2 (TFC1 and TFC2); mitochondrial genes (i.e., NQO1, NDUFB2 and UQCR11) were upregulated in T3 (High metabolic TFCs, TFC3); some TFC and blood markers were upregulated in T4 (TFC + blood) and CT-related genes such as DCN and FBLN1 were upregulated in T5 (TFC-CT) (Supplementary Data 3).Differential expression (DE) analysis between cluster T1 and the rest of the clusters showed the presence of immunoglobulins, probably due to a cross-reaction between TFCs and plasma cells located in the proximity (Supplementary Data 4, Supplementary Fig. 6a, b). Interestingly, this analysis highlighted the increased expression in TFCs of selenoprotein genes such as Glutathione Peroxidase 4 (GPX4), which are important for the correct removal of reactive oxygen species in the thyroid gland31, as well as different thymosins such as Thymosin Beta-4 (TMSB4X) and Thymosin Beta-10 (TMSB10), which contribute to cytoskeleton regulation and are upregulated in malignant cells30,32(Fig. 2f). Furthermore, ligand-receptor analysis confirmed that the number of interactions involving this AITD-associated TFCs cluster was higher than those interactions involving the other clusters in the TFCs region (Supplementary Fig. 7). When we performed a biological enrichment analysis, we showed an increased presence of autoimmune-related pathways and antigen presentation related molecules in this cluster (Supplementary Data 4). Interestingly, we observed an upregulation of the CD74/Macrophage migration inhibitory factor (MIF) system on TFCs which has been related to injury repair and antigen presentation pathways (Fig. 2f, g, Supplementary Data 4). Indeed, when we statistically analyzed differences of CD74/MIF co-expression within the TFCs spots between conditions, we found significantly higher levels of CD74/MIF in HT (median = 2.202 [2.016, 2.394]) and GD (median = 1.620 [1.092, 2.036]) compared to controls (median = 0.115 [−0.485, 0.605], p adjusted < 0.0001) (Fig. 2h). Spatial plotting focused on TFCs exposed a differential distribution of CD74/MIF co-expression which was higher in areas close to the immune infiltrate, similarly to TDS (Spearman rho test: r = −0.44, p < 0.001), suggesting again two subpopulations differentiated by their proximity to TILs (Fig. 2g). To further validate this presumption, we plotted CD74/MIF on TFCs from HT single cell data confirming that TFCs with low TDS tended to upregulate CD74 and MIF (Supplementary Fig. 5c). Immunostaining of CD74 confirmed the increased expression of this receptor on local TFCs close to infiltrating lymphocytes in HT and GD, reflecting a differential TFCs landscape in these conditions. On the other hand, MIF was expressed in the cytoplasm of both AITD and control TFCs (Fig. 2i and Supplementary Fig. 8).We then attempted to confirm the active involvement of CD74/MIF TFCs in the pathogenesis of AITD by conducting cell-cell interaction analysis within our ST data using CellChat33. The canonical CD74-MIF interaction requires the participation of CD44 as a co-receptor. This analysis identified AITD-associated spots from HT samples as MIF senders, which interacted with immune clusters. However, these TFCs were not MIF receivers (Fig. 2j). Indeed, AITD-associated TFCs showed increased interactions with immune cells related to antigen presentation pathways (HLA or MHC-II), suggesting an alternative role of CD74 beyond the CD74-CD44/MIF axis (Supplementary Fig. 9). On the other hand, AITD associated spots from GD samples were both MIF senders and receivers, with an increased communication probability for the CD74-CD44/MIF interaction compared to controls (Fig. 2j). Moreover, these TFCs did not present interactions with immune cells related to antigen presentation (Supplementary Fig. 9). Interestingly, both conditions exhibited MIF/CD74-CXCR4 interactions, considering AITD-associated spots as MIF senders and immune cell clusters as CD74-CXCR4 receptors (Fig. 2j and Supplementary Fig. 9). This interaction is probably related to T cell chemotaxis, as previously reported34.To identify differential thyroid subpopulations between both diseases, we performed a subsequent re-clustering of this AITD-associated cluster. We obtained seven clusters at 0.3 resolution (Fig. 2k). Notably, we observed a clear separation between HT and GD spots into different clusters: cluster 0 was mostly associated with GD; clusters 2–6 associated with HT, while both conditions overlapped in cluster 1 (Fig. 2k). Regarding the differences observed between HT and GD TFCs, we performed a DE analysis between both conditions. The enrichment analysis of the DE genes in GD thyrocytes showed an increased representation of pathways related to extracellular matrix, cell junctions, cell adhesion and thyroid hormone synthesis (Fig. 2k). On the other hand, in HT there was an increase of reactive oxygen species, electron transport chain (ETC), oxidative phosphorylation and pathways associated with oxidative stress-derived damage (i.e., Huntington, Parkinson, Alzheimer diseases and amyotrophic lateral sclerosis) compared to GD (Fig. 2k). Moreover, we performed a DE followed by an enrichment analysis of TFCs pseudobulk area confirming the pathways associated to AITD and the differences between both diseases already described (Supplementary Fig. 10).In concordance with ST results, these molecular signatures associated with ETC and stress/damaging pathways (Supplementary Data 2) were also confirmed with scRNA-seq data from HT. Interestingly, both were enhanced in thyrocytes characterized by a low differentiation state and an upregulation of CD74/MIF (Supplementary Fig. 5d, e). When we performed the deconvolution of HT spots, “damaged” TFCs were present in higher proportion than”non-damaged” TFCs in all clusters (Supplementary Fig. 6a).Overall, we identified a subpopulation of thyrocytes in AITD that express the CD74/MIF axis close to the immune infiltrate. In HT, this subpopulation exhibited a poor differentiated signature and damaged-associated pathways.Connective tissue profiling unveiled specific fibroblast signatures: GD-associated ADIRF+ myofibroblasts in interfollicular areas and HT inflammation-associated fibroblasts (IAFs) in connective tissue areasTo further address the transcriptional patterns of the connective tissue, spots were first clustered by unsupervised analysis. At 0.3 resolution, we obtained five clusters (Fig. 3a, b; C0-C4). We observed a high expression of thyroglobulin, probably due to an overlap of connective tissue and TFCs (Supplementary Fig. 6a). Regarding the top markers in each cluster, we observed that cluster C0 showed less TG content and was mostly expressed in HT and GD (Fig. 3a–c; Supplementary Fig. 6a); whereas subgroup C1 expressed high levels of hemoglobin genes (HBB, HBA1, HBA2), clusters C2 and C3 expressed thyrocyte-associated markers (i.e., TG, DIO2, PAX8, TPO) and C4 exhibited a high presence of immune cell-related genes (i.e., FDCSP, CD52, CD79A) (Fig. 3a–c, Supplementary Data 5).Fig. 3: Analysis of the connective tissue reveals different fibroblast subpopulations in HT and GD samples.a Spatial distribution of CT clusters at 0.3 resolution. b UMAP of the connective tissue clustering at 0.3 resolution. c Distribution per condition in the UMAP colored by condition. d Hierarchical clustering of the top markers of cluster 0. e Myofibroblasts-like signature: Spatial distribution and density plot across conditions. HT (n = 3): spots = 6403, mean = −0.046 ± 0.770, GD (n = 3): spots = 6583, mean = 0.431 ± 0.844, controls (n = 2): spots = 3592, mean = −0.035 ± 0.725. f IAFs signature: Spatial distribution and density plot across conditions. HT (n = 3): spots = 6403, mean = 0.051 ± 0.429, GD (n = 3): spots = 6583, mean = −0.023 ± 0.362, controls (n = 2): spots = 3592, mean = −0.033 ± 0.280. g Spots deconvolution of HT samples: proportion and location of myofibroblasts and IAFs subpopulations (CXCL12 IAFs and IGFBP6 IAFs). Side bars show the expected cell abundance. Dunn test for statistical analysis and adjusted two-sided p values with Holm (e, f). IAFs inflammatory-associated fibroblasts, HT Hashimoto´s thyroiditis, GD Graves´ disease, UMAP Uniform Manifold Approximation and Projection for Dimension Reduction, ECM extracellular matrix, SMCs smooth muscle cells.To obtain the exclusive molecular signatures of disease-associated connective tissue, we selected cluster C0 to perform a hierarchical clustering of its top 40 marker genes, excluding ribosomal protein related genes (Fig. 3d). We applied functional annotation of each gene based on its molecular expression into different categories: endothelial, fibroblasts associated to extracellular matrix (ECM) organization/inflammatory associated fibroblasts (IAFs), fibroblasts/macrophages, and myofibroblast and smooth muscle cells (SMCs) associated genes (Fig. 3d).Interestingly, we observed two signatures distributing disparately that clustered together across conditions (Fig. 3d): (i) The first subset of genes enhanced in HT and GD included Transgelin (TAGLN), Myosin Light Chain 9 (MYL9) and the Adipogenesis Regulatory Factor (ADIRF), which were attributed to contractile cells, such as myofibroblast-like cells (myofibroblasts or pericytes) or SMCs. (ii) The second group of genes (mostly represented in HT samples) included the Insulin Like Growth Factor Binding Protein 6 (IGFBP6), Tenascin XB (TNXB), Peptidase Inhibitor 16 (PI16), C-X-C Motif Chemokine Ligand 14 (CXCL14) and Fibulin 2 (FBLN2), genes associated to ECM organization, chemotaxis or fibrosis in fibroblasts35,36,37,38,39. This second signature was assigned to IAFs. Spatial representation of the two signatures in the thyroid sections showed a significant increased presence of myofibroblast-like cells within the connective tissue, but also in the thyrocyte area, especially in GD samples (mean = 0.431 ± 0.844) compared to HT (mean = −0.046 ± 0.770) and controls ((mean = −0.035 ± 0.725), p adjusted < 0.001 in both cases), which suggested their presence in the interfollicular space (Fig. 3e, gene signatures in Supplementary Data 2). On the other hand, IAFs were significantly overrepresented in the connective tissue of HT samples (mean = 0.051 ± 0.429) compared to GD (mean = −0.023 ± 0.362) and controls ((mean = −0.033 ± 0.280); p < 0.001 in both cases) (Fig. 3f). To validate the putative presence of myofibroblasts and IAFs and to rule out their association to other cells, we analyzed these signatures in HT single cell data and annotated them based on top markers and enrichment analysis into four cell types: SMCs, myofibroblasts, IAFs related to adhesion and chemoattraction, and IAFs related to ECM organization (Supplementary Fig. 11a, b). When we plotted our ST signatures from CT in HT single cell data, their expression matched with those cell subsets associated with SMCs/myofibroblasts and IAFs respectively (Supplementary Fig. 11c). We plotted main genes individually to further characterize these subpopulations: myofibroblasts were enriched in ACTA2 and ADIRF with low levels of MYH11, and IAFs overexpressed DCN. There were two subpopulations of IAFs, one related to adhesion and chemotaxis with increased CXCL12 levels, and a second subpopulation related to ECM reorganization with upregulated IGFBP6 (Supplementary Fig. 11d). Deconvolution of HT samples located myofibroblasts in TFCs zones close to the immune infiltrate. CXCL12 + IAFs seemed to distribute close to and even inside the immune infiltrate, whereas IGFBP6 + IAFs were restricted to connective tissue regions (Fig. 3g).To confirm the existence and the differential distribution across conditions of myofibroblast and IAFs subpopulations, we performed a multicolor immunofluorescence analysis. In control samples, CD34+ fibroblasts were present in the interstitial space surrounding thyroid follicles and did not express myofibroblast markers (α-SMA, TAGLN, and ADIRF). In contrast, in AITD CD34+ fibroblasts distributed in the interstitial space between thyroid follicles expressed high levels of α-SMA+ and TAGLN+ markers characteristic of myofibroblasts (mean α-SMA IHC score in HT 0.964 ± 0.625 and GD 1.352 ± 1.060 vs 0.266 ± 0.487 in control thyroid tissue, p = 0.004 and 0.001, respectively) (Fig. 4a–c and Supplementary Fig. 12). Interestingly, a special subpopulation of myofibroblast that expressed ADIRF+ was significantly increased in GD tissue compared to both HT and controls (mean ADIRF IHC score in GD 1.552 ± 0.797 vs 0.593 ± 0.647 in HT thyroid tissue and 0.398 ± 0.630 in control thyroid tissue, p < 0.001 in both cases, Fig. 4a, b, d). Accordingly, a strong direct correlation between the IHC score of both ADIRF and a-SMA in GD samples was observed (r = 0.8095; p < 0.0001, Fig. 4e). Interestingly, clinical parameters such as TPOAB showed a direct correlation with both ADIRF and α-SMA scores (r = 0.3243; p = 0.0297 and r = 0.4079; p = 0.0254, respectively. Figure 4e) Overall, these results confirmed that while generic myofibroblasts are increased in both HT and GD, specific ADIRF+ myofibroblasts are enriched in GD and are almost exclusively present in this condition.Fig. 4: Immunostaining of myofibroblasts markers from HT, GD, and control thyroid tissues.a Representative confocal immunofluorescence images, in healthy control, HT and GD tissue samples, of the myofibroblast markers a-SMA (magenta), and ADIRF (red) in combination with a mesenchymal marker, CD34 (green) or a smooth muscle cell marker, MYH11 (magenta). Nuclei are stained with DAPI (blue). Objective: 63X. Scale bar: 50 μm. b Immunohistochemistry of ADIRF and a-SMA in serial tissue sections form controls, HT, and GD patients. Scale 200 μm. c, d Quantification of α-SMA and ADIRF in thyroid tissue from controls, HT, and GD patients. The quantification is described in Methods section. Kruskal-Wallis test and Dunn´s test for statistical analysis and adjusted two-sided p values. e Analysis of correlation between clinical parameters and immunohistochemical scores with Spearman rho test and two-sided p values. The area of the circles shows the absolute value of corresponding correlation coefficients, Data are expressed as the arithmetic mean ± SD. Stainings were confirmed in at least seven biological replicates. HT Hashimoto’s thyroiditis, GD Graves´ disease, TSH thyroid stimulating hormone, TPOAB thyroid peroxidase autoantibodies, TGAB thyroglobulin autoantibodies, TSH-R AB thyroid stimulating hormone receptor autoantibodies.We then validated the presence of IAFs, using DCN as a general marker of this cell subset40. In addition, we studied two different subpopulations, chemoattractant IAFs that co-expressed CXCL12 and ECM reorganizers IAFs that co-expressed IGFBP6. DCN was significantly increased in HT connective tissue compared to controls (p = 0.0256). No difference was observed between HT and GD (p = 0.3147) and GD vs controls (p = 0.4615) (Fig. 5a). Both CXCL12+ and IGFBP6+ IAFs were almost exclusively distributed in HT connective tissue compared to controls and GD (Fig. 5b, c). Interestingly, cell-cell communication analysis detected an increased number of cell traffic interactions between CT clusters (senders) and TILs (receivers) in HT and GD samples (Fig. 5d). In detail, one of the strongest interactions corresponded to the CXCL12-CXCR4 pathway from CT fibroblasts to immune cells from other regions of the tissue, especially in HT (Fig. 5d, e and Supplementary Fig. 13).Fig. 5: Immunostaining of different IAF markers in HT, GD, and control thyroid tissues.a Representative confocal immunofluorescence images, in healthy control, HT and GD tissue samples, of the IAFs marker DCN (red) combined with a mesenchymal/fibroblast marker, CD34 (green) and a myofibroblast marker, α-SMA (magenta). Nuclei are stained with DAPI (blue). Arrows indicate colocalization of CD34 and DCN. Objective: 63X. Scale bar: 50 μm. Right panels, quantification of DCN in thyroid tissue samples. Quantification is described in Methods section.χ2-square test for statistical analysis (p = 0.0338) Fisher´s exact test for individual comparisons and adjusted two-sided p values. b Representative confocal immunofluorescence images, in healthy control, HT, and GD tissue samples, of CXCL12 (red) in combination with a mesenchymal/fibroblast marker, CD34 (green), and a myofibroblast marker, α-SMA (magenta). Nuclei are stained with DAPI (blue). Arrows indicate colocalization of CD34 and CXCL12. Objective: 63X. Scale bar: 50 μm. c Representative confocal immunofluorescence images, in healthy control, HT, and GD tissue samples, of IGFBP6 (green) in combination with DCN (red) and a myofibroblast marker, α-SMA (magenta). Nuclei are stained with DAPI (blue). Arrows indicate colocalization of IGFBP6 and DCN. Objective: 63X. Scale bar: 50 μm. Stainings were confirmed in at least seven biological replicates. d Chemokines sent by CT clusters and their receptors in HT and GD samples. CXCL12-CXCR4 is depicted with dashed lines, and significant interactions between regions are highlighted in yellow. Circle size corresponds to the p value, and probability scores are indicated by color. e Spatial plot of CXCL12 and CXCR4 in controls, HT, and GD samples. AITD autoimmune thyroid diseases, HT Hashimoto´s thyroiditis, GD Graves´ disease, CT connective tissue, TFCs thyroid follicular cells, TILs thyroid infiltrating lymphocytes, V vessels, SMCs smooth muscle cells, IGs immunoglobulins.Altogether, these results validated the transcriptomic data, where we confirmed the increased presence and established the precise spatial distribution of two fibroblasts subpopulations associated with inflammation in HT.Vascular molecular profiling displayed specific markers in AITD related to angiogenesis and increased permeabilityVessel-annotated spots were clustered into five different groups at 0.3 resolution (Fig. 6a). Notably, we did not find substantial differences between conditions (Fig. 6b). DE and pathway enrichment analysis among clusters showed that cluster V0 was enriched in mitochondrial activity-associated genes; V1 exhibited markers related to SMCs (ACTA2, MYH11, TAGLN, ADIRF), V2 and V3 probably presented mixed spots with TFCs, since top genes were related to thyrocytes (i.e., TPO, TG, TFF3), and with TILs (with genes such as TRBC2, CD52 or CD69), respectively. V4 displayed gene sets related to angiogenesis and endothelial development (Supplementary Data 6).Fig. 6: Analysis of ST data from vessel annotated-spots and distribution of PLVAP+ capillaries.a Clustering at 0.3 resolution of the vessels zones from spatial transcriptomic data of AITD and control samples. b Distribution of conditions across the clustering of vessel zones. Volcano plots and remarkable genes of the differential gene expression within vessels areas between (c) HT and controls and (d) GD and controls. e Representation of the expression signature of ETC complex III genes (UQCR11, UQCRQ, UQCRH genes) across conditions within vessels areas. Dunn test: HT (n = 3): spots = 404, median = 0.459 [0.106, 0.736]; GD (n = 3): spots = 211, median = −0.001 [−0.468, 0.517]; controls (n = 2): spots = 254, median = −0.266[−0.567, 0.040]; all p adjusted by Holm test <0.001. f Volcano plot of the general differential expression analysis between HT and GD. g Spatial distribution of PLVAP in control, HT and GD samples (h) Violin plots comparing the quantification of PLVAP+ capillaries between the three groups. T test with Welch´s correction for statistical analysis and adjusted two-sided p values (i) Representative confocal immunofluorescence images, in healthy control, HT, and GD tissue samples, of PLVAP (red) and α-SMA (green). Nuclei are stained with DAPI (blue). Objective: 63X. Scale bar: 50 μm, zoom: 25 μm. Stainings were confirmed in at least seven biological replicates. Adjusted p value with Bonferroni (c, d, f) ST spatial transcriptomics, HT Hashimoto´s thyroiditis, GD Graves’ disease, UMAP Uniform Manifold Approximation and Projection for Dimension Reduction, ETC electron transport chain.To further identify differences within vessels between AITD and control tissues, we performed a DE analysis between conditions using vessel pseudobulks. We showed that genes related to the ETC complex III, such as UQCR11, UQCRQ, and UQCRH, which have been reported to be essential in angiogenesis41, were upregulated in AITD compared to controls (Fig. 6c, d). Indeed, they were more highly expressed in HT (median = 0.459 [0.106, 0.736]), compared to GD (median = −0.001[−0.468, 0.517]) and controls (median = −0.266[−0.567, 0.040]) (p adjusted < 0.001; Fig. 6e). In addition, ACKR1, an endothelial marker related to High Endothelial Venules (HEVs) previously described in HT tissue27, was upregulated in HT vessels compared to control samples, a condition almost exclusively associated with this disease (Fig. 6c). Moreover, EGF Like Domain Multiple 7 (EGFL7), a blood vessel inducer42, and markers previously associated to angiogenesis and permeability in AITD such as endoglin (ENG)25 were also upregulated in HT and GD vs controls (Fig. 6c, d). It is worth highlighting the significant expression in HT and GD vessels of Plasmalemma Vesicle Associated Protein (PLVAP) (Fig. 6c, d), an endothelial cell specific protein associated with diaphragmatic fenestrated vessels43,44,45 that has not been previously reported in an AITD context.General DE analysis between HT and GD confirmed the upregulation of ACKR1 in HT; meanwhile PLVAP was enhanced in GD (Fig. 6f). Concerning spatial distribution in the slides, we observed that ACKR1 was highly expressed in vessels and connective tissue areas of HT samples, as expected, whereas ENG and EGFL7 had a wide expression in the TFCs area of AITD samples, (Supplementary Fig. 14). Interestingly, PLVAP also showed broad expression not only in vessels spots but also in TFCs and CT areas, with particularly high levels observed in GD samples (Fig. 6g). To validate these markers, we analyzed the different molecular profiles of the endothelial cells (ECs) from HT scRNA-seq data. We corroborated previous findings and identified three clusters: ACKR1 + ECs, capillary FLT1 + (Fms Related Receptor Tyrosine Kinase 1) and arterial SEMA3G + (Semaphorin 3 g) (Supplementary Fig. 15a, b). EGFL7 and ENG were exclusively and highly expressed in capillaries and ACKR1 + ECs (Supplementary Fig. 15c). Interestingly, PLVAP expression was upregulated mostly in capillaries (Supplementary Fig. 15c). This correlation might be attributed to fenestrated capillaries that increase their permeability. Moreover, the deconvolution analysis of the HT spatial samples revealed capillaries distributed in TFCs areas as suggested, whereas ACKR1 HEVs distributed outlining the TILs regions (Supplementary Fig. 6c).In agreement with ST data, immunostaining of thyroid sections showed a significantly higher number of PLVAP+ capillaries within the interstitial space surrounding thyroid follicles in AITD compared to control thyroid tissue (mean PLVAP+ capillaries in HT 7.438 ± 2.607 and GD9.714 ± 3.646 vs 4.429 ± 2.472 in controls; p = 0.003 and p = 0.002, respectively, Fig. 6h). Although we did not find significant differences between GD and HT, a trend of increased number of PLVAP+ capillaries in GD was observed (mean PLVAP+ capillaries in HT 7.438 ± 2.607 vs 9.714 ± 3.646 in GD, p = 0.064) (Fig. 6h). Interestingly, no PLVAP expression was found in a-SMA+ arteries from AITD and control tissues, confirming our previous transcriptomic results (Fig. 6i). Moreover, PLVAP score showed a strong direct correlation with both ADIRF and α-SMA scores (r = 0.7136; p = 0.0002 and r = 0.5972; p = 0.0033, respectively, Fig. 4e) and also with TGAB levels (r = 0.5192; p = 0.0473, Fig. 4e).Spatial distribution of thyroid infiltrating immune cellsTo evaluate the spatial distribution of the immune cells, we considered that immune infiltration within the thyroid gland comprises two histologically differentiated zones: TILs areas (HT and GD samples) and GC forming tertiary lymphoid organs (HT samples). After that, we clustered both zones and deconvoluted their spots using single-cell signatures: TCD4, TCD8, NKs cells, ACKR1 + HEVs, myeloid cells, germinal center B cells, proliferative B cells, naïve B cells, and plasmablasts.To address the molecular differences between TILs in AITD and their potential implications, we conducted unsupervised clustering at a 0.3 resolution through sample integration (Supplementary Data 7). Intriguingly, we identified six clusters (Fig. 7a, b): I0, predominantly found in HT samples, histologically encompassed both mantle and infiltration areas. Indeed, I0 exhibited upregulation of genes more associated with FDCs and B cells (FDCSP, CXCL13, VPREB3, CD19), as well as infiltrating T cells (TRAC, CD3D, CD3E, CCL21), suggesting mixed spots (Fig. 7c, e). We aimed to molecularly separate the mantle zone from the immune infiltrate in cluster I0 by performing a re-clustering at a 0.1 resolution. I0 was divided into a cluster closer to the GC, which upregulated FDC-associated genes (Fig. 7d, cluster 1), and another distributed over the infiltrate enriched in CCL19, FOS, and CCL21, among others (Fig. 7d, cluster 0 and Supplementary Data 8). Clusters I1 and I2 comprised both HT and GD spots and occupied borderline zones between dense immune infiltrated areas and gland stroma, characterized by immunoglobulin genes. I2 exhibited high content of thyrocyte markers (i.e., TPO, PAX8, TG). Cluster I3, represented by HT spots, was located in the borderline areas between the immune infiltrate and gland stroma, and showed upregulation of T cell-associated genes (TRAC, TRBC1, CD3D) and B cell markers (MS4A1, CR2). Finally, clusters I4 and I5 incorporated a low quantity of spots: I4 showed upregulation of myeloid markers (CXCL10, SPP1, IFI30) and I5, exclusive of HT samples, increased expression of immunoglobulin genes, myeloid (CD68) and T cell markers (CD8B, GZMA) (Fig. 7e). DE analysis of pseudobulk of TILs regions between HT and GD did not provide meaningful results (Supplementary Data 9).Fig. 7: Analysis of TILs and GC in HT and GD samples. Deconvolution of the spots using immune cells signatures.a UMAP of all the spots of the TILs regions from HT and GD samples. b UMAP colored by the different clusters at 0.3 resolution. c Spatial mapping of TILs clusters. d UMAP and spatial mapping of I0 re-clustering at a resolution of 0.1. e Dot plot of representative markers of immune and stroma cells across clusters of TILs. f Hematoxylin and eosin and GC region detail showing the proportions of deconvoluted immune cells using signatures from scRNAseq analysis. Dotted white lines delineate the mantle (between outer and inner lines) and the GC, including the light and dark zones (inner line). H&E hematoxylin & eosin, HT Hashimoto´s thyroiditis, GD Graves´ disease, UMAP Uniform Manifold Approximation and Projection for Dimension Reduction, GC Germinal Center, TILs thyroid infiltrating lymphocytes, HEVs high endothelial venules.In GC areas, the most upregulated markers corresponded to FDCs (Supplementary Data 1). Spatial location of the top five GC markers revealed a specific concentration of these genes in GC regions, especially FDCSP and TCL1A, suggesting a precise location of FDCs and B cells (Supplementary Fig. 16).To shed light into the distribution of immune cells, we performed spot deconvolution in HT samples. Focusing on the histological areas, GC was enriched in proliferative B cells, GC B cells, and TCD4 (Fig. 7f and Supplementary Fig. 6a). The canonical FDC markers were not detected in the scRNA-seq data, therefore we were not able to map the spatial localization of these cells in the deconvolution analysis. The mantle zone surrounding the GC exhibited an enhanced distribution of naive B cells, TCD4 and, to a lesser extent, TCD8 and NK cells (Fig. 7f). Other dense immune zones also exhibited a higher proportion of T cells (Fig. 7f). Remarkably, myeloid cells (DC and macrophages) were distributed surrounding the infiltrate, whereas ACKR1 + HEVs were located within TILs regions and CT areas. Additionally, plasmablasts disseminated within TFCs areas (Fig. 7f and Supplementary Fig. 6a).Regarding the unavailability of scRNAseq data from the thyroid gland of GD patients to perform spot deconvolution, we plotted different immune gene markers in our ST of GD samples to elucidate their location. Remarkably, our results showed abundant plasma cells not only in infiltrated areas, but also in TFCs regions close to the immune infiltrate and connective tissue. Moreover, we detected a significant expression of myeloid cell markers across the tissue in these samples (Supplementary Fig. 17).
