The study included 2 HT patients, 2 GD patients, and 1 control donor. All the patients and the control donor agreed to a multisample library and sequencing protocol covering all study procedures. The study protocol was approved by the ethics review board of Yantai Yuhuangding Hospital. Demographic and clinical data are presented in Supplementary Table S1 . According to the study requirement, thyroid lesions were obtained from HT and GD patients, and paranodal tissue was obtained from the control donor (who underwent surgery for suspected malignant thyroid nodules) for sequencing. RNA samples were extracted by Trizol (Thermo Fisher Scientific, Waltham, MA, USA). All sequencing was performed at Biomarker Technologies Corporation (Beijing, China).

Sample preparation and cDNA library construction were performed as described in the 10× Genomics Single Cell 3′v3.1 Kit User Guide. The cDNA product and library concentrations were detected according to the qubit 4.0 fluorescence quantitative instrument, and the insertion fragment size of cDNA library was detected by qseq400 biological analyzer. Finally, the sample library was sequenced using novaseq 6000 instrument of Illumina platform. After recognizing Casava base, the obtained original image file was converted into a sequence file and stored in fastq format. The 10× genomics official software CellRanger was then used to compare and quantify the sequencing data.

For long noncoding RNA (lncRNA) sequencing, at the request of the manufacturer, the sequencing library was generated using the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, Beverly, MA, USA), and the index code was added to the attribute sequence of each sample. The TruSeq SR Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA) was used to cluster the index coding samples on the acBot Cluster Generation System. Following cluster generation, the library preparation was sequenced on Illumina Hiseq platform and double terminal readings were generated. The raw data obtained were stored in fastq format. The transcriptome was then assembled using StringTie (20), based on the readings mapped to the reference genome. The assembled transcripts were annotated using the gffcompare program (21).

For microRNA (miRNA) sequencing, small RNA libraries were constructed and library quality was assessed according to the manufacturer’s instructions. Subsequently, the clustering of the index-coded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v4-cBot-HS (Illumina, San Diego, CA, USA). Following cluster generation, the library preparations were sequenced on an Illumina platform and single-end readings were generated. The raw data obtained were stored in fastq format. Considering Homo_sapiens.GRCh38_release95 as the reference genome, the unannotated readings were aligned with the reference genome using the Bowtie software (22).

The total RNA (1 µg) was prepared for the cDNA library using the cDNA PCR Sequencing Kit (SQK-PCS109) protocol provided by Oxford Nanopore Technologies (ONT). The cDNA library was added to FLO-MIN109 flowcells and run on the PromethION platform. The Guppy software in the MinKNOW2.2 software package was used for base calling, and the fast5 format data was converted to fastq format for storage. Subsequently, the full-length sequence was aligned with the reference genome using minimap2 software (23). Following clustering through the alignment information, the consensus sequence was obtained using pinfish software. The consensus sequences of each sample were merged and aligned with the reference genome through minimap2.

Metabolites were extracted from the samples, and supernatants were isolated for sequencing, as required for the experiment. The LC/MS system for metabolomics analysis was composed of Waters Acquity I-Class PLUS ultra-high performance liquid tandem Waters Xevo G2-XS QT of high-resolution mass spectrometer. The column used was purchased from the Waters Acquity UPLC HSS T3 column (1.8 µm, 2.1 × 100 mm).

The positive ion mode consisted of mobile phase A (0.1% formic acid aqueous solution) and mobile phase B (0.1% formic acid acetonitrile). Similarly, negative ion mode consisted of mobile phase A (0.1% formic acid aqueous solution) and mobile phase B (0.1% formic acid acetonitrile). Moreover, the injection volume was 1 μl.

Subsequently, the raw data collected with MassLynx V4.2 were processed by Progenesis QI for peak extraction, peak alignment, and other data processing operations, based on the Progenesis QI online METLIN database (24) and Biomark’s self-built library for identification, and at the same time, theoretical fragment identification and mass deviation. All were within 100 parts per million (ppm). Further analysis was performed after normalizing the original peak area information with the total peak area.

To explore the dysregulated genes of HT and GD, we performed differential expression analysis with the Limma package (25), and p < 0.05 and |log fold change (logFC)| >0.5 were considered significant.

To further explore the biological processes (BP) and pathways involved in dysregulated genes, differentially expressed genes (DEGs) were subjected to BP analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using the clusterProfiler package in R language (26); p < 0.05 were considered significant.

Moreover, using the c2.cp.kegg.v6.2.symbols.gmt in the MsigDB V6.2 database (27) as the background set, the DEGs were subjected to Gene Enrichment Analysis (GSEA) (28). GSEA was performed using clusterProfiler package; p < 0.05 were considered significant.

With reference to a study by Charoentong et al. (29), immune infiltration analysis was performed using the GSEA method to evaluate the differences in immune infiltration between GD and HT patients and control donor, with special focus on cells with high degree of similarity to AITD, such as T cells, T helper cells, and macrophages.

