Genetic Variants in the ST6GAL1 Gene Are Associated with Thyroglobulin Plasma Level in Healthy Individuals

Abstract

Background:

Thyroglobulin (Tg) is a 660 kDa iodoglycoprotein that serves as a scaffold for thyroid hormone synthesis. Although a twin study showed that variability of serum Tg levels has a substantial genetic basis, no genome-wide association study (GWAS) of serum/plasma Tg levels has been performed to date. The aim of this study was to identify genetic variants associated with plasma Tg levels among healthy individuals.

Methods:

A GWAS was conducted on two Croatian cohorts, and a combined analysis was performed. The analyses included 1094 individuals. A total of 7,597,379 variants, imputed using the 1000 Genomes reference panel, were analyzed for association. GWAS was performed under an additive model, controlling for age, sex, and relatedness within each data set. Combined analysis was conducted using the inverse-variance fixed-effects method.

Results:

Sixteen variants located on chromosome 3, within the ST6GAL1 gene, reached genome-wide significance. The lead SNP was rs4012172 ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 1.29 \times {10^{ - 10}}$$ \end{document} ), which explained 3.19% of the variance in Tg levels. ST6GAL1 belongs to the sialyltransferase protein family, which has a fundamental role in the synthesis of specific sialylated structures on various glycoproteins, including Tg. It is known that only immature Tg (poorly sialylated or desialylated) can be transferred to the bloodstream.

Conclusions:

A highly biologically plausible locus was identified that could have a role in the regulation of plasma Tg levels in healthy individuals.

Introduction

Thyroglobulin (Tg) is a 660 kDa iodoglycoprotein, the largest and most abundant protein in the thyroid gland (1). The synthesis of Tg occurs in the thyrocytes and neosynthesized Tg is subjected to glycosylation in the rough endoplasmic reticulum and in the Golgi compartment, followed by iodination and hormone formation at the follicular lumen/apical membrane boundary (2). The mature, hormone-containing Tg is then transferred to the follicular lumen where it functions as a thyroid hormone (TH) store (2,3). Thus, Tg can be considered as a prohormone, as it is a precursor of thyroxine (T4) and triiodothyronine (T3) (4). When THs are required, Tg is internalized by endocytosis in the thyrocyte and subjected to proteolysis, followed by the release of THs at the basolateral membrane into the bloodstream. In healthy individuals, only some Tg vesicles bypass intrathyroidal lysosomal proteolysis and reach the basolateral surface. This process, called transcytosis, seems to be the main pathway involved in the release of Tg into the circulation (5,6).

Serum Tg levels show a great deal of individual variation in healthy individuals (7). Several studies found an inverse association between iodine intake and serum Tg (8,9). Moreover, there is an association between thyroid volume and serum Tg (10), and the levels of serum Tg significantly change with thyroid abnormalities and dysfunction such as thyroiditis, autoimmune thyroiditis, and thyroid cancer (11 –13).

Glycosylation is a key event in Tg maturation (6,14). The glycan part of the Tg molecule contains 8–10% of total carbohydrates—mainly galactose, mannose, fucose, N-acetyl glucosamine, and sialic acid residues (15). Glycosylation of Tg affects several thyroid functions, including iodination and hormone synthesis, specific interactions with microsomes and membrane receptors, protein folding and recycling, as well as Tg immunoreactivity and the directing of Tg to subcellular and extracellular compartments (12).

A twin study suggested that there is a strong genetic component to the variability of serum Tg (16). No candidate or genome-wide association studies (GWAS) investigating the genes associated with Tg level have been reported so far.

The aim of this study was to identify loci associated with the level of Tg in plasma. A GWAS of plasma Tg levels was performed in 1094 healthy individuals from two Croatian cohorts.

Methods

Study population

This study was performed on subjects originating from two Croatian cohorts: one from the city of Split (CROATIA_Split) and the other from the island of Korcula (CROATIA_Korcula), within the 10,001 Dalmatians project (17). Characteristics of the cohorts are shown in Table 1. All participants who could have any thyroid disease according to anamnestic data and detailed biochemical findings were excluded. Individuals taking thyroid medication or who had undergone thyroid surgery, who reported thyroid disorders, or who had Tg, thyrotropin (TSH), free T3 (fT3), free T4 (fT4), Tg autoantibodies (TgAb), or thyroid peroxidase autoantibodies (TPOAb) levels falling outside of the normal reference range were excluded, leaving 1094 participants for analysis. Information on the family history of thyroid diseases was not available. Of note, the Croatian population has a sufficient iodine intake through legally mandated salt iodization (18). The study was approved by the Research Ethics Committees in Croatia and Scotland, and all subjects gave written informed consent.

Table 1.

Characteristics of the Two Cohorts

Variables	Split	Korcula
N	605	489
Women, n (%)	321 (53%)	297 (61%)
Median age, (q_L,q_U)	51 (39, 61)	56 (46,67)
Median BMI, kg/m² (q_L,q_U)	26.78 (23.97, 29.85)	27.52 (25.05, 30.38)
Median Tg, ng/mL (q_L,q_U)	9.20 (4.80, 14.50)	10.20 (6.40, 15.70)

