Xueying Liu
Yahu Miao
Qiu Zhang- Department of Endocrinology, First Affiliated Hospital of Anhui Medical University, Hefei, China
Background: Autoimmune thyroid disease (AITD) is induced by various factors, including inheritability, which regulates gene expression. Multiple loci correlated with AITD have been discovered utilizing genome-wide association studies (GWASs). Nevertheless, demonstrating the biological relevance and function of these genetic loci is difficult.
Methods: The FUSION software was utilized to define genes that were expressed differentially in AITD using a transcriptome-wide association study (TWAS) method in accordance with GWAS summary statistics from the largest genome-wide association study of 755,406 AITD individuals (30,234 cases and 725,172 controls) and levels of gene expression from two tissue datasets (blood and thyroid). Further analyses were performed such as colocalization, conditional, and fine-mapping analyses to extensively characterize the identified associations, using functional mapping and annotation (FUMA) to conduct functional annotation of the summary statistics of 23329 significant risk SNPs (P < 5 × 10−8) recognized by GWAS, together with summary-data-based mendelian randomization (SMR) for identifying functionally related genes at the loci in GWAS.
Results: There were 330 genes with transcriptome-wide significant differences between cases and controls, and the majority of these genes were new. 9 of the 94 unique significant genes had strong, colocalized, and potentially causal correlations with AITD. Such strong associations included CD247, TPO, KIAA1524, PDE8B, BACH2, FYN, FOXK1, NKX2-3, and SPATA13. Subsequently, applying the FUMA approach, novel putative AITD susceptibility genes and involved gene sets were detected. Furthermore, we detected 95 probes that showed strong pleiotropic association with AITD through SMR analysis, such as CYP21A2, TPO, BRD7, and FCRL3. Lastly, we selected 26 genes by integrating the result of TWAS, FUMA, and SMR analysis. A phenome-wide association study (pheWAS) was then carried out to determine the risk of other related or co-morbid phenotypes for AITD-related genes.
Conclusions: The current work provides further insight into widespread changes in AITD at the transcriptomic level, as well as characterized the genetic component of gene expression in AITD by validating identified genes, establishing new correlations, and uncovering novel susceptibility genes. Our findings indicate that the genetic component of gene expression plays a significant part in AITD.
1 Introduction
Autoimmune thyroid disease (AITD), which comprises Graves’ disease (GD) and Hashimoto’s thyroiditis (HT), is induced by an immune system dysfunction that results in an immune response on the thyroid. It is one of the most typical autoimmune conditions, which is genetically heterogeneous (1). According to epidemiological data regarding genetic susceptibility to AITD, AITD concordance is greater than 50% in monozygotic twins (2). The vast majority of AITD cases go undetected until the immunological response damages the thyroid gland. This is in complete contradiction to other prevalent autoimmune conditions that frequently co-occur with AITD (3), such as rheumatoid arthritis (RA), where several risk loci have been found, multiple effective treatments are obtainable, and an effort is being made on early detection and even disease prevention (4, 5). With advances in comprehending the immunological and molecular basis of autoimmune conditions, gene therapy has emerged as a potential treatment option for patients who suffer from the disease (6). Therefore, the demand for prevention and early treatment of AITD is necessary, as well as a thorough knowledge of disease pathogenesis.
Genome-wide association studies (GWASs) have been frequently used, and they constitute an important step in identifying the particular genes that underlie AITD genetic risk. Six loci with evidence of AITD linkage were revealed in a genome-wide study of 56 multigenerational AITD families (354 people) utilizing 387 microsatellite markers (7). In addition, Chu et al. (8) performed a GWAS in 1,536 patients with GD (cases) and 1,516 controls and then examined a group of associated single nucleotide polymorphisms (SNPs) in a second sample of 3,994 cases and 3,510 controls. They identified two novel susceptibility loci and corroborated four previously described loci (in the major histocompatibility complex, TSHR, CTLA4, and FCRL3) (the RNASET2-FGFR1OP-CCR6 region at 6q27 and an intergenic region at 4p14). Oryoj et al. (9) discovered a new SNP at the VAV3 locus associated with HT using a sample of 181 cases and 1363 healthy controls. Even though GWAS has allowed the discovery of SNPs that confer susceptibility to AITD, the functional significance of these genetic variants is unknown. Examining intermediary processes between the genome and the phenotype, like gene expression, may assist us in comprehending the molecular mechanisms behind AITD.
TWAS, one method for prioritizing causative genes at GWAS loci, could use expression reference panels (eQTL cohorts with expression and genotype data) to identify gene-trait associations from GWAS datasets (10–12). Furthermore, Zhu et al. (11) introduced Summary-data-based Mendelian Randomization (SMR), which may incorporate summary-level data from GWAS with eQTL data to find genes related to a complex trait due to pleiotropy at the loci found in GWAS. Using this method, previous research has recently discovered novel trait-associated genes for 33 human complex traits (11, 13). To date, these methods applied to AITD have not been well documented.
Herein, we carried out a TWAS to find genetically regulated genes correlated with AITD. Using conditional analysis, we revealed that genetically regulated expression generates several genome-wide significant signals from the AITD GWAS. Then we investigate how these genes induce dysregulation of neighboring co-expressed genes. Furthermore, Functional Mapping and Annotation (FUMA) was performed to make a functional annotation for significant risk SNPs identified by GWAS. SMR analysis was then employed to identify significant genes that were pleiotropically related to AITD. Finally, we performed a PheWAS and identified the risk of other associated or co-morbid phenotypes for AITD-related genes. This present research will advance the field by elucidating the molecular mechanisms of AITD, including several disease-relevant tissues, and identifying critical genes for the disease that might be future prospects for investigation on its therapy and etiology.
2 Materials and methods
2.1 Datasets
The analysis used 1) genome-wide summary statistics from the GWAS of AITD by Saevarsdottir et al. (14), 2) 4 SNP weight sets from 3 separate transcriptomic reference samples (GTEx v8, NTR (Netherlands Twin Register) and YFS (Young Finns Study), and 3) the 1000 Genomes Project linkage disequilibrium (LD) reference.
To begin, we obtained the GWAS data from 755,406 European participants, including 30,234 cases and 725,172 controls, and the AITD summary statistic from the deCODE Genetics SUMMARY DATA (14), containing genetic susceptibility information to AITD for 36,004,550 HapMap3 SNPs from two datasets, Iceland and the UK Biobank.
In addition, SNP weights from three tissues were utilized: peripheral blood, whole blood, and thyroid. SNP weights indicate the associations between SNPs and the expression of their associated gene. Four SNP weight sets were obtained from the TWAS FUSION official website (http://gusevlab.org/projects/fusion/#reference-functional-data). The weights were assigned to four RNA reference samples as follows: The Netherlands Twin Register (NTR) and the Young Finns Study (YFS) both yield gene expression data on blood tissue, and the GTEx Consortium, a cutting-edge study that examines expression in whole blood and thyroid separately. Finally, the FUSION website was used to obtain the 1000 Genomes Phase 3 European LD reference, which included 489 individuals (http://gusevlab.org/projects/fusion/).
2.2 Statistical analysis
Ubuntu 22/Bash (GNU Project Bourne Again Shell) and R version 4.2.1 (The R Project for Statistical Computing, Vienna, Austria Anhui, China) were utilized to conduct all statistical analyses.
2.3 Transcriptome-wide significance threshold
To calculate the transcriptome-wide significance threshold for this investigation, we employed a permutation approach from a previous TWAS (15). This method calculates a significant threshold that is modified for the number of evaluated features and takes feature correlation within and across SNP weight sets into account. A strict Bonferroni-corrected study-wise threshold was set: P = 2.38 × 10-6 (0.05/21,053) (total quantity of genes across panels).
2.4 TWAS FUSION and colocalization
FUSION is a series of programs for carrying out transcriptome-wide and regulome-wide association studies (TWAS and RWAS). FUSION develops predictive models of the genetic component of a functional/molecular phenotype and analyzes disease-associated components using GWAS summary statistics (16). TWAS analysis was carried out in this study utilizing software that followed the TWAS FUSION protocol with default settings (16). In addition, the colocalization test, which gives data of a shared causal variant between the predicted functional feature and the trait, would be an alternative to the TWAS test (which looks for a substantial correlation between the predicted functional feature and the trait). Colocalization was conducted using the coloc R package with a 0.5Mb window for all genes meeting transcriptome-wide significance (P < 0.05), and FUSION was used to evaluate the posterior probability (PP) that GWAS and TWAS associations share a causal SNP. The coloc statistic for each of the coloc hypotheses is represented by the PP0 (no association), PP1 (functional association only), PP2 (GWAS association only), PP3 (independent functional/GWAS associations), and PP4 (colocalized functional/GWAS associations) columns.
2.5 Conditional analysis
Using the FUSION software, conditional and joint analyses were conducted for transcriptome-wide significant (Bonferroni-corrected criteria) signals to confirm whether these significant signals were attributable to multiple-associated features or conditionally independent (16). The software assesses the influence of all important features within each locus by examining residual SNP correlations with AITD after accounting for the predicted expression of other features (16). When other features in the region’s predicted expression are taken into consideration, the approach may identify which features indicate independent associations, referred to as jointly significant, and which do not, referred to as marginally significant (16). For post-processing and creating numerous conditional output plots as well as summary statistics, the “FUSION.post process.R” script was utilized.
2.6 Statistical fine mapping
Mancuso N et al. (17) introduced FOCUS, a fine-mapping method that may estimate credible sets of causal genes from TWAS utilizing prediction eQTL weights, LD, and GWAS summary statistics. FOCUS assesses the posterior inclusion probability (PIP) of each feature being causal in a region of association, similar to statistical fine mapping of GWAS results, and uses the sum of PIP to construct the default 90% credible set, a set of features expected to contain the causal features. The PIP of individual features with PIP > 0.5 indicates that a feature has a higher probability of being causal than any other feature in the region. Without using the tissue prioritization option, the FOCUS fine mapping function was performed on all SNP weight panels at the same time.
2.7 Functional mapping and annotation of GWAS
2.7.1 Definition of genomic risk loci based on GWAS
The AITD meta-analysis, which included 30,234 cases and 725,172 controls, identified 23,329 SNPs with P < 5.0 × 10-8 (14). For identifying genomic risk loci for AITD susceptibility according to previous GWAS summary statistics (14), LD structure from the European ancestry of the 1000 Genomes Reference panel was computed and utilized in FUMA (18). Based on previous literature, genomic risk loci, independent significant SNPs, lead SNPs, and candidate SNPs were defined (18).
2.7.2 Annotation of candidate SNPs in genomic risk loci
FUMA is a post-GWAS annotation program that uses candidate SNPs to identify possible causal variants for AITD. The FUMA platform was used to develop a number of SNP functional annotation approaches (18). ANNOVAR (19) (‘gene-based annotation’) is used to examine the functional consequences of all SNPs on genes utilizing Ensembl genes (build 85). Note that candidate SNPs were annotated for functional consequences (3′ untranslated region, 5′ untranslated region, downstream, exonic, intergenic, intronic, ncRNA_exonic, ncRNA_intronic, ncRNA_splicing, splicing, upstream across genes) on gene functions.
2.7.3 MAGMA for gene analysis and gene-set analysis
FUMA utilizes the MAGMA tool to generate gene-based P-values (gene analysis) and gene set P-values (gene set analysis) from input AITD GWAS summary statistics (20). The major histocompatibility complex region and indels were eliminated from the analysis. For gene analysis, all SNPs within genes were mapped to 19,987 protein-coding genes. The default significant genes threshold is defined as P = 2.50 × 10−6 (the Bonferroni correction, 0.05/19,987). Individual genes that share particular biological activities were pooled and studied for correlation with disease or trait in gene set analysis. The gene set P-value is calculated utilizing the same threshold criteria (P = 0.05/15,488).
2.8 SMR analysis
SMR analyses were carried out to observe genes that were causally associated with AITD, allowing us to prioritize functionally relevant genes in the GWAS loci (11). SMR incorporates GWAS and eQTL summary statistics to assess for a pleiotropic correlation between gene expression and a trait owing to a shared variant at a locus using Mendelian Randomization (MR) principles.
The cis-eQTL genetic variants were utilized as the instrumental variables (IVs) for gene expression in the SMR analysis. For gene expression in blood and thyroid, SMR analysis was conducted separately. Meanwhile, for blood, we employed the Westra, CAGE and GTEx v8 eQTL summarized data (21–23), which comprised 3,511, 2,765, and 670 participants respectively. In addition, we utilized the version 8 release of the GTEx eQTL summarized data (23), which was for thyroid tissue with 574 participants. The blood and thyroid eQTL data may be downloaded online (https://yanglab.westlake.edu.cn/software/smr/#DataResource).
To assess the presence of LD in the identified association, the heterogeneity in dependent instruments (HEIDI) test was conducted. The genome-wide significance level for the SMR test is defined as 0.05/the number of probes. The significance values of these probes are P < 5.89 × 10-6 (peripheral blood from CAGE.sparse), P < 8.44 × 10-6 (peripheral blood from westra_eqtl_hg19), P < 7.65 × 10-6 (whole blood from GETx), and P < 5.57 × 10-6 (thyroid from GETx) respectively and those probes with little evidence of heterogeneity (PHEIDI ≥ 0.05) were retained.
2.9 Phenome-wide association studies
A phenome-wide association study (PheWAS) was undertaken utilizing publicly accessible data from the GWAS Atlas (https://atlas.ctglab.nl) to identify phenotypes associated with the AITD genes via TWAS. The top ten phenotypes were reported.
3 Result
3.1 Transcriptome-wide significant hits
We discovered 499 significant features from 330 unique genes that were differentially expressed (P < 2.38 × 10-6) across several SNP weight sets in AITD (Figure 1; Supplementary Table S1). Among the 499 significant features, 237 were upregulated, while 262 were downregulated. In comparison to the previous AITD GWAS by Saevarsdottir S et al, Chu et al, and Oryoji et al. (8, 9, 14), 296 unique genes were novel and 34 were previously implicated. The largest number of associations were from the GTEx thyroid set (209 associated features).
Figure 1 The relationship between gene expression and autoimmune thyroid disease. Manhattan-style plot of z scores for each of the tested genes, across all autosomes and tested single nucleotide polymorphism weight sets. Blue lines indicate the transcriptome-wide significance threshold. All the statistically significant genes are shown.
3.2 Colocalization analysis
We calculated the PP that the genetic and functional correlations generated by various causal SNPs (PP3) and a shared causal SNP (PP4) to evaluate the colocalization status of the gene. 159 (118 unique genes) of the 499 significant features were regarded as colocalized based on their high PP4 content (PP4 > 0.8). It was observed that 12 significant features (9 unique genes) explain all of the signals at their loci (PP4 = 1) (Supplementary Table S2).
3.3 Conditional analysis
We discovered that numerous significant features were located inside the same locus (defined as a 0.5Mb window). Conditional analysis of the 499 significant features revealed 135 jointly (129 unique genes) and 364 marginally significant features (Supplementary Table S3), indicating that many identified genes were linked to AITD owing to their co-expression with the 135 independent features.
We further examined how adjusting for gene expression in features impacted the correlations between SNPs and traits. The percentage of variance in GWAS correlations accounted for gene expression correlations at a specific location ranges from 4.94% to 100%. Genome-wide association signals on Chromosomes 1, 2, 4, 8, 11, 12, 13, 16, 17, 18 and 22 fully explained relative genetic variations (Supplementary Table S4). This might imply that genetic risk for AITD is mediated by functional feature alterations at these loci.
3.4 TWAS fine mapping
FOCUS was conducted to compute PIP for each feature (17). We found 16 significant features (16 unique genes) with PIP > 0.5 in three tissues, suggesting that these are likely causal in their correlations with AITD. Of these, 9 genes were supported by the colocalization analysis (Table 1). The genes with the highest possibility of causality were SH2B3 (YFS blood), SPATA13 (GTEx thyroid), TPO (GTEx thyroid), IRF4 (GTEx thyroid), and BACH2 (NTR blood) (PIP = 1.00) (Table 1). Furthermore, we identified multiple significant features from other tissues, such as adipose, artery, brain, liver, and so forth. The detailed information is presented in Supplementary Table S5.
Table 1 Causal posterior probabilities for genes in 90%-credible sets for AITD TWAS signals.
3.5 FUMA of GWAS
3.5.1 Functional annotation analysis by FUMA
The FUMA approach was utilized to conduct functional annotation of AITD GWAS summary statistics, which included 30,234 cases and 725,172 controls. Starting with the AITD GWAS summary statistics, 259 lead SNPs and 23,329 candidate SNPs with LD with the lead SNPs were discovered among 955 independent significant SNPs across 128 genomic loci (Figure 2; Supplementary Figure 1). We totally prioritized 583 unique genes from 128 loci through FUMA (Supplementary Table S6, S7).
Figure 2 The consequences of candidate SNPs functionally annotated. *P < 0.05.
3.5.2 Gene analysis, gene-set analysis
There were 424 significant genes found (Supplementary Table S8), and 15,488 gene sets were identified. 142 of these gene sets were found to be significant (P < 3.23 ×10-6). The most significant gene set was positive regulation of the immune system process (P = 7.18×10-22), followed by regulation of the immune system process (P = 1.19 × 10-22). Of these pathways, 50 gene sets associated with AITD were tested, and 159 were associated with other autoimmune conditions such as systemic lupus erythematosus (SLE) and type 1 diabetes mellitus (T1DM) (Supplementary Table S9).
3.6 SMR analysis
3.6.1 Basic information of the summarized data
According to the SMR analysis, The number of probes was 20,943 (8485 from CAGE.sparse, 5923 from westra_eqtl_hg19, and 6535 from GETx whole blood) and 8,972 in the analysis for blood and thyroid, respectively.
3.6.2 SMR analysis in blood
We found 71 probes that were significantly associated with AITD using CAGE, Westra, and GTEx eQTL summarized data (Table 2). Using AITD GWAS summarized data, we found a few probes tagging FCRL3/FCRH3 (β [SE] = -0.055 [0.011], PSMR = 2.05 ×10-7, PHEIDI = 0.350, CAGE eQTL study; β [SE] = -0.077 [0.015], PSMR = 2.03 ×10-7, PHEIDI = 0.082, Westra eQTL study; β [SE] = -0.158 [0.032], PSMR = 7.57 ×10-7, PHEIDI = 0.089, GTEx eQTL study) (Figures 3A, B) and BRD7 (β [SE] = 0.096 [0.017], PSMR = 2.44 ×10-8, PHEIDI = 0.078, CAGE eQTL study; β [SE] = 0.157 [0.029], PSMR = 3.87 ×10-8, PHEIDI = 0.082, Westra eQTL study; β [SE] = 0.102 [0.019], PSMR = 3.35 ×10-8, PHEIDI = 0.105, GTEx eQTL study), that showed strong pleiotropic association with AITD (Figure 4). Particularly, the most significant pleiotropic associations with AITD were detected on CYP21A2 (β [SE] = -0.339 [0.040], PSMR = 1.41 ×10-17, PHEIDI = 0.079, GTEx eQTL study).
Table 2 Significant genes showing pleiotropic association with AITD in blood.
Figure 3 Prioritizing genes around FCRL3/FCRH3 in association with AITD in blood. Top plot, gray dots represent the -log 10 (P values) for SNPs from the GWAS of AITD, and rhombuses represent the -log 10 (P values) for probes from the SMR test with solid rhombuses indicating that the probes pass HEIDI test and hollow rhombuses indicating that the probes do not pass the HEIDI test. Middle plot, eQTL results in blood for three probes, ILMN_1691693, ILMN_1797428, ILMN_1699599 (A) and ENSG00000160856 (B) probes, tagging FCRL3/FCRH3. Bottom plot, location of genes tagged by the probes. Highlighted in maroon indicates probes that pass SMR threshold. FCRL3, Fc Receptor Like 3; FCRH3, Fc Receptor Homolog 3; AITD, autoimmune thyroid disease; GWAS, genome-wide association studies; SMR, summary-data-based Mendelian randomization; HEIDI, heterogeneity in dependent instruments; eQTL, expression quantitative trait loci.
Figure 4 Prioritizing genes around BRD7 in association with AITD in blood. Top plot, gray dots represent the -log 10 (P values) for SNPs from the GWAS of AITD, and rhombuses represent the -log 10 (P values) for probes from the SMR test with solid rhombuses indicating that the probes pass HEIDI test and hollow rhombuses indicating that the probes do not pass the HEIDI test. Middle plot, eQTL results in blood for two probes, ILMN_2082810, and ILMN_1696420 tagging BRD7. Bottom plot, location of genes tagged by the probes. Highlighted in maroon indicates probes that pass SMR threshold. BRD7, Bromodomain Containing 7; AITD, autoimmune thyroid disease; GWAS, genome-wide association studies; SMR, summary-data-based Mendelian randomization; HEIDI, heterogeneity in dependent instruments; eQTL, expression quantitative trait loci.
3.6.3 SMR analysis in thyroid
The GTEx eQTL summary data revealed twenty-four genes that were pleiotropically correlated with AITD in thyroid tissue (Table 3). NR3C2 displayed the most significant pleiotropic associations with AITD (β [SE] = 0.200 [0.021], PSMR = 6.88 ×10-21, PHEIDI = 0.367), then TPO (β [SE] = -0.208 [0.025], PSMR = 8.91 ×10-17, PHEIDI = 0.886) and RNASET2 (β [SE] = 0.465 [0.061], PSMR = 3.73 ×10-14, PHEIDI = 0.107). Of note, IRF3 showed significant pleiotropic association not only in the thyroid but also in the blood (β [SE] = 0.307 [0.057], PSMR = 8.41 ×10-8, PHEIDI = 0.613; and β [SE] = 0.150 [0.025], PSMR = 3.11 ×10-9, PHEIDI = 0.156; β [SE] = 0.144 [0.023], PSMR = 2.29 ×10-10, PHEIDI = 0.074) (Figures 5A, B).
Table 3 Significant genes showing pleiotropic association with AITD in thyroid tissue.
Figure 5 Prioritizing genes IRF3 in association with AITD in blood (A) and thyroid (B). Top plot, gray dots represent the -log 10 (P values) for SNPs from the GWAS of AITD, and rhombuses represent the -log 10 (P values) for probes from the SMR test with solid rhombuses indicating that the probes pass HEIDI test and hollow rhombuses indicating that the probes do not pass the HEIDI test. Middle plot, eQTL results in blood. Bottom plot, location of genes tagged by the probes. Highlighted in maroon indicates probes that pass SMR threshold. IRF3, Interferon Regulatory Factor 3; AITD, autoimmune thyroid disease; GWAS, genome-wide association studies; SMR, summary-data-based Mendelian randomization; HEIDI, heterogeneity in dependent instruments; eQTL, expression quantitative trait loci.
3.7 The result of PheWAS
Since plenty of genes were associated with AITD, we selected 26 genes that were simultaneously present in the result of TWAS, FUMA and SMR analysis (Figure 6; Supplementary Table S10). We performed a PheWAS and identified the risk of other associated or co-morbid phenotypes for AITD-related genes, which included immunological traits (e.g., granulocyte percentage of myeloid white cells); connective tissue traits (e.g., rheumatoid arthritis); cardiovascular traits (e.g., resting heart rate). The top ten phenotypes of 26 genes are presented in Supplementary Table S11.
Figure 6 The result of an integration TWAS, FUMA and SMR analysis.
4 Discussion
AITD is a prevalent autoimmune condition that affects millions of people all over the world (2). Although recent GWAS have been effectively identifying a few risk loci for AITD, the functional significance of these correlations has remained unknown due to the difficulties in interpreting LD and inferring causality (24). TWAS takes full advantage of eQTL expression reference panels that comprise expression and genotype data to seek gene-trait correlations from GWAS datasets (10, 12, 16). The current study was the first to attempt an AITD TWAS by applying the summary statistics of 755,406 AITD individuals (30,234 cases, and 725,172 controls) from the most recent ATID GWAS to reveal the genetic components of AITD gene expression through an extensive exploration of different tissue types, and intensive characterization of identified associations representing pivotal transcriptome changes.
Here, several key findings were emphasized. First, we identified 330 genes significantly associated with AITD, the majority of which were novel (296/330). 118 of the 330 significant features were independently associated with AITD, implying that the majority of the detected features are dependent on LD and correlated predicted expression of neighboring genes. Following that, based on a set of subsequent analyses, we discovered 9 correlations with AITD due to their significant TWAS signal, colocalization (PP4), and PIP. Furthermore, using the FUMA approach, prioritized 583 unique genes from 128 loci through the MAGMA tool. Moreover, utilizing SMR analysis, we totally identified 95 genes that showed pleiotropically associated with AITD in blood and thyroid tissues. Noteworthily, we finally obtained 26 key susceptibility genes by integrating the results of TWAS, FUMA, and SMR methods. When compared to the aforementioned meta-analysis, 12 of these 26 key genes had previously been reported, while 14 were novel (14).
We confirmed previously reported genes such as BACH2, TPO, RNASET2, IRF, etc. BACH2 is a transcriptional repressor in the BACH family of basic leucine zipper transcription factors. The BACH2 protein is generated throughout B-cell maturation and is important in the development and function of both the innate and adaptive immune systems (25, 26). Thyroid peroxidase (TPO) encoded protein functions as an enzyme and is essential for thyroid gland function. AITD-associated loci, including TPO-rs11675434, and BACH2-rs10944479, were discovered in a GWAS meta-analysis including 18,297 individuals, with 1,769 TPOAb-positives and 16,528 TPOAb-negatives for TPOAb-positivity and in 12,353 individuals for TPOAb serum levels, with replication in 8,990 individuals (27). RNASET2 is a ribonuclease that is important in the innate immune response by identifying and degrading RNAs from microbial pathogens that are then detected by TLR8 (28). The gene RNASET2 has been suggested as a putative risk hotspot for GD in GWASs (8). An emerging longitudinal, case-control study suggested that RNASET2 was considerably overexpressed in GD patients in comparison to healthy controls (29). However, since the evidence is questionable due to the limited sample size, replication is necessary. The gene IRF3 showed an observably pleiotropic relationship with AITD both in blood and the thyroid. IRF3, encoding a member of the interferon regulatory transcription factor (IRF) family, functions as a key transcriptional regulator of type I interferon-dependent immune responses, which are crucial in the innate immune response (30). Previous GWASs have also identified an association between IRF3 and various autoimmune diseases such as RA, T1DM, and SLE (31–33). Additionally, coinciding with an original meta-analysis of GWAS conducted by Saevarsdottir et al. (14), we found FCRL3 was significantly associated with AITD. According to the predicted molecular structure, FCRL3 may function by encoding a membrane protein that conveys signals into cells (34). Moreover, previous studies suggested that FCRL3 might be a risk gene for RA (35, 36). The findings above, together with ours, indicated that FCRL3 plays a significant role in the immune system, emphasizing its enormous potential as a prospective target for AITD treatment and prevention.
Among these 14 novel genes identified, CYP21A2, MSH5, PPP1R18, and MMEL1 are found in or near the human major histocompatibility complex (MHC) region on chromosome 6, where these genes play a crucial function in immune response and immunological regulation and are implicated in several inflammatory and autoimmune diseases (37–42). Previous GWASs identified plenty of genetic variants in the MHC region correlated with AITD (8, 9). In blood, we found that CYP21A2 had the most significant pleiotropic associations with AITD. CYP21A2, as well as MSH5, PPP1R18, and NEU1 have been demonstrated by previous GWASs associated with other autoimmune disorders such as RA, ankylosing spondylitis (AS), T1DM, and SLE (31–33, 43, 44). In addition, the three genes MMEL1, CUTA, and CCDC88B were found that have significant relation with RA (32), the gene SUOX with RA as well as T1DM (31, 32), and UHRF1BP1 with SLE (33). Thyroid disorders have frequently occurred in patients with systemic autoimmune conditions including RA and SLE (45, 46). For example, genetic variations were presented in a research of 35 families with numerous cases of SLE accompanied by AITD, in which a gene susceptibility locus was detected in 5q14.3-q15, the major susceptibility locus for SLE, also observed in AITD (47). To some extent, such a result suggested a potential genetic link between SLE and AITD despite the small samples. The gene BRD7, in particular, was found to be significantly associated with AITD. BRD7 is a member protein of the bromodomain-containing protein family and plays a crucial role in the pathogenesis of cancers and the regulation of inflammation, metabolism, and obesity (48–50). In addition, a cross-sectional study revealed that obesity raises the risk of AITD by affecting a loop involving leptin (51). This is the first study to identify a connection between the gene BRD7 and AITD; moreover, published studies on the correlation between BRD7 and other autoimmune disorders are scarce, and more study is urgently needed.
Additionally, gene-set analysis indicated that the top five predominant pathways correlated with AITD were positive regulation of immune system process, regulation of immune system process, lymphocyte activation, T cell activation, and immune system development. There are some commonalities between AITD and other autoimmune diseases. Huber, A et al. (52) believed that all currently identified genes associated with both AITD and T1DM are involved in immune regulation, specifically in presenting antigenic peptides to T cells. Further, activated lymphocytes are often involved in the progression of autoimmune diseases such as SLE and RA (53), which has also been proven by Burbano, C et al. (54).
These findings above underline the transcriptomic alterations in AITD that probably play essential roles in disease etiology and pathogenesis, attributed to the same genetic variations as the disease and determinants, resulting in transcriptomic alterations in gene expression and traits. Such outcomes on key genes and gene-sets for AITD can guide future drug target research as well as candidate gene exploration in animal or cell research, where the consequences of molecular and its pathway changes can be easier to be explored by inducing gene knock-down or up-regulation.
Although our findings were promising, several major limitations need to be addressed. First, since only European individuals were involved in this investigation, the findings may not apply to other ethnicities. More research is required to verify the current results in specific ethnic groups, such as Asians. In addition, the validity of our findings was constrained owing to the relatively limited sample sizes of available gene expression references, suggesting that larger samples are urgently required. Furthermore, because the GWAS samples had a wide range of AITD definitions (GD or HT), future study is required to examine how restricted phenotyping influences findings at the transcriptomic level based on more accurate definitions.
In summary, we present evidence for extensive alterations in AITD at a transcriptomic level. Our study can identify new associations and characterize the transcriptomic alterations that currently recognized risk factors undergo. We also identified several genes that are likely to be essential in AITD. These results indicate that the genetic component of gene expression is crucial in AITD.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.
Author contributions
Conceptualization: XL. Methodology: XL. Formal analysis and investigation: XL and YM. Writing - original draft preparation: XL. Writing - review and editing: XL, CL, WL, QF. Funding acquisition: QZ. Resources: XL and YM. Supervision: QZ and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding
National Natural Science Foundation of China (grant No.81970703). Anhui Province Clinical Medical research transformation Project (202204295107020027).
Acknowledgments
We thank Alexander Gusev for creating the TWAS FUSION method used here. We also thank the Psychiatric Genomics Consortium Autism Spectrum Disorder Working Group for making the Autoimmune thyroid disease genome-wide association study results publicly available.
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.
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1161311/full#supplementary-material
Supplementary Figure 1 | Summary of genomic risk loci based on AITD GWAS. Y axis, Genomic risk loci are displayed by the ‘chromosome: start position-end position’. X axis, Histograms from left to right depict the size of the genomic loci, the number of candidate SNPs, mapped genes by positional mapping and eQTL mapping in the genomic loci, and the number of identified genes located within the genomic loci, respectively.
References
1. Hwangbo Y, Park YJ. Genome-wide association studies of autoimmune thyroid diseases, thyroid function, and thyroid cancer. Endocrinol Metab (Seoul) (2018) 33(2):175–84. doi: 10.3803/EnM.2018.33.2.175
2. Antonelli A, Ferrari SM, Corrado A, Di Domenicantonio A, Fallahi P. Autoimmune thyroid disorders. Autoimmun Rev (2015) 14(2):174–80. doi: 10.1016/j.autrev.2014.10.016
3. Fallahi P, Elia G, Ragusa F, Ruffilli I, Camastra S, Giusti C, et al. The aggregation between AITD with rheumatologic, or dermatologic, autoimmune diseases. Best Pract Res Clin Endocrinol Metab (2019) 33(6):101372. doi: 10.1016/j.beem.2019.101372
4. Aletaha D, Smolen JS. Diagnosis and management of rheumatoid arthritis: a review. Jama (2018) 320(13):1360–72. doi: 10.1001/jama.2018.13103
5. Cope AP. Considerations for optimal trial design for rheumatoid arthritis prevention studies. Clin Ther (2019) 41(7):1299–311. doi: 10.1016/j.clinthera.2019.04.014
6. Shu SA, Wang J, Tao MH, Leung PS. Gene therapy for autoimmune disease. Clin Rev Allergy Immunol (2015) 49(2):163–76. doi: 10.1007/s12016-014-8451-x
7. Tomer Y, Barbesino G, Greenberg DA, Concepcion E, Davies TF. Mapping the major susceptibility loci for familial graves’ and hashimoto’s diseases: evidence for genetic heterogeneity and gene interactions. J Clin Endocrinol Metab (1999) 84(12):4656–64. doi: 10.1210/jcem.84.12.6216
8. Chu X, Pan CM, Zhao SX, Liang J, Gao GQ, Zhang XM, et al. A genome-wide association study identifies two new risk loci for graves’ disease. Nat Genet (2011) 43(9):897–901. doi: 10.1038/ng.898
9. Oryoji D, Ueda S, Yamamoto K, Yoshimura Noh J, Okamura K, Noda M, et al. Identification of a hashimoto thyroiditis susceptibility locus via a genome-wide comparison with graves’ disease. J Clin Endocrinol Metab (2015) 100(2):E319–24. doi: 10.1210/jc.2014-3431
10. Gamazon ER, Wheeler HE, Shah KP, Mozaffari SV, Aquino-Michaels K, Carroll RJ, et al. A gene-based association method for mapping traits using reference transcriptome data. Nat Genet (2015) 47(9):1091–8. doi: 10.1038/ng.3367
11. Zhu Z, Zhang F, Hu H, Bakshi A, Robinson MR, Powell JE, et al. Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets. Nat Genet (2016) 48(5):481–7. doi: 10.1038/ng.3538
12. Barbeira AN, Dickinson SP, Bonazzola R, Zheng J, Wheeler HE, Torres JM, et al. Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics. Nat Commun (2018) 9(1):1825. doi: 10.1038/s41467-018-03621-1
13. Pavlides JM, Zhu Z, Gratten J, McRae AF, Wray NR, Yang J. Predicting gene targets from integrative analyses of summary data from GWAS and eQTL studies for 28 human complex traits. Genome Med (2016) 8(1):84. doi: 10.1186/s13073-016-0338-4
14. Saevarsdottir S, Olafsdottir TA, Ivarsdottir EV, Halldorsson GH, Gunnarsdottir K, Sigurdsson A, et al. FLT3 stop mutation increases FLT3 ligand level and risk of autoimmune thyroid disease. NATURE (2020) 584(7822):619–23. doi: 10.1038/s41586-020-2436-0
15. Pain O, Pocklington AJ, Holmans PA, Bray NJ, O’Brien HE, Hall LS, et al. Novel insight into the etiology of autism spectrum disorder gained by integrating expression data with genome-wide association statistics. Biol Psychiatry (2019) 86(4):265–73. doi: 10.1016/j.biopsych.2019.04.034
16. Gusev A, Ko A, Shi H, Bhatia G, Chung W, Penninx BW, et al. Integrative approaches for large-scale transcriptome-wide association studies. Nat Genet (2016) 48(3):245–52. doi: 10.1038/ng.3506
17. Mancuso N, Freund MK, Johnson R, Shi H, Kichaev G, Gusev A, et al. Probabilistic fine-mapping of transcriptome-wide association studies. Nat Genet (2019) 51(4):675–82. doi: 10.1038/s41588-019-0367-1
18. Watanabe K, Taskesen E, van Bochoven A, Posthuma D. Functional mapping and annotation of genetic associations with FUMA. Nat Commun (2017) 8(1):1826. doi: 10.1038/s41467-017-01261-5
19. Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res (2010) 38(16):e164. doi: 10.1093/nar/gkq603
20. de Leeuw CA, Neale BM, Heskes T, Posthuma D. The statistical properties of gene-set analysis. Nat Rev Genet (2016) 17(6):353–64. doi: 10.1038/nrg.2016.29
21. Westra HJ, Peters MJ, Esko T, Yaghootkar H, Schurmann C, Kettunen J, et al. Systematic identification of trans eQTLs as putative drivers of known disease associations. Nat Genet (2013) 45(10):1238–43. doi: 10.1038/ng.2756
22. Lloyd-Jones LR, Holloway A, McRae A, Yang J, Small K, Zhao J, et al. The genetic architecture of gene expression in peripheral blood. Am J Hum Genet (2017) 100(2):228–37. doi: 10.1016/j.ajhg.2016.12.008
23. Consortium G. The GTEx consortium atlas of genetic regulatory effects across human tissues. SCIENCE (2020) 369(6509):1318–30. doi: 10.1126/science.aaz1776
24. Gallagher MD, Chen-Plotkin AS. The post-GWAS era: from association to function. Am J Hum Genet (2018) 102(5):717–30. doi: 10.1016/j.ajhg.2018.04.002
25. Swaminathan S, Duy C, Müschen M. BACH2-BCL6 balance regulates selection at the pre-b cell receptor checkpoint. Trends Immunol (2014) 35(3):131–7. doi: 10.1016/j.it.2013.11.002
26. Igarashi K, Kurosaki T, Roychoudhuri R. BACH transcription factors in innate and adaptive immunity. Nat Rev Immunol (2017) 17(7):437–50. doi: 10.1038/nri.2017.26
27. Medici M, Porcu E, Pistis G, Teumer A, Brown SJ, Jensen RA, et al. Identification of novel genetic loci associated with thyroid peroxidase antibodies and clinical thyroid disease. PloS Genet (2014) 10(2):e1004123. doi: 10.1371/journal.pgen.1004123
28. Greulich W, Wagner M, Gaidt MM, Stafford C, Cheng Y, Linder A, et al. TLR8 is a sensor of RNase T2 degradation products. CELL (2019) 179(6):1264–1275.e13. doi: 10.1016/j.cell.2019.11.001
29. Gallo D, De Vito A, Roncoroni R, Bruno A, Piantanida E, Bartalena L, et al. A potential role of human RNASET2 overexpression in the pathogenesis of graves’ disease. ENDOCRINE (2023) 79(1):55–9. doi: 10.1007/s12020-022-03207-4
30. Liu S, Cai X, Wu J, Cong Q, Chen X, Li T, et al. Phosphorylation of innate immune adaptor proteins MAVS, STING, and TRIF induces IRF3 activation. SCIENCE (2015) 347(6227):aaa2630. doi: 10.1126/science.aaa2630
31. Bradfield JP, Qu HQ, Wang K, Zhang H, Sleiman PM, Kim CE, et al. A genome-wide meta-analysis of six type 1 diabetes cohorts identifies multiple associated loci. PloS Genet (2011) 7(9):e1002293. doi: 10.1371/journal.pgen.1002293
32. Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, et al. Genetics of rheumatoid arthritis contributes to biology and drug discovery. NATURE (2014) 506(7488):376–81. doi: 10.1038/nature12873
33. Julià A, López-Longo FJ, Pérez Venegas JJ, Bonàs-Guarch S, Olivé À, Andreu JL, et al. Genome-wide association study meta-analysis identifies five new loci for systemic lupus erythematosus. Arthritis Res Ther (2018) 20(1):100. doi: 10.1186/s13075-018-1604-1
34. Davis RS, Dennis G Jr., Odom MR, Gibson AW, Kimberly RP, Burrows PD, et al. Fc receptor homologs: newest members of a remarkably diverse fc receptor gene family. Immunol Rev (2002) 190:123–36. doi: 10.1034/j.1600-065x.2002.19009.x
35. Simmonds MJ. GWAS in autoimmune thyroid disease: redefining our understanding of pathogenesis. Nat Rev Endocrinol (2013) 9(5):277–87. doi: 10.1038/nrendo.2013.56
36. Thabet MM, Wesoly J, Slagboom PE, Toes RE, Huizinga TW. FCRL3 promoter 169 CC homozygosity is associated with susceptibility to rheumatoid arthritis in Dutch caucasians. Ann Rheum Dis (2007) 66(6):803–6. doi: 10.1136/ard.2006.064949
37. Concolino P, Mello E, Patrosso MC, Penco S, Zuppi C, Capoluongo E. p.H282N and p.Y191H: 2 novel CYP21A2 mutations in Italian congenital adrenal hyperplasia patients. Metabolism (2012) 61(4):519–24. doi: 10.1016/j.metabol.2011.08.008
38. Bánlaki Z, Doleschall M, Rajczy K, Fust G, Szilágyi A. Fine-tuned characterization of RCCX copy number variants and their relationship with extended MHC haplotypes. Genes Immun (2012) 13(7):530–5. doi: 10.1038/gene.2012.29
39. Valdes AM, Thomson G. Several loci in the HLA class III region are associated with T1D risk after adjusting for DRB1-DQB1. Diabetes Obes Metab (2009) 11(Suppl 1):46–52. doi: 10.1111/j.1463-1326.2008.01002.x
40. Fernando MM, Freudenberg J, Lee A, Morris DL, Boteva L, Rhodes B, et al. Transancestral mapping of the MHC region in systemic lupus erythematosus identifies new independent and interacting loci at MSH5, HLA-DPB1 and HLA-G. Ann Rheum Dis (2012) 71(5):777–84. doi: 10.1136/annrheumdis-2011-200808
41. Carmel M, Michaelovsky E, Weinberger R, Frisch A, Mekori-Domachevsky E, Gothelf D, et al. Differential methylation of imprinting genes and MHC locus in 22q11.2 deletion syndrome-related schizophrenia spectrum disorders. World J Biol Psychiatry (2021) 22(1):46–57. doi: 10.1080/15622975.2020.1747113
42. Zhu M, Xu K, Chen Y, Gu Y, Zhang M, Luo F, et al. Identification of novel T1D risk loci and their association with age and islet function at diagnosis in autoantibody-positive T1D individuals: based on a two-stage genome-wide association study. Diabetes Care (2019) 42(8):1414–21. doi: 10.2337/dc18-2023
43. Li Z, Akar S, Yarkan H, Lee SK, Çetin P, Can G, et al. Genome-wide association study in Turkish and Iranian populations identify rare familial Mediterranean fever gene (MEFV) polymorphisms associated with ankylosing spondylitis. PloS Genet (2019) 15(4):e1008038. doi: 10.1371/journal.pgen.1008038
44. Barrett JC, Clayton DG, Concannon P, Akolkar B, Cooper JD, Erlich HA, et al. Genome-wide association study and meta-analysis find that over 40 loci affect risk of type 1 diabetes. Nat Genet (2009) 41(6):703–7. doi: 10.1038/ng.381
45. Conigliaro P, D’Antonio A, Pinto S, Chimenti MS, Triggianese P, Rotondi M, et al. Autoimmune thyroid disorders and rheumatoid arthritis: a bidirectional interplay. Autoimmun Rev (2020) 19(6):102529. doi: 10.1016/j.autrev.2020.102529
46. Watad A, Mahroum N, Whitby A, Gertel S, Comaneshter D, Cohen AD, et al. Hypothyroidism among SLE patients: case-control study. Autoimmun Rev (2016) 15(5):484–6. doi: 10.1016/j.autrev.2016.01.019
47. Namjou B, Kelly JA, Kilpatrick J, Kaufman KM, Nath SK, Scofield RH, et al. Linkage at 5q14.3-15 in multiplex systemic lupus erythematosus pedigrees stratified by autoimmune thyroid disease. Arthritis Rheum (2005) 52(11):3646–50. doi: 10.1002/art.21413
48. Shu S, Wu HJ, Ge JY, Zeid R, Harris IS, Jovanović B, et al. Synthetic lethal and resistance interactions with BET bromodomain inhibitors in triple-negative breast cancer. Mol Cell (2020) 78(6):1096–1113.e8. doi: 10.1016/j.molcel.2020.04.027
49. Zhao R, Liu Y, Wang H, Yang J, Niu W, Fan S, et al. BRD7 plays an anti-inflammatory role during early acute inflammation by inhibiting activation of the NF-кB signaling pathway. Cell Mol Immunol (2017) 14(10):830–41. doi: 10.1038/cmi.2016.31
50. Park SW, Herrema H, Salazar M, Cakir I, Cabi S, Basibuyuk Sahin F, et al. BRD7 regulates XBP1s’ activity and glucose homeostasis through its interaction with the regulatory subunits of PI3K. Cell Metab (2014) 20(1):73–84. doi: 10.1016/j.cmet.2014.04.006
51. Marzullo P, Minocci A, Tagliaferri MA, Guzzaloni G, Di Blasio A, De Medici C, et al. Investigations of thyroid hormones and antibodies in obesity: leptin levels are associated with thyroid autoimmunity independent of bioanthropometric, hormonal, and weight-related determinants. J Clin Endocrinol Metab (2010) 95(8):3965–72. doi: 10.1210/jc.2009-2798
52. Huber A, Menconi F, Corathers S, Jacobson EM, Tomer Y. Joint genetic susceptibility to type 1 diabetes and autoimmune thyroiditis: from epidemiology to mechanisms. Endocr Rev (2008) 29(6):697–725. doi: 10.1210/er.2008-0015
53. Rosetti F, Madera-Salcedo IK, Rodríguez-Rodríguez N, Crispín JC. Regulation of activated T cell survival in rheumatic autoimmune diseases. Nat Rev Rheumatol (2022) 18(4):232–44. doi: 10.1038/s41584-021-00741-9
54. Burbano C, Villar-Vesga J, Vásquez G, Muñoz-Vahos C, Rojas M, Castaño D. Proinflammatory differentiation of macrophages through microparticles that form immune complexes leads to T- and b-cell activation in systemic autoimmune diseases. Front Immunol (2019) 10:2058. doi: 10.3389/fimmu.2019.02058
Keywords: autoimmune thyroid disease, GWAS, TWAS, FUMA, SMR
Citation: Liu X, Miao Y, Liu C, Lu W, Feng Q and Zhang Q (2023) Identification of multiple novel susceptibility genes associated with autoimmune thyroid disease. Front. Immunol. 14:1161311. doi: 10.3389/fimmu.2023.1161311
Received: 08 February 2023; Accepted: 20 April 2023;
Published: 01 May 2023.
Edited by:
Giuseppe Murdaca, University of Genoa, ItalyReviewed by:
Devis Benfaremo, Marche Polytechnic University, ItalyYoshiyuki Ban, Teikyo University Chiba Medical Center, Japan
Copyright © 2023 Liu, Miao, Liu, Lu, Feng and Zhang. 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: Qiu Zhang, zhangqiu@ahmu.edu.cn