The IntegrateData function of the Seurat package for R language (30) was used to merge single-cell data, and cell clustering analysis was performed according to the default parameters. The clustering results were uniformly reduced and visualized using PCA and T distributions and Stochastic Nearest-Neighbor Embedding (t-SNE), projected onto a 2-dimensional image, defined as a single-cell atlas. Moreover, the FindMarkers function of the Seurat package was used to identify DEGs for each cluster, and p < 0.05 and | logFC| >0.5 were considered significant. In addition, cell-type annotation was performed according to the known markers ( Supplementary Table S2 ).

Ligand-receptor binding is one of the main forms of intercellular signal transduction. Therefore, we performed the cell communication analysis with the iTALK package in R (31). iTALK works by identifying highly or differentially expressed genes in cell clusters, with subsequent matching and pairing of these genes through a ligand receptor database to find important intercellular communication events.

The DEGs involved in the HT and GD regulatory networks were validated by quantitative polymerase chain reaction (qPCR). Trizol (Thermo Fisher Scientific, Waltham, MA, USA) was used to lyse tissues and extract total RNA. The cDNA was then synthesized according to the instructions of the kit (Tiangen, Beijing, China), and the qPCR reaction system was prepared by the qPCR Kit (Gene Star, Beijing, China), in which 2×q-PCR Mix 10 µl, upstream and downstream primers 0.5 µl each, replenish water to 20 µl. Finally, qPCR was performed using Applied Biosystems Stepone plus real-time quantitative PCR instrument (Thermo Fisher Scientific, Waltham, MA, USA) with the following conditions: predenaturation at 95°C for 1 min, followed by 40 cycles of 95°C for 15 s and 60°C for 15 s. The primers used in this study are shown in Supplementary Table S3 .