BMI, body mass index; Tg, thyroglobulin.

Genotyping and imputation

Supplementary Table S1 shows cohort summary information on genotyping, imputation, and quality-control procedures. The final number of single nucleotide polymorphisms (SNPs) included in the analysis was 8,777,560 for the CROATIA_Split sample and 9,182,797 for the CROATIA_Korcula sample, with 7,597,379 overlapping SNPs present in both cohorts.

Biochemical measurements

Plasma concentrations of THs and antibodies were measured by immunoassay methods with the Liaison XL Biomedica Chemiluminescence Analyzer, using in vitro assays for quantitation of Tg (REF 311861), TSH (REF 311211), fT3 (REF 311531), fT4 (REF 311611), TgAb (REF 311711), and TPOAb (REF 311701). The following reference ranges were used for the study population: Tg 0.2–50 ng/mL, TSH 0.3–3.6 mIU/L, fT3 3.39–6.47 pmol/L, fT4 10.29–21.88 pmol/L, TgAb 5–100 IU/mL, and TPOAb levels 1–16 IU/mL (19). Biochemical measurements were performed in the Biochemistry Laboratory in the Department of Nuclear Medicine at the University Hospital Split.

Statistical analysis

A GWAS was performed within each data set, and then the results were combined.

Prior to analysis, Tg levels were adjusted for age and sex using linear regression analysis, and inverse normal transformation was used to achieve a normal distribution for the calculated residuals. GWAS was performed on transformed residuals using a linear mixed model under the assumption of additive genetic effects. An inverse normal transformation was used on Tg levels. Betas (β) represent the change in outcome for each additional effect allele under an additive model, expressed in standard deviation units/allele. To obtain more interpretable effect estimates, a linear regression analysis was also performed on the original scale of Tg levels. The association analysis for the CROATIA_Split sample was carried out using a combination of the R package “GenABEL” and “SNPTEST” software, while the R packages “GenABEL” and “VariABEL” were used for the CROATIA_Korcula sample (20 –22).

A combined analysis of the individual results was performed using a fixed-effect inverse-variance weighted model, as implemented in the R package “MetABEL” (23). To visualize the results, a Manhattan plot and a quantile–quantile (QQ) plot were generated using the R package “qqman,” regional association plots for genomic regions within 500 Kb of the top hit using Locus Zoom software, and a forest plot for the significant SNP association using “MetABEL” (24,25). Cluster plots were visually inspected for quality for all suggestively associated SNPs using the Illumina GenomeStudio software package.

To identify independently associated SNPs, conditional association analysis was performed on the lead SNP using “SNPTEST” and R software. Post hoc power calculations were performed using Quanto, assuming an additive genetic model and genome-wide significance level \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\left( {p = 5 \times {{10}^{ - 8}}} \right) ( 26 )$$ \end{document} .

Results

GWAS analysis was performed in the CROATIA_Split and CROATIA_Korcula populations separately. Four SNPs within the ST6 Beta-Galactoside Alpha-2,6-Sialyltransferase 1 (ST6GAL1) gene reached genome-wide significance in the CROATIA_Split population, and all of them were replicated in the CROATIA_Korcula population with a p-value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$< 0.0125$$ \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( 0.05 / 4$$ \end{document} ; Supplementary Table S2).

Cohort-level results were combined. There was little evidence for population stratification at a cohort-level ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${{ \rm{ \lambda }}_{{ \rm{Split}}}} = 0.99$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${{ \rm{ \lambda }}_{{ \rm{Korcula}}}} = 1.01$$ \end{document} ) or combined analysis level (λ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$= 1.004 )$$ \end{document} .

The combined analysis results are shown in Figure 1A. The QQ plot showed an excess of significant p-values (Fig. 1B). Sixteen SNPs in a region on chromosome 3 containing the ST6GAL1 gene reached genome-wide significance ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p < 5 \times {10^{ - 8}}$$ \end{document} ). The regional association plot is shown in Figure 2. The most significant SNP was rs4012172 ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 1.29 \times {10^{ - 10}}$$ \end{document} ), which explained 3.19% of the variance in Tg levels. The C allele of the rs4012172 SNP was associated with a higher Tg level (effect size [β] = 0.30, standard error [SE] = 0.05). Effect sizes were in the same direction in both cohorts (Fig. 3). The rs4012172 SNP resides within intron 2 of ST6GAL1. Of 16 identified SNPs, all variants are intronic except the rs4686837 SNP, which is a 5′ UTR variant (Supplementary Table S3). The effect sizes on the original scale of Tg for significant SNPs ranged between 1.89 and 2.45 ng/mL per associated allele (Supplementary Table S3). No secondary independent association signals with genome-wide significance were found in this region by conditional association analysis.

FIG. 1.

