- 1School of Basic Medicine and Clinical Pharmacy, China Pharmaceutical University, Nanjing, Jiangsu, China
- 2Department of Pharmacy, National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
- 3Department of Medical Oncology, National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
Objective: By integrating whole-genome resequencing (WGR) with longitudinal transcriptomic profiling, this study aimed to unravel the genetic–transcriptional regulatory network underlying thyroid immune-related adverse events (irAEs) in non-small cell lung cancer (NSCLC) patients treated with sintilimab. A key objective was to identify molecular biomarkers with predictive and therapeutic relevance.
Methods: This prospective study included NSCLC patients receiving sintilimab, from whom peripheral blood samples were collected at three time points: baseline, post-first treatment, and post-second treatment. RNA sequencing (RNA-seq) and 30× WGR were performed. Differential gene expression analysis was conducted on the RNA-seq data, followed by longitudinal consistency filtering using the Longitudinal Concordant Gene Intersection (LCGI) algorithm to identify robust differentially expressed genes (DEGs). These DEGs underwent downstream integration with protein–protein interaction (PPI) network analysis and cis-expression quantitative trait loci (cis-eQTL) mapping to pinpoint key genes and regulatory single-nucleotide polymorphisms (SNPs) associated with thyroid irAEs.
Results: The LCGI algorithm identified 13 DEGs exhibiting sustained directional shifts across treatment timepoints. Integration with conventional DEG signatures revealed a functionally cohesive module, with C1QA/B/C, FLT1, TEK, PDGFRB, SPP1, and HLA-DPB1 emerging as central regulators of thyroid irAEs. Cis-eQTL mapping identified 500 SNPs with significant cis-regulatory effects on 153 genes. A “C3 complement-matrix axis” was uncovered as a pivotal node, promoting macrophage polarization toward a pro-inflammatory phenotype. Based on the refined PPI network, we proposed a cascading pathological model in which a self-sustaining feedback loop drove progressive and irreversible thyroid autoimmunity.
Conclusion: This study established a genetic–transcriptional regulatory framework for sintilimab-induced thyroid irAEs and identified a candidate gene set with biomarker potential. Our findings highlighted the central role of complement-driven mechanisms, providing a foundation for precision risk prediction and targeted intervention strategies that preserve antitumor efficacy while mitigating autoimmune toxicity.
Introduction
Immune checkpoint inhibitors (ICIs) have revolutionized the therapeutic landscape of non-small cell lung cancer (NSCLC), offering durable clinical benefit in a subset of patients (1, 2). Sintilimab, a fully human IgG4 monoclonal antibody targeting PD-1, exhibits high-affinity binding and sustained blockade of the PD-1/PD-L1 pathway (3), leading to significantly improved survival outcomes in NSCLC (4). However, thyroid immune-related adverse events (irAEs) have emerged as some of the most prevalent toxicities associated with ICI-based regimens, with incidence rates approaching 50% in combination therapies (5). These endocrine irAEs span a clinical spectrum that includes hypothyroidism, hyperthyroidism, and destructive thyroiditis, often necessitating treatment discontinuation and substantially impairing patient quality of life (6, 7). Despite its proven clinical efficacy, the genetic and transcriptomic underpinnings of sintilimab-induced thyroid irAEs in NSCLC remain largely undefined.
The advancement of gene sequencing technologies has enabled major breakthroughs in elucidating the genetic basis of disease pathogenesis (8, 9). Genome-wide association studies (GWAS) have been widely used to identify genetic loci associated with various traits and diseases, but they generally require large cohorts to achieve statistically robust and reproducible genome-wide associations (10). As a result, for studies involving small sample sizes, such as investigations into adverse drug reactions, GWAS often prove insufficient as a stand-alone analytical strategy. Moreover, even in large-scale GWAS, although numerous trait-associated loci have been successfully identified, the mechanistic interpretation of these variants remains an ongoing challenge: the specific molecular pathways through which these variants influence biological functions and disease phenotypes are often poorly understood (11).
Expression quantitative trait loci (eQTL) studies illuminate the functional consequences of genetic variants by linking them to transcriptional changes, offering critical mechanistic insights into GWAS-identified loci and the etiology of complex diseases (12). The integration of multiple expression modalities, such as total gene expression and allele-specific expression, significantly enhances both the sensitivity and specificity of eQTL detection, particularly in heterogeneous tissues like tumors. Advanced mapping strategies such as pTReCASE, which incorporate tumor purity estimates and distinguish between tumor-derived and normal cellular expression, offer a more refined and accurate framework for detecting tumor-specific regulatory effects (13). These methodological advances not only improve the identification of disease-associated genes and actionable therapeutic targets but also deepen our understanding of the molecular circuitry driving disease progression (14). In the context of immunological research, eQTL analyses hold distinct advantages by capturing cell- and tissue-specific regulation, temporal expression dynamics, and multicellular interactions, thus providing powerful tools to dissect the genetic basis of immune-related disorders and inform novel therapeutic strategies (15).
Previous investigations into irAEs have largely been restricted to static comparisons of gene expression between patient groups following a single ICI administration. However, immune responses are inherently dynamic and temporally regulated. To better capture this complexity, the current study assessed transcriptomic profiles at three sequential time points: pre-treatment, post-first treatment, and post-second treatment. Differential expression and pathway enrichment analyses were performed to systematically characterize the evolving transcriptomic landscape in response to sintilimab administration.
This study leveraged transcriptomic data to perform an unbiased, genome-wide eQTL analysis aimed at uncovering single-nucleotide polymorphisms (SNPs) linked to gene expression changes implicated in thyroid irAEs. By integrating SNP profiles with longitudinal transcriptomic data, we reconstructed the genetic-transcriptional regulatory network underlying thyroid irAEs, thus addressing a key limitation in this field, namely, the constraints of small sample sizes in rare adverse event studies. This multi-omics approach illuminated the mechanistic connections between genetic variation and the emergence of thyroid irAEs, revealing candidate molecular pathways that might drive their development. Our findings also pointed to a set of clinically actionable biomarkers with potential utility for dynamic risk stratification in NSCLC patients undergoing sintilimab treatment. Through this comprehensive framework, we provided new insights into the genetic architecture of immune toxicity, with broader implications for precision medicine in cancer immunotherapy.
Materials and methods
Study design and participants
This study was designed as a prospective, real-world investigation. Eligible participants met the following criteria: (1) hospitalization at the National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences between January 2021 and December 2024; (2) histological or clinical confirmation of NSCLC; (3) completion of at least two cycles of sintilimab administered every 3 weeks; (4) age between 18 and 80 years; and (5) availability of complete and evaluable baseline (BT) data prior to treatment initiation.
Exclusion criteria included: (1) history of thyroid surgery or radiotherapy to the cervical region; (2) preexisting thyroid disease or dysfunction; and (3) severe cardiac, hepatic, renal, or autoimmune disorders.
The study protocol was approved by the Ethics Committee of the National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences (approval number: 2020122509100302). Written informed consent was obtained from all participants before enrollment.
Whole-Genome Resequencing
Peripheral blood samples were collected, and genomic DNA was isolated using the TIANamp Genomic DNA Kit. DNA concentration was quantified with a Qubit® 3.0 Fluorometer (Life Technologies, CA, USA), while integrity and purity were assessed via 1% agarose gel electrophoresis (120 V, 45 min). Library preparation was performed following the TruSeq DNA Sample Preparation Guide (Illumina, 15026486 Rev. C). Library quality was verified by qPCR using the Bio-Rad iQ SYBR Green Kit.
Sequencing was conducted on the Illumina NovaSeq 6000 platform using paired-end (PE) protocols, generating 150-bp PE reads with an average coverage depth of 30×. Clean reads were aligned to the UCSC hg38 human reference genome using BWA. Variant calling was subsequently performed using the Genome Analysis Toolkit (GATK) to extract candidate polymorphic SNP loci.
RNA sequencing
High-quality total RNA was extracted from peripheral blood samples. RNA purity was assessed using a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and RNA integrity was evaluated using the Agilent 2100/5400 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Ribosomal RNA was depleted from total RNA prior to library construction.
Fragmented RNA was used as a template for first-strand cDNA synthesis primed with random hexamers. Second-strand synthesis was performed using DNA Polymerase I, RNase H, and a dNTP mix in which dUTP replaced dTTP. The resulting double-stranded cDNA underwent purification, end repair, adenylation, and adapter ligation. Uracil-containing second strands were selectively digested with USER enzyme (New England Biolabs) to ensure strand specificity. Final libraries were purified using AMPure XP beads (Beckman Coulter). Libraries that passed quality control were sequenced on the Illumina NovaSeq X Plus platform using a 150-bp PE protocol.
Statistical analyses
The overall analytical workflow is illustrated in Figure 1. Differential gene expression analysis was conducted using the DESeq2 package (version 1.34.0), while eQTL mapping was performed using the Matrix eQTL package (version 2.3) in R (version 4.3.3).
Figure 1. Research process for investigating immune-related thyroid dysfunction in NSCLC patients receiving sintilimab treatment.
Differential expression analyses
RNA-seq quantification data were analyzed using the DESeq2 package (version 1.34.0) in R (version 4.3.3). Raw gene counts were imported from a preprocessed expression matrix, accompanied by metadata specifying experimental groups and covariates (sex and age). Genes with low expression, defined as fewer than 10 counts in fewer than k samples, where k corresponds to the sample size of the smallest comparison group, were excluded to reduce noise and retain biologically meaningful signals.
Categorical variables (e.g., sex and group) were encoded as factors, while continuous variables (e.g., age) were standardized via z-score transformation (mean-centered and scaled to unit variance). Differential expression testing was performed using a negative binomial generalized linear model that accounted for sex and age as covariates. Dispersion estimates were moderated using empirical Bayes shrinkage prior to Wald significance testing.
Differentially expressed genes (DEGs) were defined based on dual criteria: an absolute log2 fold change >1 and an adjusted p-value (padj, Benjamini–Hochberg correction) <0.05.
Gene expression can be influenced by a variety of environmental and biological factors, including diet and lifestyle (16), chemical exposure (17), hormonal fluctuations (18, 19), and cellular states such as proliferation and differentiation (20). To mitigate confounding from such variables, we implemented the Longitudinal Concordant Gene Intersection (LCGI) method, a robust analytic strategy designed to reduce false positives while prioritizing sustained biological signals. This approach incorporates longitudinal RNA-seq data from two post-treatment time points, filtering genes by fold-change thresholds (A1T >1, A2T >0.5), p-value significance (<0.05), and directional consistency across time. The final DEG set was defined as the intersection of genes meeting all criteria, thereby enriching for genes with persistent, treatment-associated expression changes.
Gene ontology functional annotation and kyoto encyclopedia of genes and genomes pathway enrichment analysis
To minimize baseline confounding and achieve a comprehensive characterization of pathways associated with thyroid irAEs, we integrated enrichment analyses from two gene sets: (1) concordantly regulated DEGs identified via the LCGI method; and (2) combined DEGs from the post-treatment stages A1T and A2T, excluding those differentially expressed at BT.
GO functional annotation and KEGG pathway enrichment analyses were performed using the DAVID online platform (https://davidbioinformatics.nih.gov/home.jsp). Visualization of enrichment results was achieved through Sankey diagrams generated with Bioladder (https://www.bioladder.cn) and bubble plots created via Wei Sheng Xin (http://www.bioinformatics.com.cn).
Construction of protein-protein interaction networks
DEGs across the three time points were compared and visualized using Venn diagrams. PPI networks were constructed based on these DEGs using the STRING database. The resultant interaction data were imported into Cytoscape software to analyze network topology and facilitate visualization. Key genes within the network were identified using the CytoHubba plugin employing the maximal clique centrality (MCC) algorithm, while functionally relevant gene clusters were detected using the MCODE plugin. Integration of these analyses enabled the identification of significantly connected and functionally important gene modules.
eQTL analysis
Cis-acting genetic associations between SNPs and gene expression levels were identified using a linear model framework implemented in the MatrixEQTL R package (version 2.3). SNP-gene pairs located within a 1-Mb window flanking each gene’s transcription start site were tested for association. All significant cis-eQTLs were defined by a false discovery rate (FDR) threshold of <0.05 and retained for subsequent functional annotation.
Integrative analysis
Genes identified through cis-eQTL mapping were intersected with DEGs derived from transcriptomic analyses. The overlapping gene sets were then analyzed via the STRING online platform, applying a network interaction confidence cutoff set to “confidence” to delineate high-confidence PPIs. Integrating these data with prior knowledge from the literature enabled the inference of putative molecular mechanisms driving thyroid irAEs.
Results
Patient characteristics
A total of 29 patients were enrolled in this study. Within this cohort receiving combined sintilimab and chemotherapy, 21 patients (72.41%) maintained normal thyroid function indices at BT, including T3, T4, FT3, FT4, and TSH, and exhibited persistently normal TSH levels following treatment. In contrast, eight patients (27.59%) with normal BT thyroid function developed abnormal TSH levels post-treatment (Table 1). Importantly, no statistically significant differences were observed between these two groups with respect to sex, age, tumor histological subtype, disease stage, concomitant medication regimens, or BT thyroid function parameters.
Table 1. BT characteristics and clinical parameters comparison between patients with normal thyroid function and thyroid dysfunction after sintilimab combined with chemotherapy treatment.
Differential expression analyses and dynamic expression patterns
Differential expression analysis was performed to identify genes exhibiting significant expression differences between the thyroid irAEs group and controls across three distinct timepoints (BT, A1T, A2T). Volcano plots (Figure 2) illustrate gene expression changes at each stage. Raw p-values were adjusted using the Benjamini–Hochberg procedure to control the FDR, generating padj. In these plots, red dots denote genes meeting the criteria of padj < 0.05 and absolute log2 fold change > 1, while green dots indicate genes with padj < 0.05 and absolute log2 fold change > 1.
Figure 2. Volcano plot highlighting significant gene expression variations across the 3 distinct stages: BT, A1T, A2T.
To assess the impact of sintilimab treatment on gene expression dynamics, we applied the LCGI method, identifying genes with consistent directional changes, either upregulated or downregulated, at both A1T and A2T stages, while rigorously excluding genes exhibiting similar BT expression changes. This approach yielded 13 DEGs, comprising two upregulated and 11 downregulated genes (Figure 3).
To capture progressive transcriptional adaptations, we retained a primary threshold of p < 0.05 and |log2 fold change| > 1 (highlighted as green dots on volcano plots) across all comparisons, while relaxing the fold change cutoff to |log2FC| > 0.5 specifically for the adjacent A2T phase to detect subtler expression shifts. Heatmaps (Figures 4, 5) visualize critical gene expression patterns across samples post-first and second treatment cycles, with hierarchical clustering dendrograms at the top illustrating sample similarity.
GO functional annotation and pathway enrichment analysis
GO enrichment results are visualized in the bubble plot (Figure 6). Key enriched terms, considering both enrichment ratio and p-value, converged on complement activation as a central biological process, implicating neuro-immune crosstalk alongside perturbations in metabolic and transport pathways that collectively contribute to thyroid irAEs.
Subsequent KEGG pathway enrichment analysis, performed via the DAVID database, identified 18 DEGs enriched across multiple pathways. Eleven significantly enriched pathways were visualized using a bubble chart (Figure 7), and a Sankey diagram (Figure 8) was constructed to systematically map the gene-pathway relationships. These enrichment patterns emphasized the functional interplay between the complement and coagulation cascades (hsa04610) and the MAPK signaling pathway (hsa04010).
Construction of PPI networks
DEGs were analyzed for PPIs using the STRING database, and the resulting network was visualized with Cytoscape software. After filtering out disconnected nodes based on network topology, the refined network comprised 29 nodes and 45 edges (Figure 8).
Key regulatory nodes were identified using the CytoHubba plugin with the MCC algorithm, highlighting genes with high connectivity or essential functional roles. Notably, PDGFRB, TEK, FLT1, C1QB, SPP1, C1QC, and C1QA emerged as top candidates (MCC ≥ 10).
Subnetwork detection via MCODE, using a K-core threshold of 2, revealed modules of densely interconnected nodes. The highest-scoring subnetwork included PDGFRB, EFNA1, C1QB, TMEM119, C1QA, TEK, HLA-DPB1, FLT1, and C1QC, connected by 16 edges. This module surpassed the degree cutoff threshold, designating it as the most functionally pivotal cluster within the PPI network.
eQTL and integrative analysis
Cis-eQTL analysis identified 514 significant associations involving 512 unique SNPs and 135 genes (Figure 9). Notably, the gene ADARB2 was found at the intersection of both the eQTL target gene set and DEGs. Three SNPs, rs12763675, rs2209625, and rs10904542, were implicated as putative eQTLs modulating ADARB2 expression, suggesting a potential role in the pathogenesis of thyroid irAEs.
Further integrative analysis using the STRING database revealed 173 distinct eQTLs regulating 37 genes that directly interact with DEGs. These 37 genes, in turn, influenced 39 downstream target genes, highlighting a complex regulatory network potentially contributing to thyroid irAE development.
Discussion
The LCGI approach robustly filters out transient fluctuations by demanding consistent gene expression changes across multiple timepoints, thereby elevating confidence in detecting biologically meaningful and sustained signals. By requiring genes to exhibit concordant directional dynamics, either persistent upregulation or downregulation, across independent temporal measurements, LCGI effectively acts as an intrinsic experimental replicate. This design dramatically reduces the joint false-positive rate (<0.0025, derived from 0.05 × 0.05), markedly outperforming conventional single-timepoint FDR-controlled analyses, a feature especially advantageous in studies constrained by limited sample sizes.
Importantly, given that irAEs reflect progressive disturbances in immune homeostasis, such as sustained T-cell activation and cumulative inflammatory cytokine release, LCGI selectively enriches for genes with durable transcriptional responses intimately linked to irAE pathogenesis, while efficiently excluding ephemeral stress-induced gene expression. This targeted focus on stable molecular signatures not only enhances mechanistic insight but also improves the identification of robust biomarkers predictive of irAE onset.
Among the 13 genes identified via LCGI, several emerged as central hub nodes within the PPI network and have been experimentally validated for their pivotal roles in immune regulation. Notably, SPP1 (secreted phosphoprotein 1, also known as Osteopontin) orchestrates both innate and adaptive immune responses by modulating immune cell activation, proliferation, and differentiation, thereby finely tuning the magnitude of immune reactions (21). Through interactions with integrins (e.g., αvβ3) and the CD44 receptor, SPP1 facilitates immune cell trafficking to sites of inflammation (22, 23). Its chemotactic properties likely promote T-cell migration into healthy tissues, potentially driving the pathogenesis of irAEs.
Functionally, SPP1 modulates the secretion of key pro-inflammatory cytokines, including IL-1β, TNF-α, and IL-6, while concurrently dampening cytotoxic CD8+ T cell activity (23). Prior studies have demonstrated a strong positive correlation between SPP1 expression and the infiltration of diverse immune cell subsets in cancer, implicating it in the establishment of an immunosuppressive tumor microenvironment (24, 25). Moreover, genetic variants of SPP1 have been significantly linked to susceptibility and clinical severity in autoimmune diseases such as systemic lupus erythematosus (SLE) and rheumatoid arthritis (RA), underscoring its potential as a biomarker for irAEs (21).
Within the PPI network, C1QA, C1QB, and C1QC formed a prominent subnetwork identified through CytoHubba and MCODE analyses, implicating their involvement in thyroid irAEs. C1q, a hexameric protein complex composed of these three distinct subunits, plays a central role in maintaining immune homeostasis by mediating immune complex clearance, regulating phagocytosis, balancing cytokine production, and modulating T-cell subsets. Dysfunction of C1q has been directly linked to the onset and progression of autoimmune diseases such as SLE and RA, underscoring its crucial function in immune tolerance and disease susceptibility (26).
Elevated C1q expression correlates strongly with disease activity and inflammatory severity in autoimmune conditions, including RA (27) and Takayasu arteritis (TA). Moreover, fluctuations in C1q levels have emerged as a promising biomarker for monitoring the therapeutic efficacy of immunosuppressive treatments (28). Notably, during the A1T stage, transcripts of C1QA, C1QB, and C1QC were significantly upregulated relative to controls, suggesting that immune hyperactivation following sintilimab administration might increase susceptibility to thyroid irAEs.
These observations illuminated the dual role of C1q in autoimmune pathophysiology and treatment response, emphasizing the necessity for careful endocrine irAE monitoring in patients undergoing ICI therapy.
Based on the constructed PPI network and the outcomes of GO and KEGG enrichment analyses, we proposed a pathogenic cascade underpinning thyroid irAEs, driven by coordinated gene functional disruptions. Downregulation of FLT1 undermines vascular endothelial stability (29), which, in concert with decreased TEK expression, perturbs ANGPT/TIE2 signaling homeostasis. This destabilizes VE-cadherin-mediated adherens junctions, compromising endothelial barrier integrity and markedly increasing vascular permeability (30). Simultaneously, upregulation of PDGFRB exacerbates barrier dysfunction, facilitating vascular leakage and inflammatory cell infiltration.
Concurrently, elevated expression of C1QA, C1QB, and C1QC activates the classical complement cascade, promoting formation of C3 convertase (31) and membrane attack complexes (MAC) that directly induce lysis of thyroid follicular cells (32). Downregulation of keratins further destabilizes the cellular cytoskeleton (33), synergizing with complement-mediated damage to release autoantigens such as thyroglobulin (Tg) and thyroid peroxidase (TPO).
Reduced SPP1 expression impairs macrophage-mediated clearance of apoptotic debris (34), intensifying antigen exposure, while disrupting the Treg/Th17 balance (35), thereby diminishing immunosuppressive control and amplifying pro-inflammatory responses. The resulting autoantigens are efficiently presented by antigen-presenting cells with upregulated HLA-DPB1, triggering activation of autoreactive CD4+ T cells. These T cells secrete interferon-gamma (IFN-γ), which further upregulates HLA-DPB1 expression (36) and complement production, establishing a self-perpetuating feedback loop of antigen presentation and inflammatory amplification. This cascade culminates in relentless thyroid autoimmune injury (Figure 10).
Figure 10. Integrated pathogenic cascade of thyroid irAEs: from vascular dysfunction to autoimmune amplification.
Leveraging STRING-based analysis of eQTL-targeted genes, we identified 12 genes that directly interacted with 11 candidate genes (FLT1, PDGFRB, TEK, C1QA, C1QB, C1QC, HLA-DPB1, SPP1, KRT72, KRT73, KRT18) implicated in thyroid irAE pathogenesis. These 12 genes were linked to 75 SNPs. Notably, C3 emerged as a central bridging gene, engaging directly with four key mechanistic genes, SPP1, C1QA, C1QB, and C1QC, while genetically associated with three specific SNPs (rs189966229, rs9749508, and rs148292769). This revealed a novel genetic-immunological axis centered on C3-associated eQTLs, orchestrating complement-mediated tissue injury as a core driver of thyroid irAEs.
The identified SNPs modulated C3 expression, positioning it as a pivotal node that translated genetic susceptibility into downstream effector pathways. Through its interactions with SPP1 and the complement initiators C1QA/B/C, C3 formed a “complement-matrix axis” that promoted macrophage polarization toward a pro-inflammatory phenotype, thereby amplifying pathogenic immune responses.
Conclusion
In summary, our longitudinal transcriptomic investigation delineated dynamic gene expression profiles associated with thyroid irAEs during sintilimab therapy. By integrating genomic variants with transcriptomic data via eQTL analysis, we constructed a comprehensive genetic framework underlying thyroid irAEs, illuminating the intricate links between genetic variation and clinical phenotypes. Notably, this study is the first to identify C3 eQTLs as a central regulatory nexus bridging genetic susceptibility and complement-driven thyroid irAEs. Therapeutic targeting of this pathway holds promise for dissociating antitumor efficacy from autoimmune toxicity, offering a refined strategy for precision immunotherapy.
Despite these scientifically and clinically impactful findings, limitations remain. Restricted by sample size and follow-up duration, validation in larger, longitudinal clinical cohorts is essential. Furthermore, mechanistic insights into SNP functionality require dedicated in vitro and in vivo studies to unravel the causal pathways linking genetic variation to irAE pathogenesis.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the the corresponding authors.
Ethics statement
The studies involving humans were approved by Ethics Committee of National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
WC: Writing – review & editing, Formal analysis, Methodology, Writing – original draft. MZ: Data curation, Writing – review & editing. TL: Writing – review & editing. BS: Writing – review & editing. HS: Formal analysis, Writing – review & editing. YS: Writing – review & editing. YL: Writing – review & editing, Supervision, Conceptualization. FY: Methodology, Supervision, Investigation, Writing – review & editing, Conceptualization. GL: Conceptualization, Writing – review & editing, Formal analysis, Supervision, Methodology.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was supported by Beijing Hope Run Special Fund of Cancer Foundation of China (Grant/Award Number: LC2020L03).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1. Luciani A, Ghidini A, Borgonovo K, Parati MC, and Petrelli F. Outcome of non-small-cell lung cancer with driven mutations treated with anti-PD-(L)1 agents: A systematic review. Tumori. (2023) 109(5):442–9. doi: 10.1177/03008916221122601
2. Greillier L, Tomasini P, and Barlesi F. The clinical utility of tumor mutational burden in non-small cell lung cancer. Transl Lung Cancer Res. (2018) 7:639–46. doi: 10.21037/tlcr.2018.10.08
3. Zhang L, Mai W, Jiang W, and Geng Q. Sintilimab: A promising anti-tumor PD-1 antibody. Front Oncol. (2020) 10:594558. doi: 10.3389/fonc.2020.594558
4. Xie J, Wu X, Wu J, Huang F, and Xu L. Meta−analysis of the efficacy and safety of sintilimab for treating advanced non−small cell lung cancer. Oncol Lett. (2022) 24(6):425. doi: 10.3892/ol.2022.13545
5. Zhang Q, Jiao X, and Lai X. Clinical characters and influence factors of immune checkpoint inhibitor-related thyroid dysfunction. J Clin Endocrinol Metab. (2023) 108:2916–23. doi: 10.1210/clinem/dgad260
6. Muir CA, Clifton-Bligh RJ, Long GV, Scolyer RA, Lo SN, Carlino MS, et al. Thyroid immune-related adverse events following immune checkpoint inhibitor treatment. J Clin Endocrinol Metab. (2021) 106:E3704–13. doi: 10.1210/clinem/dgab263
7. Lu D, Yao J, Yuan G, Gao Y, Zhang J, and Guo X. Immune checkpoint inhibitor-related new-onset thyroid dysfunction: A retrospective analysis using the US FDA adverse event reporting system. Oncologist. (2022) 27:E126–32. doi: 10.1093/oncolo/oyab043
8. Hou YCC, Neidich JA, Duncavage EJ, Spencer DH, and Schroeder MC. Clinical whole-genome sequencing in cancer diagnosis. Hum Mutat. (2022) 43:1519–30. doi: 10.1002/humu.24381
9. Wang G, Xu Y, Wang Q, Chai Y, Sun X, Yang F, et al. Rare and undiagnosed diseases: From disease-causing gene identification to mechanism elucidation. Fundam Res. (2022) 2:918–28. doi: 10.1016/j.fmre.2022.09.002
10. Barsh GS, Copenhaver GP, Gibson G, and Williams SM. Guidelines for genome-wide association studies. PloS Genet. (2012) 8:7–8. doi: 10.1371/journal.pgen.1002812
11. Abbas M, Diallo A, Goodney G, and Gaye A. Leveraging the transcriptome to further our understanding of GWAS findings: eQTLs associated with genes related to LDL and LDL subclasses, in a cohort of African Americans. Front Genet. (2024) 15:1345541. doi: 10.3389/fgene.2024.1345541
12. Cano-Gamez E and Trynka G. From GWAS to function: using functional genomics to identify the mechanisms underlying complex diseases. Front Genet. (2020) 11:424. doi: 10.3389/fgene.2020.00424
13. Wilson DR, Ibrahim JG, and Sun W. Mapping tumor-specific expression QTLs in impure tumor samples. J Am Stat Assoc. (2020) 115:79–89. doi: 10.1080/01621459.2019.1609968
14. Peters JE, Lyons PA, Lee JC, Richard AC, Fortune MD, Newcombe PJ, et al. Insight into Genotype-Phenotype Associations through eQTL Mapping in Multiple Cell Types in Health and Immune-Mediated Disease. PloS Genet. (2016) 12(3):e1005908. doi: 10.1371/journal.pgen.1005908
15. Barbara E and Stranger PLDJ. Coordinating GWAS results with gene expression in a systems immunologic paradigm in autoimmunity. Curr Opin Immunol. (2012) 24:544–51. doi: 10.1016/j.coi.2012.09.002.Coordinating
16. Bogaards FA, Gehrmann T, Beekman M, Lakenberg N, Suchiman HED, de Groot CPGM, et al. Secondary integrated analysis of multi-tissue transcriptomic responses to a combined lifestyle intervention in older adults from the GOTO nonrandomized trial. Nat Commun. (2023) 15(1):7013. doi: 10.1038/s41467-024-50693-3
17. Godoy P and Reif R. Transcription factors controlling responses to toxic chemicals. Arch Toxicol. (2013) 87:3–4. doi: 10.1007/s00204-012-0981-5
18. Itoi K, Motoike I, Liu Y, Clokie S, Iwasaki Y, Uchida K, et al. Genome-Wide Analysis of Glucocorticoid-Responsive Transcripts in the Hypothalamic Paraventricular Region of Male Rats. Endocrinology. (2019) 160(1):38–54. doi: 10.1210/en.2018-00535
19. Fuxjager MJ, Lee JH, Chan TM, Bahn JH, Chew JG, Xiao X, et al. Research resource: Hormones, genes, and athleticism: Effect of androgens on the avian muscular transcriptome. Mol Endocrinol. (2016) 30:254–71. doi: 10.1210/me.2015-1270
20. Belmonte-Mateos C and Pujades C. From cell states to cell fates: how cell proliferation and neuronal differentiation are coordinated during embryonic development. Front Neurosci. (2022) 15:781160. doi: 10.3389/fnins.2021.781160
21. Jung S, Ha J, Park JH, and Yoo KH. Decoding SPP1 regulation: Genetic and nongenetic insights into its role in disease progression. Mol Cells. (2025) 48:100215. doi: 10.1016/j.mocell.2025.100215
22. White FJ, Burghardt RC, Hu J, Joyce MM, Spencer TE, and Johnson GA. Secreted phosphoprotein 1 (osteopontin) is expressed by stromal macrophages in cyclic and pregnant endometrium of mice, but is induced by estrogen in luminal epithelium during conceptus attachment for implantation. Reproduction. (2006) 132:919–29. doi: 10.1530/REP-06-0068
23. Zhao Y, Huang Z, Gao L, Ma H, and Chang R. Osteopontin/SPP1: a potential mediator between immune cells and vascular calcification. Front Immunol. (2024) 15:1395596. doi: 10.3389/fimmu.2024.1395596
24. Nedjadi T, Ahmed ME, Ansari HR, Aouabdi S, and Al-Maghrabi J. Identification of SPP1 as a prognostic biomarker and immune cells modulator in urothelial bladder cancer: A bioinformatics analysis. Cancers (Basel). (2023) 15:5704. doi: 10.3390/cancers15235704
25. Gao W, Liu D, Sun H, Shao Z, Shi P, Li T, et al. SPP1 is a prognostic related biomarker and correlated with tumor-infiltrating immune cells in ovarian cancer. BMC Cancer. (2022) 22:1–22. doi: 10.1186/s12885-022-10485-8
26. Thielens NM, Tedesco F, Bohlson SS, Gaboriaud C, and Tenner AJ. C1q: A fresh look upon an old molecule. Mol Immunol. (2017) 89:73–83. doi: 10.1016/j.molimm.2017.05.025
27. Olsen NJ, Ho E, and Barats L. Clinical correlations with serum C1q levels in patients with rheumatoid arthritis. Arthritis Rheumatol. (1991) 34:187–91. doi: 10.1002/art.1780340209
28. Chen S, Luan H, He J, Wang Y, Zeng X, Li Y, et al. Serum C1q concentration is associated with disease activity in Chinese Takayasu arteritis patients: A case-control study. Heal Sci Rep. (2021) 4. doi: 10.1002/hsr2.252
29. Otowa Y, Moriwaki K, Sano K, Shirakabe M, Yonemura S, Shibuya M, et al. Flt1/VEGFR1 heterozygosity causes transient embryonic edema. Sci Rep. (2016) 6:1–8. doi: 10.1038/srep27186
30. Cleuren A and Molema G. Organotypic heterogeneity in microvascular endothelial cell responses in sepsis—a molecular treasure trove and pharmacological Gordian knot. Front Med. (2023) 10:1252021. doi: 10.3389/fmed.2023.1252021
31. Seifert L, Zahner G, Meyer-Schwesinger C, Hickstein N, Dehde S, Wulf S, et al. The classical pathway triggers pathogenic complement activation in membranous nephropathy. Nat Commun. (2023) 14(1):473. doi: 10.1038/s41467-023-36068-0
32. Zhao C, Yu Y, Liu J, Lu G, Li T, Gao Y, et al. Diversity of complement activation in different thyroid diseases. Int Immunopharmacol. (2022) 106:108636. doi: 10.1016/j.intimp.2022.108636
33. Uchiumi A, Yamashita M, and Katagata Y. Downregulation of keratins 8, 18 and 19 influences invasiveness of human cultured squamous cell carcinoma and adenocarcinoma cells. Exp Ther Med. (2012) 3:443–8. doi: 10.3892/etm.2011.413
34. Yang X, Liu Z, Zhou J, Guo J, Han T, Liu Y, et al. SPP1 promotes the polarization of M2 macrophages through the Jak2/Stat3 signaling pathway and accelerates the progression of idiopathic pulmonary fibrosis. Int J Mol Med. (2024) 54:1–12. doi: 10.3892/ijmm.2024.5413
35. Zheng Y, Zhao L, Xiong Z, Huang C, Yong Q, Fang D, et al. Ursolic acid targets secreted phosphoprotein 1 to regulate Th17 cells against metabolic dysfunction-associated steatotic liver disease. Clin Mol Hepatol. (2024) 30:449–67. doi: 10.3350/cmh.2024.0047
Keywords: immune-related adverse events (irAEs), sintilimab, eQTL analysis, complement system, longitudinal transcriptomics
Citation: Chen W, Zhang M, Li T, Shang B, Su H, Shi Y, Liu Y, Yu F and Li G (2025) Mapping the genetic–transcriptional landscape of thyroid irAEs in sintilimab therapy: toward biomarker-guided immunotoxicity prediction. Front. Immunol. 16:1671594. doi: 10.3389/fimmu.2025.1671594
Received: 23 July 2025; Accepted: 07 October 2025;
Published: 27 October 2025.
Edited by:
Andrea Corsello, Agostino Gemelli University Polyclinic (IRCCS), ItalyReviewed by:
Dmitry Aleksandrovich Zinovkin, Gomel State Medical University, BelarusHesong Wang, Fourth Hospital of Hebei Medical University, China
Copyright © 2025 Chen, Zhang, Li, Shang, Su, Shi, Liu, Yu and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yutao Liu, bGl1eXV0YW9AY2ljYW1zLmFjLmNu; Feng Yu, eXVmZW5nQGNwdS5lZHUuY24=; Guohui Li, bGdoMDYwM0BjaWNhbXMuYWMuY24=
†These authors have contributed equally to this work
Mingyu Zhang2†