In this study, all analyses were performed based on the Bioinforcloud platform (http://www.bioinforcloud.org.cn).

The study flowchart is illustrated in Figure 1 . To identify aberrantly expressed lncRNAs and miRNAs in HT and GD patients, differential expression analysis of ONT sequencing and RNA sequencing data was performed ( Supplementary Table S4 ). On analysis, 6,153 lncRNAs and 741 miRNAs were found to be differentially expressed among HT patients ( Figure 2A ). Similarly, 5,492 lncRNAs and 165 miRNAs were observed to be differentially expressed among GD patients ( Figure 2B ). Moreover, the findings of correlation analysis validated the robustness of both sequencing methods used in this study ( Figures 2C, D ). Finally, we could extract genes with high confidence of dysregulated expression among HT and GD patients for follow-up analysis ( Figures 2E, F ).

To identify the molecular dysregulation mechanisms associated with HT and GD, we extracted high-confidence dysregulated genes expressed in HT and GD for enrichment analysis ( Figure 1 ). HT was found to have significant leukocyte–cell adhesion, T-cell activation and differentiation, cytokine metabolic processes, and immune response-activated signal transduction of BP ( Figure 3A ). While, in GD, leukocyte–cell adhesion, hormone metabolic processes, and cell chemotaxis of BP were found to be significantly enriched ( Figure 3B ). Moreover, the KEGG signaling pathway associated with thyroid disease was found to be enriched in both HT ( Figure 3C ) and GD ( Figure 3D ). These include parathyroid hormone synthesis, secretion of thyroid hormone, and action of thyroid hormone signaling pathways and autoimmune thyroid diseases. Simultaneously, several other immune-inflammatory and metabolic signaling-related pathways were enriched in both HT and GD, including cytokine–cytokine receptor interactions, chemokine signaling pathways, Th1 and Th2 cell differentiation, NF-kappa B signaling pathways, phenylalanine metabolism, tryptophan metabolism, and tyrosine metabolism.

To further determine the role of metabolites in HT and GD, we identified differential metabolites in HT ( Figure 3E ) and GD ( Figure 3F ) relative to controls. The differential metabolites were subjected to metabolic pathway enrichment analysis, which revealed that the metabolic pathways of phenylalanine, tyrosine and tryptophan biosynthesis, phenylalanine metabolism, and steroid biosynthesis were significantly activated in HT [ Figure 3G(a) ], while the metabolic pathways of linoleic acid metabolism and steroid hormone biosynthesis were inhibited [ Figure 3G(b) ]. Moreover, steroid hormone biosynthesis, caffeine metabolism, and arachidonic acid metabolic pathways were activated in GD [ Figure 3G(c) ], while phenylalanine, tyrosine and tryptophan biosynthesis, and linoleic acid metabolic pathways were inhibited [ Figure 3G(d) ].

Given the important role of immune responses in HT and GD, we further evaluated the abundance of immune cells in the samples ( Figure 3H ), with varying degrees of infiltration of CD4⁺ T, CD8⁺ T, Mac, Th1, and Th2 cells.

To further understand the effect of dysregulated gene on disease at the cellular level, the abundance of CD4⁺ T cells in HT, GD, and control samples was explored, which are illustrated in Figure 4A . Moreover, we identified the genes that were co-overexpressed with high confidence dysregulation at the bulk level and in CD4⁺ T cells and explored the biological signatures involved in these co-disordered genes. In HT, these co-disordered genes were significantly involved in the thyroid hormone signaling pathway, T-cell receptor signaling pathway, Th1 and Th2 cell differentiation, Th17 cell differentiation, cytokine–cytokine receptor interactions, and NF-kappa B signaling pathway ( Figure 4B ). While, in GD, parathyroid hormone synthesis, secretion and action, ECM–receptor interactions, purine metabolism, cAMP signaling pathway, NF-kappa B signaling pathway, cytokine–cytokine receptor interactions, and T-cell receptor signaling pathway were significantly enriched ( Figure 4C ). In HT, GSEA demonstrated significant enrichment of cytokine–cytokine receptor interactions and T-cell receptor signaling pathways ( Figure 4D ), while no enrichment was observed in GD. Moreover, some dysregulated metabolic signaling pathways including amino acid biosynthesis and purine metabolism were found to be significantly enriched in both HT ( Supplementary Figure S1A ) and GD ( Supplementary Figure S1B ).

The abundance of CD8⁺ T cells in HT, GD, and control is illustrated in Figure 5A . We identified genes that were co-overexpressed with high confidence dysregulation, both at the bulk level and in CD8⁺ T cells, and evaluated the biological signals involved in these co-disordered genes. In both HT ( Figure 5B ) and GD ( Figure 5C ), these co-disordered genes were substantially involved in antigen processing and presentation, Th1 and Th2 cell differentiation, Th17 cell differentiation, cytokine–cytokine receptor interactions, and NF-kappa B signaling pathway. In particular, synthesis, secretion, and action of parathyroid hormone as well as leukocyte migration pathways across the endothelium were enriched in HT. In GD, purine metabolism and cAMP signaling pathways were substantially enriched. In HT, GSEA revealed significant enrichment of cytokine–cytokine receptor interactions and NF-kappa B signaling pathways ( Figure 5D ), while similar enrichment was not observed in GD. However, purine metabolic signaling pathways were enriched in both HT ( Figure S2A ) and GD ( Figure S2B ).

In particular, we found significant differences in the abundance of Mac cells in HT, GD, and control ( Figure 6A ). By exploring the biological signatures involved in genes with high confidence of dysregulated expression at the bulk level and in Mac cells, in both HT ( Figure 6B ) and GD ( Figure 6C ), we observed significant involvement of these co-disordered genes in synthesis, secretion, and action of parathyroid hormone, antigen processing and presentation, Th1 and Th2 cell differentiation, Th17 cell differentiation, NOD-like receptor signaling pathway, the NF-kappa B signaling pathway, and cytokine–cytokine receptor interactions. While in HT, the thyroid hormone signaling pathways and Toll-like receptor signaling pathways were significantly enriched. Phagosomal, phenylalanine, tyrosine, and tryptophan biosynthetic pathways were found to be significantly enriched in GD. In HT, GSEA revealed significant enrichment of cytokine–cytokine receptor interactions, NF-kappa B signaling pathway, Th1 and Th2 cell differentiation, and Th17 cell differentiation pathways ( Figure 6D ), while similar enrichment was absent in GD. Moreover, the purine and phenylalanine metabolic signaling pathways were enriched in HT ( Supplementary Figure S3A ) and GD ( Supplementary Figure S3B ), respectively.

To determine the mechanism of intercellular regulation of the microenvironment of HT and GD, we analyzed the intercellular communication of CD4⁺ T, CD8⁺ T, and Mac cells. This led us to explore the differentially expressed ligands and receptors in different cell types, thereby allowing us to construct a comprehensive regulatory network that might be the mechanism leading to HT and GD. iTALK analysis revealed the presence of a significant difference in the abundance of ligand–receptor pairs in both HT ( Figure 7A ) and GD ( Figure 7B ). Finally, we mapped the regulatory network of dysregulated gene expression between microenvironment T cells and macrophages in HT and GD. In HT, microenvironmental T cells target macrophages via the IL16-CCR5 axis, FGF10-FGFR1 axis, CXCL13-CXCR3 axis, and CD70-CD27 axis, activating cytokine–cytokine receptor interactions and T-cell receptor signaling pathways ( Figure 7C ), whereas, in GD, microenvironmental T cells target macrophages via the CXCR3-CXCL10 axis, PKM-CD44 axis, and CCL21-CCR7 axis, activating EMC receptor interaction and NF-kappa B signaling ( Figure 7D ). Furthermore, qPCR confirmed that some of the genes involved in the HT ( Figure 7E ) and GD ( Figure 7F ) regulatory networks were differentially expressed compared to the control.