Manhattan plot and quantile–quantile (QQ) plot of combined analysis for thyroglobulin (Tg) levels. (A) Manhattan plot of single nucleotide polymorphisms (SNP) for Tg levels. The x-axis represents chromosomes, and the y-axis represents −log₁₀ (p-values) obtained by combined analysis. The red line indicates the genome-wide threshold for significance \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\left( {p = 5 \times {{10}^{ - 8}}} \right)$$ \end{document} , while the blue line indicates the suggestive threshold of significance \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\left( {p = 5 \times {{10}^{ - 6}}} \right)$$ \end{document} . (B) The QQ plot shows the divergence between the observed and expected −log₁₀ (p-value) of the combined analysis results (black dots). The red line indicates the expected distributions under the null hypothesis. The gray lines represent the 95% confidence interval. Color images are available online.

FIG. 2.

Regional association plot of ST6GAL1 region. The most significant SNP (rs4012172) is shown in purple. Otherwise, the colors of the circles denote their correlations (LD r ²) with the top SNP. Color images are available online.

FIG. 3.

Forest plot of the association between the ST6GAL1 gene polymorphism rs4012172 and Tg levels.

In addition to the genome-wide significant association identified for SNPs within ST6GAL1, several suggestive associations were also identified with a p-value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$< 5 \times {10^{ - 6}}$$ \end{document} (Table 2). These are rs35464485 near the OR5B12 gene ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 8.01 \times {10^{ - 7}}$$ \end{document} ), rs140183574 within the GABRB3 and GABRA5 genes ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 1.90 \times {10^{ - 6}}$$ \end{document} ), rs34017077 near the KAT14 gene ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 2.85 \times {10^{ - 6}}$$ \end{document} ), rs2966248 within the HSD17B2 gene ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 2.89 \times {10^{ - 6}}$$ \end{document} ), rs79641712 within the GADD45G gene ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 2.89 \times {10^{ - 6}}$$ \end{document} ), rs9842143 within the CCDC54 gene ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 3.53 \times {10^{ - 6}}$$ \end{document} ), and rs117636523 within the POM121L12 gene ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p = 3.98 \times {10^{ - 6}}$$ \end{document} ). Association results for all suggestively associated SNPs are shown in Supplementary Table S3.

Table 2.

Associations of Top Single Nucleotide Polymorphisms \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\left( {p < 5 \times {{10}^{ - 6}}} \right)$$ \end{document} ) with Thyroglobulin Levels

SNP	Genome context	Chr	Position	Nearest gene	EA/OA	EAF	β	SE	p-Value
rs4012172	Intron	3	186741511	ST6GAL1	C/T	0.65	0.30	0.05	1.29 × 10⁻¹⁰
rs35464485	Intergenic	11	58224123	OR5B12	T/A	0.97	0.61	0.12	8.01 × 10⁻⁷
rs140183574	Intron	15	27137583	GABRB3, GABRA5	C/G	0.83	0.41	0.09	1.90 × 10⁻⁶
rs34017077	Intergenic	20	18177700	KAT14	G/A	0.76	0.30	0.06	2.85 × 10⁻⁶
rs2966248	3′ UTR variant	16	82126916	HSD17B2	C/A	0.62	0.21	0.05	2.89 × 10⁻⁶
rs79641712	Intron	9	92305846	GADD45G	T/C	0.05	0.48	0.10	3.19 × 10⁻⁶
rs9842143	Intron	3	107052083	CCDC54	A/G	0.51	0.20	0.04	3.53 × 10⁻⁶
rs117636523	Intergenic	7	52787062	POM121L12	G/A	0.27	0.23	0.05	3.98 × 10⁻⁶

SNP, single nucleotide polymorphism.

Discussion

In this GWAS of plasma Tg levels in healthy populations, the ST6GAL1 gene on chromosome 3 was identified as a novel significant locus. The strongest association was observed for the ST6GAL1 gene polymorphism rs4012172 (p = 1.29 × 10⁻¹⁰), with the C allele being associated with higher Tg levels. The identified SNP accounts for 3.19% of population variance in plasma Tg. ST6GAL1 belongs to the sialyltransferase family, Golgi membrane-bound enzymes that catalyze the biosynthesis of sialylated oligosaccharides, and therefore has a fundamental role in the synthesis of specific sialylated structures (27). The highest expression of ST6GAL1 is found in liver, thyroid, spleen, lymph node, and prostate tissue (28). According to the Human Splicing Finder (29), there is no information on the predictive effect of the rs4012172 SNP in ST6GAL1 protein expression or activity. Among the remaining 15 significant polymorphisms, the rs3872724 variant was predicted to create a donor splice site within the intron, while for other polymorphisms, there was no available prediction information. As specified in GTEx (30), the strongest eQTL association signals reported in thyroid tissue for this gene are for rs3821819 and rs967367 (p = 8.2 × 10⁻¹³). Both variants are also significantly associated with the level of Tg in this data set (p = 1.9 × 10⁻⁹ and p = 2.42 × 10⁻⁹, respectively), and are in moderate LD (r ² = 0.74, 1000 Genomes Phase III) with the lead rs4012172 SNP.

Since all of the GWAS hits are located in non-coding sites, either untyped rare variants or more indirect mechanisms could be responsible for the variation of plasma Tg levels.

Besides iodine availability, sialylation of Tg plays an important role in TH synthesis (2). Sialylation of Tg is associated with three aspects of Tg biology: recycling, autoregulation, and immunoreactivity (12).

Tg stored in the follicular lumen is internalized in thyrocytes and broken down to release THs. A portion of the internalized Tg (with no or low hormone content) that is poorly sialylated and iodinated cannot be proteolyzed by lysosomes (6). The thyroid gland rather tries to protect and recycle the immature Tg through a two-step transport system: first, through the Golgi compartment where novel oligosaccharides are added, followed by transport back to the follicular lumen where additional iodination and hormone formation occurs (2,6). Another fraction of the immature Tg molecules is transported via transcytotic vesicles to the basolateral thyrocyte membrane and then released into the bloodstream (2,31).

One of the potential mechanisms of Tg autoregulation and endocytosis is placed at the apical membrane of thyrocytes where there is direct contact of signal transduction receptors with Tg from the follicular lumen. Up to now, three candidate signal transduction receptors have been described: the asialoglycoprotein (ASGP), N-acetylglucosamine receptor, and megalin. Sialylation of Tg is considered necessary for the release of Tg into thyrocytes via receptor-mediated endocytosis. The ASGP and N-acetylglucosamine receptors are probably involved in sorting immature Tg molecules for recycling to the follicular lumen. Poorly sialylated Tg located in the periphery of the follicular lumen has a greater affinity for the ASGP receptor than normally sialylated Tg (2,32). Sialylation could also have an important role in the variability of plasma Tg levels in healthy individuals, since it seems to play an important role in regulating Tg secretion into the bloodstream (33).

Sialic acids of native Tg have a protective role in the autoimmunization process. Several studies showed that desialylation may unmask the Tg isoantigen (34,35). There is now comprehensive evidence that the trigger for an autoimmune reaction is increased synthesis of desialylated Tg followed by the initiation of the autoimmunity process (35).

The polypeptide chain of the human Tg molecule contains 20 putative N-linked glycosylation sites, of which 16 are glycosylated in the mature Tg molecule (36,37). Sialylation is a late post-translational modification of human Tg. It has been reported that normal human Tg has a sialic acid content of 7.3 μg/mg (38). Sialyltransferases catalyze the transfer of sialic acid from cytidine monophosphosialic acid to the terminal position of carbohydrate chains of glycolipids and glycoproteins (39). ST6GAL1 encodes sialyltransferase 1, also known as β-galactiside α-2,6-sialyltransferase (39). Sialic acids are attached to the carbohydrate units of human Tg primarily with α-2,6 linkages established by the activity of sialyltransferase 1 (38). Alpha 2,6-bound sialic acids were also described on the Tg secreted by rat thyroid cell line FRTL-5 and porcine Tg (40,41). Furthermore, human thyroid tissue showed high mRNA levels of ST6GAL1 (39). Severely hyposialylated Tg was noticed in a patient with congenital goiter with hypothyroidism, and the authors observed that the α-2,6 sialyltransferase activity was deficient in that family (38). Therefore, based on this evidence and the present finding, it is suggested that ST6GAL1 variations may alter Tg sialylation and consequently influence the release of Tg into the bloodstream.

Genetic defects that affect the biochemical pathways of glycan synthesis can cause severe glycan deficiencies and lead to a variety of diseases (42). Several genome-wide studies have already revealed significant associations with variants in the ST6GAL1 gene and different traits and diseases. So far, the most significant association discovered with a variant in ST6GAL1 was identified for plasma immunoglobulin G (IgG) glycosylation, with rs11710456 as the leading SNP (p = 6.12 × 10⁻⁷⁵). Lauc et al. discovered 37 genome-wide significant SNPs within ST6GAL1 that were associated with 14 different IgG glycosylation/sialylation traits in a meta-analysis of European populations (43). Highly sialylated IgG molecules are considered as one of the key elements in immune homeostasis, as well as in the prevention and regulation of autoimmune and inflammatory diseases (44).

In a Han Chinese population, the ST6GAL1 gene was associated with immunoglobulin A nephropathy (rs7634389; p = 7 × 10⁻¹⁰) (45). Drug-induced liver injury was also associated with the ST6GAL1 gene (rs10937275; p = 1.4 × 10⁻⁸) in Europeans (46). These GWAS indicate that ST6GAL1 is a plausible cofactor in the susceptibility to immuno-inflammatory events.

Kooner et al. found an association of the ST6GAL1 gene (rs16861329, p = 3 × 10⁻⁸) with type 2 diabetes in a South Asian population (47). A study by Santisteban et al. showed that insulin and insulin growth factor 1, along with TSH, increase Tg synthesis at the mRNA levels in FRTL-5 thyroid cells (48). Changes in rat hepatocyte cell surface sialic acid content have also been observed in insulin resistance associated with diabetes mellitus (49).

Sialylation regulates immune cell function, cell signaling, activation, and differentiation (50 –53), thus supporting the findings of GWAS on the pleiotropy of ST6GAL1. The results of the present study suggest that ST6GAL1 is also likely to have a role in the regulation of the plasma Tg level.

Although no signals other than the ST6GAL1 locus reached genome-wide significance, several candidate loci showed suggestive evidence of association. Of particular interest is the rs140183574 variant located on chromosome 15 (p = 1.9 × 10⁻⁶), within the GABRB3 and GABRA5 genes. GABRB3 and GABRA5 encode proteins that serve as the receptors for gamma-aminobutyric acid (GABA), a major inhibitory neurotransmitter of the nervous system. GABA plays a role in (43) the homeostasis of the normal thyroid and has an inhibitory role in TSH-stimulated TH secretion from thyroid follicular cells (54,55).

Another interesting finding is the suggestive association for a locus near the POM121L12 gene on chromosome 7 (p = 4 × 10⁻⁶). Variants near this gene were also suggestively associated with N-glycosylation of human IgG (p = 5 × 10⁻⁶) (43). Likewise, several autoimmune diseases (type 1 diabetes, type 2 diabetes, psoriasis, and rheumatoid arthritis) have shown suggestive associations (p = 10⁻⁵–10⁻⁷) within the same region (56 –58), indicating that the POM121L12 gene could also have a pleiotropic effect.

A major strength of this study is the careful selection of healthy individuals based on detailed thyroid function assessment (including measurements of Tg, TSH, fT3, fT4, TgAb, and TPOAb) and information on self-reported thyroid disease and thyroid medications used. Other strengths include the ethnically homogeneous sample of Croatian individuals and the comprehensive set of genetic variants examined. Although both cohorts consist predominantly of older subjects, all participants included in the study were healthy individuals with normal thyroid function. Additionally, since thyroid function decreases with age (59), Tg levels were adjusted for age effects in the analysis. The main limitation of the study is the modest sample size. Thus, the study was underpowered to detect small effects at genome-wide significance. However, there were sufficient data to detect an association for the identified ST6GAL1 locus confidently, since the combined analysis had a 91% power to detect the top associated SNP rs4012172 with an effect size of 0.30 and minor allele frequency of 0.35 at the genome-wide level of significance.

In conclusion, this study reports the first GWAS on plasma Tg levels and contributes to the understanding of genetic factors underlying Tg levels in healthy individuals. The novel genome-wide significant ST6GAL1 locus revealed by this study is highly biologically plausible and likely to have a role in the regulation of plasma Tg levels. ST6GAL1 encodes sialyltransferase 1 involved in the sialylation of Tg (37,39). It is known that only desialylated or poorly sialylated Tg can be transferred into the bloodstream. The biological role of sialyltransferase 1 suggests a possible influence of ST6GAL1 variations on plasma Tg levels.

Footnotes

Acknowledgments

The Croatian Science Foundation funded this work under the project “Identification of new genetic loci implicated in regulation of thyroid and parathyroid function” (grant no. 1498). The 10,001 Dalmatians project was funded by grants from the Medical Research Council (United Kingdom), European Commission Framework 6 project EUROSPAN (contract no. LSHG-CT-2006-018947), the Republic of Croatia Ministry of Science, Education and Sports research grant (216-1080315-0302), the Croatian Science Foundation (grant 8875), CEKOM (Ministry of Economy, Entrepreneurship and Crafts), and the Research Centre of Excellence in Personalized Medicine (Ministry of Science and Education). C.H. is supported by an MRC University Unit Program Grant MC_PC_U127592696.

We would like to thank all participants of this study and acknowledge the invaluable support of the local teams in Zagreb and Split, especially that of the Institute for Anthropological Research, Zagreb, Croatia. We wish to thank Shelly Pranic, PhD, for proofreading the manuscript.

Author Disclosure Statement

No competing financial interests exist.

Supplementary Material

Supplementary Table S1

Supplementary Table S2

Supplementary Table S3

References

Malthiery

, Marriq

, Berge-Lefranc

, Franc

, Henry

, Lejeune

, Ruf

, Lissitzky

. 1989. Thyroglobulin structure and function: recent advances. Biochimie, 71:195–209.

Rousset

, Dupuy

, Miot

, Dumont

. 2000. Thyroid hormone synthesis and secretion. In: De Groot

, Chrousos

, Dungan

, Feingold

, Grossman

, Hershman

, Koch

, Korbonits

, McLachlan

, New

, Purnell

, Rebar

, Singer

, Vinik

, (eds) Endotext. Vol. MDText.com, Inc., South Dartmouth, MA.

Dunn

, Dunn

. 1999. The importance of thyroglobulin structure for thyroid hormone biosynthesis. Biochimie, 81:505–509.

Yang

, Pollock

, Rawitch

. 1996. Glycosylation in human thyroglobulin: location of the N-linked oligosaccharide units and comparison with bovine thyroglobulin. Arch Biochem Biophys, 327:61–70.

Marino

, Chiovato

, Mitsiades

, Latrofa

, Andrews

, Tseleni-Balafouta

, Collins

, Pinchera

, McCluskey

. 2000. Circulating thyroglobulin transcytosed by thyroid cells in complexed with secretory components of its endocytic receptor megalin. J Clin Endocrinol Metab, 85:3458–3467.

Hatipoglu

, Schneider

. 1999. Selective endocytosis of thyroglobulin: a review of potential mechanisms for protecting newly synthesized molecules from premature degradation. Biochimie, 81:549–555.

Feldt-Rasmussen

, Hyltoft Petersen

, Blaabjerg

, Horder

. 1980. Long-term variability in serum thyroglobulin and thyroid related hormones in healthy subjects. Acta Endocrinol (Copenh), 95:328–334.

Knudsen

, Bulow

, Jorgensen

, Perrild

, Ovesen

, Laurberg

. 2001. Serum Tg—a sensitive marker of thyroid abnormalities and iodine deficiency in epidemiological studies. J Clin Endocrinol Metab, 86:3599–3603.

Vejbjerg

, Knudsen

, Perrild

, Laurberg

, Carle

, Pedersen

, Rasmussen

, Ovesen

, Jorgensen

. 2009. Thyroglobulin as a marker of iodine nutrition status in the general population. Eur J Endocrinol, 161:475–481.

10.

Cahoon

, Rozhko

, Hatch

, Polyanskaya

, Ostroumova

, Tang

, Nadirov

, Yauseyenka

, Savasteeva

, McConnell

, Pfeiffer

, Brenner

. 2013. Factors associated with serum thyroglobulin levels in a population living in Belarus. Clin Endocrinol, 79:120–127.

11.

Indrasena

BSH

. 2017. Use of thyroglobulin as a tumour marker. World J Biol Chem, 8:81–85.

12.

Xavier

, Maciel

, Vieira

, Dias-da-Silva

, Martins

. 2016. Insights into the posttranslational structural heterogeneity of thyroglobulin and its role in the development, diagnosis, and management of benign and malignant thyroid diseases. Arch Endocrinol Metab, 60:66–75.

13.

Malandrino

, Latina

, Marescalco

, Spadaro

, Regalbuto

, Fulco

, Scollo

, Vigneri

, Pellegriti

. 2011. Risk-adapted management of differentiated thyroid cancer assessed by a sensitive measurement of basal serum thyroglobulin. J Clin Endocrinol Metab, 96:1703–1709.

14.

Arima

, Spiro

. 1972. Studies on the carbohydrate units of thyroglobulin. Evaluation of their microheterogeneity in the human and calf proteins. J Biol Chem, 247:1825–1835.

15.

Venkatesh

, Deshpande

. 1999. A comparative review of the structure and biosynthesis of thyroglobulin. Comp Biochem Physiol C Pharmacol Toxicol Endocrinol, 122:13–20.

16.

Premawardhana

, Lo

, Phillips

, Prentice

, Rees Smith

. 1994. Variability of serum thyroglobulin levels is determined by a major gene. Clin Endocrinol (Oxf), 41:725–729.

17.

Rudan

, Marusic

, Jankovic

, Rotim

, Boban

, Lauc

, Grkovic

, Dogas

, Zemunik

, Vatavuk

, Bencic

, Rudan

, Mulic

, Krzelj

, Terzic

, Stojanovic

, Puntaric

, Bilic

, Ropac

, Vorko-Jovic

, Znaor

, Stevanovic

, Biloglav

, Polasek

. 2009. “10 001 Dalmatians”: Croatia launches its national biobank. Croat Med J, 50:4–6.

18.

Kusic

, Jukic

, Rogan

, Juresa

, Dabelic

, Stanicic

, Boric

, Lukinac

, Mihaljevic

, Punda

, Smokvina

, Topalovic

, Katalenic

. 2012. Current status of iodine intake in Croatia—the results of 2009 survey. Coll Antropol, 36:123–128.

19.

Baric

, Brcic

, Gracan

, Skrabic

, Brekalo

, Simunac

, Lovric

, Anic

, Barbalic

, Zemunik

, Punda

, Boraska Perica

. 2018. Thyroglobulin antibodies are associated with symptom burden in patients with Hashimoto's thyroiditis: a cross-sectional study. Immunol Invest, 48:198–209.

20.

Aulchenko

, Ripke

, Isaacs

, Van Duijn

. 2007. GenABEL: an R library for genome-wide association analysis. Bioinformatics, 23:1294–1296.

21.

Marchini

, Howie

, Myers

, McVean

, Donnelly

. 2007. A new multipoint method for genome-wide association studies by imputation of genotypes. Nat Genet, 39:906–913.

22.

Struchalin

, Amin

, Eilers

PHC

, van Duijn

, Aulchenko

. 2012. An R package “VariABEL” for genome-wide searching of potentially interacting loci by testing genotypic variance heterogeneity. BMC Genet, 13:4.

23.

R Core Team 2004 R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.

24.

Pruim

, Welch

, Sanna

, Teslovich

, Chines

, Gliedt

, Boehnke

, Abecasis

, Willer

. 2010. LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics, 26:2336–2337.

25.

Turner

. 2014 qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. biorXiv DOI: 10.1101/005165.

26.

Gauderman WJ

. 2006. QUANTO 1.1: A Computer Program for Power and Sample Size Calculations for Genetic–Epidemiology Studies. University of Southern California, Los Angeles, California.

27.

Harduin-Lepers

, Vallejo-Ruiz

, Krzewinski-Recchi

M-A

, Samyn-Petit

, Julien

, Delannoy

. 2001. The human sialyltransferase family. Biochimie, 83:727–737.

28.

The Human Protein Atlas. ST6GAL1. Available at: https://www.proteinatlas.org/ENSG00000073849-ST6GAL1/tissue (accessed October 26, 2017 ).

29.

Desmet

, Hamroun

, Lalande

, Collod-Beroud

, Claustres

, Beroud

. 2009. Human Splicing Finder: an online bioinformatics tool to predict splicing signals. Nucleic Acids Res, 37:e67.

30.

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

31.

Ohtaki

, Moriya

, Suzuki

, Moriuchi

. 1967. Nonhormonal iodine escape from the normal and abnormal thyroid gland. J Clin Endocrinol Metab, 27:728–740.

32.

Sellitti

, Suzuki

. 2014. Intrinsic regulation of thyroid function by thyroglobulin. Thyroid, 24:625–638.

33.

Eggo

, Burrow

. 1983. Glycosylation of thyroglobulin—its role in secretion, iodination, and stability. Endocrinology, 113:1655–1663.

34.

Fenouillet

, Fayet

, Hovsepian

, Bahraoui

, Ronin

. 1986. Immunochemical evidence for a role of complex carbohydrate chains in thyroglobulin antigenicity. J Biol Chem, 261:15153–15158.

35.

Salabe

, Dominici

, Salabe

. 1976. Immunological properties of Tg carbohydrates: enhancement of Tg immunoreaction by removal of sialic acid. Clin Exp Immunol, 25:234–243.

36.

Yang

, Pollock

, Rawitch

. 1996. Glycosylation in human thyroglobulin: location of the N-linked oligosaccharide units and comparison with bovine thyroglobulin. Arch Biochem Biophys, 327:61–70.

37.

Zabczynska

, Kozlowska

, Pochec

. 2018. Glycosylation in the thyroid gland: vital aspects of glycoprotein function in thyrocyte physiology and thyroid disorders. Int J Mol Sci, 19.

38.

Grollman

, Doi

, Weiss

, Ashwell

, Wajchenberg

, Medeiros-Neto

. 1992. Hyposialylated thyroglobulin in a patient with congenital goiter and hypothyroidism. J Clin Endocrinol Metab, 74:43–48.

39.

Kiljanski

, Ambroziak

, Pachucki

, Jazdzewski

, Wiechno

, Stachlewska

, Gornicka

, Bogdanska

, Nauman

, Bartoszewicz

. 2005. Thyroid sialyltransferase mRNA level and activity are increased in Graves' disease. Thyroid, 15:645–652.

40.

Grollman

, Saji

, Shimura

, Lau

JTY

, Ashwell

. 1993. Thyrotropin regulation of sialic-acid expression in rat-thyroid cells. J Biol Chem, 268:3604–3609.

41.

Patwa

, Zhao

, Anderson

, Simeone

, Lubman

. 2006. Screening of glycosylation patterns in serum using natural glycoprotein microarrays and multi-lectin fluorescence detection. Anal Chem, 78:6411–6421.

42.

Freeze

. 2006. Genetic defects in the human glycome. Nat Rev Genet, 7:537–551.

43.

Lauc

, Huffman

, Pucic

, Zgaga

, Adamczyk

, Muzinic

, Novokmet

, Polasek

, Gornik

, Kristic

, Keser

, Vitart

, Scheijen

, Uh

, Molokhia

, Patrick

, McKeigue

, Kolcic

, Lukic

, Swann

, van Leeuwen

, Ruhaak

, Houwing-Duistermaat

, Slagboom

, Beekman

, de Craen

, Deelder

, Zeng

, Wang

, Hastie

, Gyllensten

, Wilson

, Wuhrer

, Wright

, Rudd

, Hayward

, Aulchenko

, Campbell

, Rudan

. 2013. Loci associated with N-glycosylation of human immunoglobulin G show pleiotropy with autoimmune diseases and haematological cancers. PLoS Genet, 9:e1003225.

44.

Lauc

, Pezer

, Rudan

, Campbell

. 2016. Mechanisms of disease: the human N-glycome. Biochim Biophys Acta, 1860:1574–1582.

45.

, Foo

, Wang

, Low

, Tang

, Toh

, Yin

, Khor

, Goh

, Irwan

, Xu

, Andiappan

, Bei

, Rotzschke

, Chen

, Cheng

, Sun

, Jiang

, Wong

, Lin

, Aung

, Liao

, Saw

, Ye

, Ebstein

, Chen

, Shi

, Chew

, Chen

, Zhang

, Li

, Xu

, Tai

, Wang

, Chen

, Zhang

, Zeng

, Zhang

, Liu

, Yu

, Liu

. 2015. Identification of new susceptibility loci for IgA nephropathy in Han Chinese. Nat Commun, 6:7270.

46.

Daly

, Donaldson

, Bhatnagar

, Shen

, Pe'er

, Floratos

, Daly

, Goldstein

, John

, Nelson

, Graham

, Park

, Dillon

, Bernal

, Cordell

, Pirmohamed

, Aithal

, Day

, Study

, International

SAEC

. 2009. HLA-B*5701 genotype is a major determinant of drug-induced liver injury due to flucloxacillin. Nat Genet, 41:816–819.

47.

Kooner

, Saleheen

, Sim

, Sehmi

, Zhang

, Frossard

, Been

, Chia

, Dimas

, Hassanali

, Jafar

, Jowett

, Li

, Radha

, Rees

, Takeuchi

, Young

, Aung

, Basit

, Chidambaram

, Das

, Grundberg

, Hedman

, Hydrie

, Islam

, Khor

, Kowlessur

, Kristensen

, Liju

, Lim

, Matthews

, Liu

, Morris

, Nica

, Pinidiyapathirage

, Prokopenko

, Rasheed

, Samuel

, Shah

, Shera

, Small

, Suo

, Wickremasinghe

, Wong

, Yang

, Zhang

, Diagram

MuTher

, Abecasis

, Barnett

, Caulfield

, Deloukas

, Frayling

, Froguel

, Kato

, Katulanda

, Kelly

, Liang

, Mohan

, Sanghera

, Scott

, Seielstad

, Zimmet

, Elliott

, Teo

, McCarthy

, Danesh

, Tai

, Chambers

. 2011. Genome-wide association study in individuals of South Asian ancestry identifies six new type 2 diabetes susceptibility loci. Nat Genet, 43:984–989.

48.

Santisteban

, Kohn

, Di Lauro

. 1987. Thyroglobulin gene expression is regulated by insulin and insulin-like growth factor I, as well as thyrotropin, in FRTL-5 thyroid cells. J Biol Chem, 262:4048–4052.

49.

Salhanick

, Amatruda

. 1988. Role of sialic acid in insulin action and the insulin resistance of diabetes mellitus. Am J Physiol, 255:E173–179.

50.

Jenner

, Kerst

, Handgretinger

, Muller

. 2006. Increased alpha2,6-sialylation of surface proteins on tolerogenic, immature dendritic cells and regulatory T cells. Exp Hematol, 34:1212–1218.

51.

Woodard-Grice

, McBrayer

, Wakefield

, Zhuo

, Bellis

. 2008. Proteolytic shedding of ST6Gal-I by BACE1 regulates the glycosylation and function of alpha4beta1 integrins. J Biol Chem, 283:26364–26373.

52.

Brennan

, Saouaf

, Van Dyken

, Marth

, Li

, Bhandoola

, Greene

. 2006. Sialylation regulates peripheral tolerance in CD4+ T cells. Int Immunol, 18:627–635.

53.

Daniels

, Hogquist

, Jameson

. 2002. Sweet ‘n’ sour: the impact of differential glycosylation on T cell responses. Nat Immunol, 3:903–910.

54.

Ahren

. 1989. GABA inhibits thyroid hormone secretion in the mouse. Thyroidology, 1:105–108.

55.

Wiens

, Trudeau

. 2006. Thyroid hormone and gamma-aminobutyric acid (GABA) interactions in neuroendocrine systems. Comp Biochem Physiol A Mol Integr Physiol, 144:332–344.

56.

Wellcome Trust Case Control Consortium. 2007. Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature, 447:661–678.

57.

Okada

, Terao

, Ikari

, Kochi

, Ohmura

, Suzuki

, Kawaguchi

, Stahl

, Kurreeman

, Nishida

, Ohmiya

, Myouzen

, Takahashi

, Sawada

, Nishioka

, Yukioka

, Matsubara

, Wakitani

, Teshima

, Tohma

, Takasugi

, Shimada

, Murasawa

, Honjo

, Matsuo

, Tanaka

, Tajima

, Suzuki

, Iwamoto

, Kawamura

, Tanii

, Okazaki

, Sasaki

, Gregersen

, Padyukov

, Worthington

, Siminovitch

, Lathrop

, Taniguchi

, Takahashi

, Tokunaga

, Kubo

, Nakamura

, Kamatani

, Mimori

, Plenge

, Yamanaka

, Momohara

, Yamada

, Matsuda

, Yamamoto

. 2012. Meta-analysis identifies nine new loci associated with rheumatoid arthritis in the Japanese population. Nat Genet, 44:511–516.

58.

Capon

, Bijlmakers

, Wolf

, Quaranta

, Huffmeier

, Allen

, Timms

, Abkevich

, Gutin

, Smith

, Warren

, Young

, Worthington

, Burden

, Griffiths

, Hayday

, Nestle

, Reis

, Lanchbury

, Barker

, Trembath

. 2008. Identification of ZNF313/RNF114 as a novel psoriasis susceptibility gene. Hum Mol Genet, 17:1938–1945.

59.

Vadiveloo

, Donnan

, Murphy

, Leese

. 2013. Age- and gender-specific TSH reference intervals in people with no obvious thyroid disease in Tayside, Scotland: The Thyroid Epidemiology, Audit, and Research Study (TEARS). J Clin Endocrinol Metab, 98:1147–1153.