Multivariate genome-wide analysis of immunoglobulin G N-glycosylation identifies new loci pleiotropic with immune function

Abstract

The N-glycosylation of immunoglobulin G (IgG) affects its structure and function. It has been demonstrated that IgG N-glycosylation patterns are inherited as complex quantitative traits. Genome-wide association studies identified loci harboring genes encoding enzymes directly involved in protein glycosylation as well as loci likely to be involved in regulation of glycosylation biochemical pathways. Many of these loci could be linked to immune functions and risk of inflammatory and autoimmune diseases. The aim of the present study was to discover and replicate new loci associated with IgG N-glycosylation and to investigate possible pleiotropic effects of these loci onto immune function and the risk of inflammatory and autoimmune diseases. We conducted a multivariate genome-wide association analysis of 23 IgG N-glycosylation traits measured in 8090 individuals of European ancestry. The discovery stage was followed up by replication in 3147 people and in silico functional analysis. Our study increased the total number of replicated loci from 22 to 29. For the discovered loci, we suggest a number of genes potentially involved in the control of IgG N-glycosylation. Among the new loci, two (near RNF168 and TNFRSF13B) were previously implicated in rare immune deficiencies and were associated with levels of circulating immunoglobulins. For one new locus (near AP5B1/OVOL1), we demonstrated a potential pleiotropic effect on the risk of asthma. Our findings underline an important link between IgG N-glycosylation and immune function and provide new clues to understanding their interplay.

Introduction

Immunoglobulin G (IgG) is the most abundant immunoglobulin isotype and the most prevalent glycoprotein in human plasma (1,2). IgG has a broad spectrum of activity and participates in multiple immune reactions including phagocytosis, complement activation and antibody-dependent cell-mediated cytotoxicity (1,3,4). All IgG molecules are N-glycosylated at a highly conserved asparagine 297 residue (also referred to as CH2-84.4) located in the Fc region, which interacts with effector cells of the immune system and some proteins of the complement system (5). Additionally, glycosylation of variable regions of light or heavy chains can occur in approximately 20% of IgG molecules (6). The attachment of different glycans generates a wide range of IgG glycoforms. The spectrum of glycans linked to IgG molecules is called the IgG N-glycome.

N-glycosylation of IgG contributes to the maintenance of its structure, stability and solubility. No less importantly, N-glycosylation drives IgG’s effector functions, thus playing an important regulatory role in the immune system (4,5,7). For example, addition of fucose to the IgG glycan core structure decreases the affinity of IgG to FcγRIIIA and FcγRIIIB receptors and diminishes antibody-dependent cell-mediated cytotoxicity (ADCC). On the other hand, increased levels of IgG sialylation and terminal galactosylation are generally considered to have an anti-inflammatory effect (reviewed in (4)). IgG N-glycome composition is affected in various physiological and pathological states and is associated with a number of diseases, such as inflammatory, infectious, autoimmune and alloimmune diseases, as well as different types of cancer (4). IgG glycosylation patterns are also considered as biomarkers of ageing (8).

Although biochemical reactions underlying protein glycosylation are well studied, the biological networks regulating these reactions remain poorly understood. Elucidating the genetic factors influencing IgG N-glycome variation helps to uncover molecules and pathways involved in establishing IgG glycosylation patterns and facilitates research linking IgG N-glycome and human complex diseases. To date, four genome-wide association studies (GWASes) of the total plasma IgG N-glycome have been published (9–12). The largest study was conducted by Klarić et al. (12), which identified and replicated 19 loci and brought the number of loci robustly associated with the IgG N-glycome composition to 22. This work suggested how these loci form a functional network and proposed new regulators of IgG glycosylation.

In addition to the loci that contain glycosyltransferase genes (MGAT3, ST6GAL1, B4GALT1, FUT6/FUT3, FUT8), previous analyses identified loci that contain genes encoding transcription factors known to play central roles in lymphocyte maturation and differentiation: RUNX1, RUNX3, IKZF1 and IKZF3 (9,11,12). For several IgG N-glycosylation loci, colocalization analysis suggested pleiotropic effects onto autoimmune and inflammatory conditions (12). Notably, for the most highly replicated and significant asthma locus at 17q12-21 (near ORMDL3-GSDMB-IKZF3 genes) (13), the asthma risk alleles were associated with the decreased levels of afucosylated, agalactosylated and monogalactosylated glycans and increased levels of fucosylated glycans (12). This locus was also associated with the risk of autoimmune and inflammatory diseases such as rheumatoid arthritis, Crohn’s disease, ulcerative colitis and primary biliary cirrhosis, and the autoimmune disease risk allele was also associated with increased levels of afucosylated agalactosylated and monogalactosylated glycans and decreased levels of fucosylated glycans (12). Consistent with these observations, it is known that although the same single nucleotide polymorphisms (SNPs) are associated with both asthma and autoimmune disease, the risk alleles associated with autoimmune disease are opposite of those associated with asthma risk (13). Similar to observations in the ORMDL3-GSDMB-IKZF3 locus, at another locus near the IRF1-SLC22A4 genes, the Crohn’s disease risk allele was demonstrated to be associated with the increased levels of afucosylated agalactosylated glycans and decreased levels of fucosylation (12).

The above observations, in general, suggest that a more proinflammatory and pro-ADCC glycan profile is associated with increased risk of autoimmune and inflammatory diseases. However, the evidence is still scattered, and identification of additional IgG N-glycosylation loci that overlap with the loci affecting the disease risk may help in understanding the relationship between glycosylation and immune function. Application of multivariate methods in genetic association studies of correlated phenotypes can increase the statistical power to discover novel loci (10,14,15) as well as facilitate the interpretation of the results (16). The efficacy of a multivariate approach in investigating the genetic architecture of IgG N-glycosylation traits has recently been demonstrated by Shen et al. (10). At the time when only nine loci were shown to be associated with IgG N-glycosylation (9), of which two were replicated, the study employing multivariate methodology identified and replicated five novel loci despite using a smaller discovery sample size.

In the present study, we aimed to provide further insights into the genetics underpinning IgG N-glycosylation by applying a multivariate approach to the largest currently available dataset linking genotypes with IgG glycan measurements [N = 8090; Klarić et al.’s study (12)] and replicating our results using data from 3147 independent samples.

Results

Phenotype grouping

Following the study by Shen et al. (10), we grouped 23 directly measured IgG N-glycan traits (described in Supplementary Material, Table S1) into nine groups based on the general structural and chemical properties of glycans, including one group that combined all 23 IgG N-glycosylation traits. The description of these groups is provided in Figure 1.

Figure 1

Directly measured IgG N-glycan traits and trait grouping for multivariate analysis. IGP1-23 traits (Supplementary Material, Table S1) correspond to UPLC peaks and approximate the percentage of a certain glycan (the most abundant glycan in the peak) in total IgG glycans (18, 19). Sugar residues are graphically represented according to Essentials of Glycobiology, 2nd edition (20). Glycan names are given according to the Oxford nomenclature (described in (21)). GlcNAc, N-Acetylglucosamine.

Discovery and replication of loci associated with IgG N-glycosylation trait groups

Multivariate analysis in the discovery cohort of 8090 samples revealed 32 loci associated with at least one group of IgG N-glycosylation traits at a genome-wide significant level of P < 5.6 × 10⁻⁹ (Table 1, quantile–quantile plots for all trait groups are given in Supplementary Material, Fig. S1 and the matrix of correlations between 23 IgG N-glycosylation traits in the discovery cohort is presented in Supplementary Material, Fig. S2a). These loci include the major histocompatibility complex (MHC) region and two loci (rs11847263 near FUT8 and rs10445337 near MAPT), which map in regions of extended linkage disequilibrium (LD) and therefore do not meet the locus definition ‘±500 kb around the lead SNP’ (regional association plots for these loci can be found in Supplementary Material, Fig. S3). Twenty-six loci were associated with ≥2 trait groups. The median number of trait groups associated with the identified loci was 6 (range 1–9). The most significantly associated glycan trait groups (‘top trait groups’) are presented in Table 1. All trait groups revealed in the multivariate analysis are listed in Supplementary Material, Table S2.

Table 1

Loci associated with IgG N-glycosylation trait groups in the multivariate analysis (discovery cohort)

Lead SNP	Chr: positiona	EA/RAb	Topc group of IgG N-glycosylation traits	β (SE)	P-value	EAFd	Ne	Number of traits groupsf	Nearest gene(s)g
Novel loci identified in a multivariate analysis
rs12049042	1:246288812	C/T	Galactosylation	0.02 (0.003)	1.20 × 10−09	0.01	4417	1	SMYD3
rs11895615	2:26113120	C/T	BisectingGlcNAc	0.01 (0.001)	5.69 × 10−10	0.70	8024	1	ASXL2
rs1372288	3:142901537	C/T	N-glycosylation	0.01 (0.002)	8.73 × 10−11	0.23	8090	2	CHST2, SLC9A9
rs12635457	3:196203979	A/G	N-glycosylation	0.02 (0.002)	1.61 × 10−13	0.64	8016	1	RNF168
rs479844	11:65551957	A/G	N-glycosylation	0.02 (0.002)	1.97 × 10−13	0.42	8089	3	OVOL1
rs4561508	17:16848750	C/T	N-glycosylation	0.01 (0.002)	1.38 × 10−10	0.90	8002	2	TNFRSF13B
Loci identified but not replicated in Klarić et al. (12)
rs6964421	7:6520676	C/T	Fucosylation	0.01 (0.002)	2.04 × 10−12	0.67	8090	5	KDELR2
rs10096810	8:103545436	A/G	BisectingGlcNAc	0.01 (0.001)	4.69 × 10−09	0.63	8090	1	ODF1
rs540719	16:23456578	C/T	Fucosylation	0.01 (0.002)	8.20 × 10−11	0.19	8087	2	COG7
rs7257072	19:19267990	C/T	Fucosylation	0.01 (0.002)	5.83 × 10−14	0.49	8090	7	MEF2B, BORCS8-MEF2B
rs2745851	20:17829280	A/G	Sialylation	0.02 (0.002)	4.19 × 10−20	0.38	8087	9	MGME1, BANF2, SNX5
Loci identified and replicated in at least one previous study (9–12)
rs311438	1:25320236	A/G	Digalactosylation	0.01 (0.001)	1.66 × 10−09	0.09	8070	1	MIR4425, RUNX3
rs6764279	3:186723975	A/C	Galactosylation	0.22 (0.006)	3.51 × 10−302	0.28	8086	9	ST6GAL1
rs3777184	5:95262044	C/G	N-glycosylation	0.03 (0.002)	1.22 × 10−27	0.26	8084	8	ELL2
rs4705950	5:131793286	C/T	N-glycosylation	0.01 (0.002)	1.99 × 10−09	0.58	8010	1	IRF1-AS1, IRF1, SLC22A5
rs3099844h	6:31448976	A/C	Fucosylation	0.03 (0.003)	2.81 × 10−20	0.14	8090	8	MICB-DT, MICB, HCG26
rs9385856	6:139625074	C/T	Fucosylation	0.02 (0.002)	1.16 × 10−17	0.43	8089	8	TXLNB
rs7758383	6:143169723	A/G	Fucosylation	0.02 (0.002)	4.86 × 10−28	0.52	8020	6	HIVEP2
rs7789913	7:50352695	C/T	Galactosylation	0.03 (0.002)	3.63 × 10−36	0.61	8089	8	IKZF1
rs7812088	7:150919829	A/G	N-glycosylation	0.03 (0.003)	6.32 × 10−35	0.12	8089	8	ABCF2
rs10124479	9:33136233	G/T	Fucosylation	0.06 (0.003)	7.52 × 10−76	0.32	8070	9	B4GALT1
rs481080	11:114442265	A/G	N-glycosylation	0.02 (0.002)	5.43 × 10−14	0.48	8090	5	NXPE4, NXPE2
rs11847263i	14:65775695	G/T	N-glycosylation	0.06 (0.003)	6.92 × 10−71	0.35	8044	9	MIR4708, FUT8
rs4074453	14:105998544	C/T	Fucosylation	0.04 (0.003)	1.55 × 10−51	0.24	8089	7	TMEM121
rs4795400	17:38067020	C/T	Monogalactosylation	0.01 (0.001)	3.41 × 10−13	0.53	7978	3	IKZF3, GSDMB, ORMDL3, LRRC3C, ZPBP2
rs10445337i	17:44067400	C/T	N-glycosylation	0.02 (0.002)	1.62 × 10−22	0.21	8076	6	MAPT
rs11650354	17:45822092	C/T	Galactosylation	0.01 (0.002)	3.34 × 10−10	0.83	8090	2	TBX21
rs9319617	17:79192446	C/T	Fucosylation	0.02 (0.002)	5.95 × 10−19	0.50	7750	6	CEP131
rs12986368	19:5830922	C/T	N-glycosylation	0.02 (0.002)	3.45 × 10−21	0.63	7717	4	FUT6
rs7281587	21:36565278	A/G	N-glycosylation	0.02 (0.002)	5.59 × 10−13	0.25	8090	5	RUNX1, RUNX1-IT1
rs17630758	22:24136542	A/G	N-glycosylation	0.03 (0.003)	3.31 × 10−41	0.15	8041	8	SMARCB1
rs5750830	22:39840828	A/C	N-glycosylation	0.08 (0.004)	8.11 × 10−101	0.74	7867	8	TAB1, MGAT3

β (SE), P-values and N are given for associations with top group of IgG N-glycosylation traits. GlcNAc, N-Acetylglucosamine. Associated loci were defined as regions within ±500 kb around the lead SNPs.

^aChromosome: position on chromosome according to GRCh37.p13 assembly.

^bEffect allele/reference allele.

^cTrait group associated with SNP at the highest level of statistical significance.

^dEffect allele frequency (discovery cohort).

^eNumber of subjects in a multivariate analysis.

^fTotal number of traits groups associated with the locus.

^gNearest gene(s) according to the dbSNP (https://www.ncbi.nlm.nih.gov/snp/).

^hMHC locus (chromosome 6, 25–35 Mb).

ⁱFor these loci, regional association plots were manually inspected, and neighboring loci were merged (three loci merged into one for the locus marked by rs10445337; two loci merged into one for the locus marked by rs11847263; Supplementary Material, Figure S3).

The six loci we discovered were not reported as genome-wide significant in any previous IgG N-glycome GWAS and were therefore considered novel; although the locus rs4561508 near TNFRSF13B reached a suggestive level of significance in the most recent GWAS by Klarić et al. (12). Their regional association plots are presented in Supplementary Material, Figure S4. Twenty six other loci overlap with those identified previously (which means that the lead SNPs corresponding to these 26 loci reside within the genome-wide significant loci reported in the studies (9–12)). Of those, 21 were replicated in at least one previous GWAS (9–12), and 5 were first discovered in Klarić et al.’s GWAS (12), but their associations were not replicated in the UPLC-measured replication cohort of that study.

Replication, which is a reproduction of statistically significant associations in an independent set of samples, is an important step in GWAS. Replication reduces the likelihood that the observed genotype–phenotype association is a chance finding or an analysis artifact (17) and is, therefore, essential before undertaking complex and expensive experimental follow-up studies. In our replication analysis, we included six novel loci discovered at the genome-wide significance level in this study and the five loci identified, but not replicated previously. Our replication analysis used 3147 samples from three cohorts. From 11 loci analyzed, 7 [5 novel and 2 discovered by Klarić et al. (12)] passed our criteria for replication: (1) were associated at P < 0.0045 with the discovery top trait group in the MANOVA replication test and (2) were associated at P < 0.0045 with the linear combinations of glycan traits (constructed in the replication cohort using the coefficients estimated in the discovery analysis), with the beta sign for the effect allele (EA) consistent with that observed in the discovery analysis (Table 2). Details are provided in Supplementary Material, Table S2, Supplementary Material, Figure S2b (correlations between the traits in the replication cohort) and Supplementary Material, Table S3 (partial regression coefficients for individual IgG N-glycosylation traits estimated for the top trait groups in the discovery and replication cohorts). Thus, our study identified and replicated seven additional loci associated with IgG N-glycome, which constitutes nearly one-third of the number of loci that are already known to affect IgG N-glycosylation (n = 22). For these loci, we performed an in silico functional analysis (see below).

Table 2

Replication of loci associated with IgG N-glycosylation trait groups

Lead SNP	Chr: positiona	EA/RAb	Topc group of IgG N-glycosylation traits	PMANOVAd	βdiscover (SE)e	βlin.com (SE)f	Plin.comf	Ng	Nearest gene(s)h	Genes prioritized by in silico methodsi or suggested based on literature dataj
Novel loci identified in a multivariate analysis
rs12049042	1:246288812	C/T	Galactosylation	0.45	0.02 (0.003)	−0.13 (0.08)	0.11	2572	SMYD3	Not tested
rs11895615	2:26113120	C/T	BisectingGlcNAc	3.01 × 10−08	0.01 (0.001)	0.16 (0.02)	1.24 × 10−11	3147	ASXL2	ASXL2 (D)
rs1372288	3:142901537	C/T	N-glycosylation	2.06 × 10−13	0.01 (0.002)	0.17 (0.03)	5.06 × 10−10	3147	CHST2, SLC9A9	SLC9A9 (L)
rs12635457	3:196203979	A/G	N-glycosylation	1.13 × 10−04	0.02 (0.002)	0.08 (0.02)	4.77 × 10−04	3147	RNF168	RNF168 (L)
rs479844	11:65551957	A/G	N-glycosylation	2.69 × 10−14	0.02 (0.002)	0.20 (0.02)	5.62 × 10−18	3147	OVOL1	AP5B1 (D), OVOL1 (S), EFEMP2 (S*)
rs4561508	17:16848750	C/T	N-glycosylation	3.27 × 10−04	0.01 (0.002)	0.13 (0.04)	2.91 × 10−04	3147	TNFRSF13B	TNFRSF13B (V, D)
Loci identified but not replicated in Klarić et al. (12)
rs6964421	7:6520676	C/T	Fucosylation	7.20 × 10−03	0.01 (0.002)	0.10 (0.02)	5.35 × 10−05	3147	KDELR2	Not tested
rs10096810	8:103545436	A/G	BisectingGlcNAc	0.09	0.01 (0.001)	0.07 (0.02)	2.75 × 10−03	3147	ODF1	Not tested
rs540719	16:23456578	C/T	Fucosylation	0.12	0.01 (0.002)	0.03 (0.03)	0.27	3147	COG7	Not tested
rs7257072	19:19267990	C/T	Fucosylation	5.61 × 10−05	0.01 (0.002)	0.13 (0.02)	8.19 × 10−08	3147	MEF2B, BORCS8-MEF2B
rs2745851	20:17829280	A/G	Sialylation	1.65 × 10−04	0.02 (0.002)	0.20 (0.04)	7.91 × 10−06	1127	MGME1, BANF2, SNX5

The loci which were not replicated are shown in italics. More details are provided in Supplementary Material, Table S2. GlcNAc, N-Acetylglucosamine.

^aChromosome: position on chromosome according to GRCh37.p13 assembly.

^bEffect allele/reference allele.

^cTrait group associated with SNP at the highest level of statistical significance in the discovery cohort.

^dP-value for the association with the top trait group in the replication cohort (MANOVA test).

^eBeta (SE) for the association with IgG N-glycosylation trait groups in the multivariate analysis in the discovery cohort.

^fBeta (SE) and P-value for the association with the linear combination of glycan traits in the replication cohort (constructed using the coefficients estimated in the discovery cohort; each coefficient corresponds to the trait in the discovery top trait group).

^gNumber of subjects in the replication analysis.

^hNearest gene(s) according to the dbSNP (https://www.ncbi.nlm.nih.gov/snp/).

ⁱGenes prioritized in the VEP analysis (V; Supplementary Material, Table S4), SMR/HEIDI analysis (S, likely shared, P_HEIDI ≥ 0.05; S^*, possibly shared, 0.05 > P_HEIDI ≥ 0.001; Supplementary Material, Table S5a) and DEPICT analysis (D; Supplementary Material, Table S6a).

^jGenes that can be considered potentially causal based on published data on their function (L).

Functional annotation of replicated loci

Variant Effect Predictor analysis

The results of the Ensembl Variant Effect Predictor (VEP) analysis (22) are shown in Supplementary Material, Table S4. The analysis identified two missense SNPs: rs34562254 in the TNFRSF13B (TNF receptor superfamily member 13B) gene and rs3796129 in the RNF168 (ring finger protein 168) gene. These polymorphisms are in high LD with the lead SNPs rs4561508 and rs12635457 associated with the trait group ‘N-glycosylation’ in our multivariate analysis. The missense variant rs34562254 A is correlated with the allele rs4561508 T (r² = 0.89 in European ancestry populations according to LDlink (23), https://analysistools.nci.nih.gov/LDlink/), and the missense allele rs3796129 T is correlated with the allele rs12635457 A (r² = 0.94 in European ancestry populations). SNP rs3796129 is triallelic (G >C, T; Pro401Arg or Pro401Gln) with the allele C being extremely rare, and LD between rs3796129 C and rs12635457 alleles is unknown. The rs34562254 (Pro205Leu) was classified as possibly damaging by the PolyPhen algorithm and as tolerated by the SIFT tool. The rs3796129 was classified as benign/tolerated by both algorithms. Based on these data, we only prioritized the gene TNFRSF13B bearing the polymorphism rs34562254.

Colocalization analysis

Colocalization analysis was performed for five out of seven loci identified and replicated in our study (rs11895615, rs479844, rs4561508, rs7257072 and rs2745851). The loci rs1372288 and rs12635457 were not included in the analysis since they did not reach a genome-wide significance in the meta-analysis of univariate GWAS results for the 23 IgG N-glycosylation traits obtained in the discovery and replication cohorts. The results of the SMR/HEIDI (Summary data-based Mendelian Randomization analysis/Heterogeneity in Dependent Instruments test) colocalization analysis (24) with expression quantitative trait loci (eQTLs) are presented in Supplementary Material, Table S5a. The analysis suggested two genes whose expression levels may be associated with the same functional polymorphisms that affect IgG N-glycosylation: OVOL1 (ovo-like transcriptional repressor 1) and EFEMP2 (EGF containing fibulin extracellular matrix protein 2). The strongest result (nominal P_SMR = 6 × 10⁻⁸, Bonferroni-corrected P_SMR = 5 × 10⁻⁵, P_HEIDI = 0.4) was obtained for the locus marked by rs479844 where increased expression of OVOL1 in Epstein–Barr virus (EBV)-transformed lymphocytes was found to be associated with the same allele as associated with the increased percentage of FA2FG2S1 glycan in total IgG glycans (IGP19 UPLC peak). For the same locus, we also observed colocalization with the expression of the EFEMP2 gene in peripheral blood; however, the SMR test significance was weaker (Bonferroni-corrected P_SMR = 0.003), and the sharing was characterized as ‘possible’ (P_HEIDI = 0.03). Based on these observations, we consider OVOL1 and EFEMP2 as possible candidate genes, which may influence IgG N-glycosylation or its regulation.

We also used colocalization analysis to check for overlap between the IgG N-glycosylation loci and the loci affecting the risk of autoimmune and inflammatory diseases (asthma, rheumatoid arthritis, primary biliary cirrhosis, ankylosing spondylitis and inflammatory bowel disease, including Crohn’s disease and ulcerative colitis). For the locus marked by rs479844 (near OVOL1), we detected colocalization with self-reported asthma (the UK Biobank phenotype). The same allele was found to be associated with the increased percentage of FA2FG2S1 glycan in total IgG glycans (IGP19 UPLC peak) and decreased asthma risk (Supplementary Material, Table S5b). For other diseases and loci analyzed, colocalization was not established.

Data-driven Expression Prioritized Integration for Complex Traits analysis

The genes prioritized by the Data-driven Expression Prioritized Integration for Complex Traits (DEPICT) (25) analysis are listed in Supplementary Material, Table S6a. Statistically significant results (FDR < 0.05) were only observed when the P-value threshold for input SNPs was set at 5.0 × 10⁻⁸. DEPICT prioritized three genes in three loci replicated in our study: ASXL2 (ASXL transcriptional regulator 2; locus marked by rs11895615), AP5B1 (adaptor related protein complex 5 subunit beta 1; locus marked by rs479844) and TNFRSF13B (locus marked by rs4561508).

DEPICT gene set enrichment analysis provided results with FDR < 0.05 only when we included all SNPs associated with IgG N-glycosylation at P < 1.0 × 10⁻⁵. The only enriched category identified was ‘decreased mature B cell number’ (Supplementary Material, Table S6b). DEPICT tissue/cell type set enrichment analysis yielded results with FDR < 0.05 only when the P-value threshold for input SNPs was set at 5.0 × 10⁻⁸. Nineteen tissues/cell types were revealed, and almost all of them were related to the hemic and immune systems: blood, palatine tonsil, lymphocytes, leukocytes, antibody producing cells, plasma cells, natural killer cells, etc. (Supplementary Material, Table S6c).

Summary of in silico gene prioritization

In silico functional analysis allowed us to suggest a total of five potentially causal genes for the three loci marked by rs11895615, rs479844 and rs4561508 (Table 2). For the locus rs11895615, a single ASXL2 gene (nearest to the lead SNP) was proposed by the DEPICT tool. For the loci rs479844 and rs4561508, three genes and one gene were proposed, respectively, of which the gene TNFRSF13B (for rs4561508) was prioritized using two approaches (eQTL colocalization analysis and VEP analysis), AP5B1 (for rs479844) was proposed only by the DEPICT tool and OVOL1 and EFEMP2 (for rs479844) were suggested only by the eQTL colocalization analysis. Colocalization with eQTL, on the one hand, can provide evidence for expression-mediated effects on the analyzed trait and, on the other hand, is not able to distinguish causation from pleiotropy and in some cases indicates a group of genes whose expression may be affected by the same regulatory locus, or even a group of functional SNPs in very tight LD. We therefore avoid making definite conclusions on the best candidate genes for the locus rs479844.

Discussion

Identification of new loci associated with IgG N-glycosylation can provide clues towards understanding the mechanisms of in vivo regulation of this important biological process intimately connected to immune system function. These leads and hypotheses should be followed up in subsequent targeted functional experiments. However, before undertaking lengthy and expensive experimental work, much effort should be invested in excluding the possibility of false-positive associations by performing replication in independent samples as well as conducting in silico functional annotation.

In this study, we performed a multivariate genome-wide association analysis of 23 plasma IgG N-glycosylation traits having combined them into nine groups based on general structural characteristics of IgG-attached glycans (see Fig. 1; one group included all 23 traits). A multivariate approach can increase the statistical power compared to standard univariate GWAS (10,14,15), and here we demonstrate this by revealing six novel genome-wide significant loci, which have not been identified in the study of Klarić et al. (12), providing us with the GWAS summary statistics for individual glycan traits. Of note, 26 out of 27 loci revealed in that study (12) (all except rs12341905 near SPINK4) also reached genome-wide significance in our work (Table 1, Supplementary Material, Table S2). We replicated five out of six of the novel loci as well as an additional two loci that have not been replicated previously. Thus, our study increases the number of established IgG N-glycosylation loci from 22 to 29.

It may be interesting to discuss why in our analysis the association with rs12341905 near SPINK4—a locus not only detected but also replicated by Klarić et al. (12)—was not genome-wide significant. First, we should mention that while a multivariate analysis increases power for identification of genetic associations with glycan traits, this is a statistical relationship that holds asymptotically but is not guaranteed to hold for every specific test. This is especially true for scenarios under which the gain of power by the multivariate analysis is minimal (10,16). Note also that in the study of Klarić et al. (12) the association with rs12341905 was the weakest among 27 others and had barely passed the genome-wide significance threshold. Therefore, it may be not too surprising that in our study that used another analysis method this SNP has just missed the study-specific genome-wide significance threshold: the P-value for the trait ‘monosialylation’ was 1.46 × 10⁻⁸, which is greater than 5.6 × 10⁻⁹. Still, we should emphasize that, given replication of this locus by Klarić et al. (12), we consider this locus as established, and we attribute the fact that we could not detect this locus at genome-wide significance using multivariate analysis to the stochastic variation in power of the methods.

For the seven newly established loci, we have performed an extensive in silico follow-up analysis. None of these loci contain genes encoding glycosyltransferases. For five of the seven replicated loci, we suggested seven candidate genes based on the analysis of literature (two genes) and/or through use of several bioinformatics tools and methods (five genes) (Table 2), with the TNFRSF13B gene (locus rs4561508) being prioritized using two independent approaches.

Among the genes prioritized as possible candidates, several could be linked to N-glycosylation and/or immune function based on current knowledge.

In the locus tagged by rs1372288, the top associated SNP is in LD with rs4839604 (r² = 0.79 in European ancestry populations according to LDlink) associated with plasma protein N-glycosylation (26,27). IgG is among the most abundant plasma glycoprotein (1,2), and we assume that changes in its glycosylation patterns may be responsible for the effects observed for the total plasma protein N-glycome. These previous studies (26,27) did not apply in silico methods to find potentially causal genes but prioritized the gene SLC9A9 (carrier family 9 member A9) based on its possible involvement in the regulation of endosomal pH (28).

The TNFRSF13B gene, which was suggested by both the VEP and DEPICT tools, encodes a lymphocyte-specific member of the tumor necrosis factor receptor superfamily known as the TACI protein (the transmembrane activator and calcium modulator and cyclophilin ligand interactor). Signaling through this protein activates several transcription factors, including NFAT, AP-1 and NF-kappa-B (29). TACI plays a crucial role in humoral immunity acting in synergy with other signaling networks to promote B cell differentiation into plasma cells and immunoglobulin production (30,31). Coding variants in TNFRSF13B are associated with common variable immunodeficiency and selective IgA deficiency (32), rare TNFRSF13B mutations can increase the risk of asthma symptoms (33) and several common variants in this gene were reported to be associated with circulating IgG levels in GWASes conducted by Liao et al. (34) and Jonsson et al. (35). One of the lead SNPs in the study of Liao et al., rs4792800, is in LD with our lead N-glycosylation-associated SNP rs4561508 (r² = 0.95 in European ancestry populations according to LDlink). These findings indicate a possible pleiotropic action of TNFRSF13B with respect to N-glycosylation of IgG and immune function.

The gene RNF168 (ring finger protein 168) bearing rs12635457 in its intron encodes an E3 ubiquitin ligase protein involved in DNA double-strand break repair and immunoglobulin class switch recombination (36). Mutations in this gene lead to a RIDDLE syndrome (an immunodeficiency and radiosensitivity disorder), the features of which include decreased immunoglobulin production (37). An association between SNPs at/near RNF168 and IgM levels was previously reported (35). As far as we are aware, there is no reported evidence for RNF168 involvement in glycosylation processes, but if so, its pleiotropic action on immune-related and glycosylation traits may be assumed.

Among IgG N-glycosylation loci established in this work, the locus tagged by rs479844 near OVOL1 is the most interesting in our opinion. We demonstrated the potential pleiotropic effect of this locus on asthma (Supplementary Material, Table S5b); the same locus was identified in GWAS for the combined eczema plus asthma phenotype (the atopic march) (38), combined phenotypes of asthma/hay fever/eczema (39) and atopic dermatitis (40,41). This locus was also found to have pleiotropic effects on dermatological and allergic diseases (42) and to be differentially methylated in childhood asthma (43). Three genes, AP5B1, OVOL1 and EFEMP2, were proposed for the rs479844 locus by our bioinformatic analysis, with colocalization results for OVOL1 being statistically strong and AP5B1 being prioritized by the DEPICT gene prioritization tool. AP5B1 (also known as DKFZp761E198) is a subunit of an adaptor protein complex AP-5 (44) likely involved in facilitating late endosome-to-Golgi retrieval—a pathway contributing to a normal lysosomal function (45). OVOL1 genes encode a putative zinc finger containing transcription factor, which is involved in the development and differentiation of epithelial tissues and germ cells (46–48) and plays a role in a skin barrier function (49). Finally, EFEMP2 encodes the EGF-containing fibulin-like extracellular matrix protein 2 (also known as fibulin-4). Rare mutations in this gene lead to cutis laxa type IB, an autosomal recessive disorder characterized by the presence of severe systemic connective tissue abnormalities [phenotype number in the Mendelian Inheritance in Man (MIM, https://www.omim.org/) database: 614437] (50,51). A potential link between the genes AP5B1, OVOL1, EFEMP2 and IgG N-glycosylation is not apparent. Therefore, we conclude that further research is needed to define the causal gene in the rs479844 locus.

Besides searching for causal genes and colocalization analysis, the functional annotation in our study involved enrichment analysis and assessment of colocalization with a set of inflammatory and autoimmune diseases. We observed an enrichment of tissue/cell/gene set categories related to the hemic and immune systems (Supplementary Material, Table S6b, S6c).

An intriguing pattern that emerges from this and previous studies is the overlap between loci that regulate IgG N-glycosylation and immune function. In previous works (9–12), a total of 22 IgG N-glycome-associated loci have been identified and replicated. Of these, five contained genes encoding N-glycosyltransferases. Of the other 17, 11 loci contained candidate genes implicated in the immune response. These loci, genes and proteins are as follows: the human leukocyte antigen (HLA) region, which contains multiple immune-related genes; HIVEP2 involved in cellular immunity (52,53); transcription factors RUNX1, RUNX3, IKZF1 and IKZF3 that play a central role in B-cell maturation and differentiation; VPREB3, which likely participates in the transport and assembly of the pre-B cell receptor (54); a transcription factor TBX21, mouse homolog of which regulates the expression of IFNγ and controls Th1 lineage commitment (55); a transcription elongation factor ELL2, which directs immunoglobulin secretion in plasma cells (56); a locus near the IRF1 (interferon regulatory factor 1) gene and a locus near IGH that codes for immunoglobulin heavy chains. Not surprisingly, the enrichment analysis highlighted the gene sets that are known to be connected to various processes involved in immunity, such as B-cell life cycle and antibody production (12). Moreover, for two loci, pleiotropic effects were demonstrated on the pro-inflammatory pattern of IgG N-glycosylation and the risk of inflammatory and autoimmune diseases such as Crohn’s (locus near IRF1-SLC22A4), ulcerative colitis, rheumatoid arthritis and primary biliary cirrhosis (locus near ORMDL3-GSDMB-IKZF3). Interestingly, the latter locus also exhibited pleiotropy with the risk of asthma, although the asthma risk allele was ‘protective’ for autoimmune disease. This is a well-known phenomenon for the 17q12-21 locus (13).

Among five novel loci detected and replicated in the present study, three can be connected to the immune system’s function. Given that this work has detected a second locus that is likely to be shared between IgG N-glycosylation and asthma, it may be tempting to speculate about biological mechanisms linking these traits. However, we believe that it is too early to draw strong conclusions, and more data should be accumulated by more powerful GWAS. Several observations suggest that the mechanism linking asthma and IgG N-glycome via the loci near ORMDL3-GSDMB-IKZF3 and OVOL1 are different and support a cautious approach to interpretation of the results. First, different IgG N-glycan traits are associated with these loci. Second, the mechanisms by which these loci alter asthma risk are not likely to be the same. For example, the ORMDL3-GSDMB-IKZF3 locus was found to be more strongly associated with asthma without hay fever/eczema than to hay fever/eczema without asthma; the reverse is true for the locus near OVOL1 (39), which is also strongly implicated in atopic dermatitis (40,41).

Our findings are consistent with the conclusion drawn in the previous studies (9,10,12) that the loci and genes associated with IgG N-glycosylation are also related to immune function and immune-related conditions. In general, our results provide additional information characterizing the interplay between glycosylation and immune function.

In conclusion, our study establishes seven new loci associated with plasma IgG N-glycosylation and suggests a number of genes potentially involved in this process. Some loci and genes (TNFRSF13B, RNF168 and OVOL1/AP5B1) are likely to have pleiotropic effects on the traits related to the immune system and inflammatory diseases. For the locus near OVOL1, we demonstrate colocalization with variants affecting the risk of asthma, and the genes TNFRSF13B and RNF168 have previously been implicated in immunity disorders and immunoglobulin production. Future studies should focus on increasing statistical power to provide better resolution for the in silico studies and on providing mechanistic insights into the observed associations by targeted functional follow-up. Besides this, it will be fruitful to investigate causative relationships between alteration of IgG N-glycosylation and immune-related conditions.

Materials and Methods

Study cohorts

Ethics statement

All studies providing genotype and phenotype data were approved by local Ethics Committees/Institutional Review Boards. Ethics approvals were given in compliance with the principles expressed in the Declaration of Helsinki. All study participants have signed written informed consent.

Discovery cohort

In the discovery stage of our analysis, we used GWAS summary statistics for plasma IgG N-glycosylation traits obtained in the study by Klarić et al. (12) (discovery dataset, European ancestry individuals, N = 8090; GWAS summary statistics can be downloaded from (57)). In that study, GWAS was conducted for a total of 77 quantitative IgG N-glycosylation traits—23 traits directly measured by ultra-performance liquid chromatography and 54 derived traits combining directly measured glycans in a biologically meaningful way. The full list of traits is given in Supplementary Material, Table S1. Description of the analyzed cohorts and details on IgG isolation from plasma samples, IgG N-glycome quantification, genotyping, glycan and genetic data pre-processing, quality control, imputation, genome-wide association analyses and meta-analysis are provided in Klarić et al. (12).

Replication cohort

Replication was performed using GWAS data obtained for a total of 3147 samples from four cohorts of European descent: EGCUT (N = 575), FINRISK (N = 552), CROATIA-Korcula2 (N = 941) and VIKING (N = 1079). Description of the EGCUT and FINRISK cohorts as well as details on phenotyping- and genotyping-related procedures can be found in the study of Klarić et al. (12).

The Viking Health Study—Shetland (VIKING) is a cross-sectional volunteer study from the isolated archipelago of Shetland in Northern Scotland, with high kinship and considerable genetic drift (58).

The CROATIA-Korcula2 study is a population-based study that sampled participants, aged 18–92 years, from the Adriatic island of Korčula in 2012 as an extension of the 10 001 Dalmatians Biobank (59). All subjects visited the clinical research center in the region, where they were examined in-person and where fasting blood was drawn and stored for future analyses. Many biochemical and physiological measurements were performed, and questionnaires of medical history as well as lifestyle and environmental exposures were collected.

GWAS results for replication cohorts were combined by the inverse-variance-weighted meta-analysis (fixed-effects model). Analysis was performed using the GWAS-MAP platform (60). SNPs with MAF < 1% were excluded. Before meta-analysis, a genomic control (61) was applied to the summary statistics of each study.

Phenotypic correlations

Phenotypic correlations between 23 IgG N-glycosylation traits were estimated as correlations between the vectors of Z-statistics obtained in GWAS for each trait (method described by Stephens (15)). Pearson correlations were computed using the R statistical software (62) version 3.4.4.

Multivariate analysis

Multivariate analysis was performed using the summary statistics-based method described by Ning et al. (63). This method enables multivariate analysis of variance (MANOVA) using the results of multiple GWAS for individual traits. Calculations were performed using the MultiABEL package version 1.1–9 (https://cran.r-project.org/src/contrib/Archive/MultiABEL/, tutorial: https://github.com/xiashen/MultiABEL/; ‘MultiSummary’ function) for R software environment (64).

Following Shen et al. (10), the statistical significance threshold for multivariate analysis was set at P-value < 5.6 × 10⁻⁹ (5.0 × 10⁻⁸/9, where 9 is the number of glycan groups).

Associated loci were defined as regions within ±500 kb around the SNPs most significantly associated with IgG N-glycosylation trait groups (lead SNPs). SNPs located between 25 and 35 Mb on the chromosome 6 were merged into a single MHC locus. Additionally, for neighboring loci, regional association plots were manually inspected, and the loci were merged if they fell into one LD block.

Replication

Replication was performed according to a protocol similar to that described in Shen et al. (10). First, summary statistics obtained in the meta-analysis of replication cohorts was used to get MANOVA results for each lead SNP and each ‘top’ trait group (the group of IgG N-glycosylation traits most significantly associated with the lead SNP in the discovery cohort). Analysis was performed using the MultiABEL R package version 1.1-9 (‘MultiSummary’ function). Association was considered replicated if the association between the lead SNP and the corresponding trait group in the replication MANOVA test passed the Bonferroni-corrected threshold of P < 0.0045 (0.05/11, where 11 is the number of loci in the replication analysis; the first criterion of replication). Next, for each lead SNP-top trait group pair, we extracted the coefficients of linear combinations of individual traits in the discovery cohort and used them to construct the corresponding trait group in the replication cohort. Associations between the lead SNPs and the constructed linear combinations were tested as described in Tsepilov et al. (64) and McDaid et al. (65). In brief, the effects and corresponding standard errors (SEs) were calculated as follows:

$$\begin{array}{rcl} {\beta}_{\mathrm{new}}&=&\frac{B\ {}^{\ast} a}{\sqrt{\mathrm{sum}\left({\Sigma}_{\mathrm{ph}}\ {}^{\ast} \left(a\otimes a\right)\right)}},\\ \\{\mathrm{SE}}_{\mathrm{new}}&=&\frac{\sqrt{\left({\mathrm{SE}}_1^2+\frac{\beta_1^2}{N}\right)\ {}^{\ast} \mathrm{sum}\left({\Sigma}_{\mathrm{ph}}\ {}^{\ast} \left(a\bigotimes a\right)\right)-\frac{{\left(B\ {}^{\ast} a\right)}^2}{N}}}{\sqrt{\mathrm{sum}\left({\Sigma}_{\mathrm{ph}}\ {}^{\ast} \left(a\bigotimes a\right)\right)}},\mathrm{respectively}; \end{array}$$

where B—matrix, where each column β_i—effect of each SNP in GWAS_i; a—vector of linear combination coefficients estimated in the discovery cohort; Σ_ph—matrix of correlations between 23 IgG N-glycosylation traits in the replication cohort; |$\bigotimes$|—outer product; SE₁—standard error of each SNP’s effect in the first GWAS₁; N—minimal N out of all GWASes for each SNP. If the association of the EA with the constructed linear combinations had the same direction of effect as in the discovery cohort and passed the Bonferroni-corrected threshold of P < 0.0045, the locus was considered replicated (the second criterion of replication).

In silico functional annotation

Prediction of SNP effects

We used the Ensembl VEP (22) to analyze the potential effects of SNPs and indels in high LD (r² > 0.8) with the replicated polymorphisms. LD was calculated using PLINK 1.9 (66) (--show-tags option) and genotype data for 503 European ancestry individuals (1000 Genomes phase 3 version 5 data). The analysis was performed using software available online (https://www.ensembl.org/info/docs/tools/vep/index.html).

Colocalization analysis

In order to find potentially causal genes involved in IgG N-glycosylation and to investigate pleiotropic effects on autoimmune and inflammatory conditions, we analyzed the colocalization between the loci identified and replicated in our multivariate analysis, eQTLs from selected tissues and loci associated with asthma, rheumatoid arthritis, primary biliary cirrhosis, ankylosing spondylitis, inflammatory bowel disease, Crohn’s disease and ulcerative colitis [these traits were selected in accordance with the similar analysis performed in the study of Klarić et al. (12)]. Colocalization was examined using SMR analysis followed by the HEIDI test (24). This approach allows to identify the genes whose expression level may be affected by the same causal SNP that is associated with the analyzed trait as well as to reveal pleiotropic loci influencing multiple traits or diseases. The SMR test indicates whether two traits (e.g. IgG N-glycosylation and gene expression level in a certain tissue or the disease) are associated with the same locus, and HEIDI test specifies whether both traits are affected by the same underlying functional SNP.

Results of multivariate GWAS cannot be directly used for a colocalization analysis. Therefore, we have meta-analyzed univariate GWAS results for the 23 IgG N-glycosylation traits obtained in discovery and replication cohorts using the inverse-variance-weighted meta-analysis method (fixed-effects model). For each locus included in the colocalization analysis, we selected the ‘primary’ IgG N-glycosylation GWAS from the resulting meta-analyses. The selection criterion was the lowest P-value for the association between the lead SNP marking the locus and IgG N-glycosylation trait for which the meta-analysis was conducted. Loci where associations did not reach P-value <5.0 × 10⁻⁸ (rs1372288 and rs12635457) were excluded from the analysis.

Summary statistics for gene expression levels in tissues/cell types were obtained from the Blood eQTL study (67) (http://cnsgenomics.com/software/smr/#eQTLsummarydata), the CEDAR project (68) (http://cedar-web.giga.ulg.ac.be/) and the GTEx project version 7 (69) (https://gtexportal.org). In total, we used data for eight tissues/cell types: peripheral blood (from Blood eQTL study); six circulating immune cell types—CD4+ T lymphocytes, CD8+ T lymphocytes, CD19+ B lymphocytes, CD14+ monocytes, CD15+ granulocytes and platelets (from CEDAR); and the EBV-transformed lymphocytes (from GTEx).

Summary statistics for autoimmune/inflammatory diseases were obtained from the following resources: for asthma (self-reported), we used the UK Biobank (70) GWAS summary statistics provided by the Neale Lab (http://www.nealelab.is/; 39 049 case and 298 110 control individuals; data were downloaded in 2017); for rheumatoid arthritis, GWAS summary statistics from Okada et al.’s study (71); for primary biliary cirrhosis, from Cordell et al.’s study (72); for ankylosing spondylitis, from Cortes et al.’s study (73) and for inflammatory bowel disease and its subtypes, from Liu et al.’s study (74).

SMR/HEIDI analysis was performed according to the protocol described by Zhu et al. (24). We used sets of SNPs having the following properties: (1) being located within ±250 kb from the lead SNPs replicated in the present study; (2) being present in both the primary GWAS and eQTL data/GWAS for the disease; (3) having MAF ≥ 0.03 in both datasets; (4) having squared Z-test value ≥10 in the primary GWAS. Those SNPs that met criteria (1), (2), (3) and (4) had the lowest P-value in the primary GWAS and were in high LD (r² > 0.8) with the lead SNPs were used as instrumental variables to elucidate the relationship between gene expression/disease and IgG N-glycosylation (we define them as ‘top SNPs’). It should be noted that SMR/HEIDI analysis does not identify a causative SNP affecting both traits. It can be either the top SNP or any other SNP in strong LD.

After defining the set of eligible SNPs for each locus, we made the ‘target’ and ‘rejected’ SNP sets and added the top SNP to the ‘target’ set. Then, we performed the following iterative procedure of SNP filtration: if the SNP from the eligible SNP set with the lowest P_SMR had r² > 0.9 with any SNP in the ‘target’ SNP set, it was added to the ‘rejected’ set; otherwise, it was added to the ‘target’ set. The procedure was repeated until eligible SNP set was exhausted, or the ‘target’ set had 20 SNPs. If we were unable to select three or more SNPs, the HEIDI test was not conducted. HEIDI statistics was calculated as |${T}_{\mathrm{HEIDI}}={\sum}_i^m{z}_{d(i)}^2,$| where m is the number of SNPs selected for analysis,

$$ {z}_{d(i)}=d_i\!\Big/ \!{\mathrm{SE}}_{({d}_{i})}$$

 

$$ {d}_i={\beta}_{{\mathrm{SMR}}_i}-{\beta}_{\mathrm{SMR}\ \left(\mathrm{top}\ \mathrm{SNP}\right)} $$

All SNP effects and SEs were standardized. Analysis was performed using the GWAS-MAP platform (60).

In the analysis of pleiotropy, we explored colocalization of regional associations of IgG N-glycosylation and selected diseases and gene expression traits in selected tissues/cell types. Similar to the previous works (12,75), the results of the SMR test were considered statistically significant if P_SMR < 6.2 × 10⁻⁵ (0.05/801, where 801 is a total number of tests corresponding to analyzed loci and gene expression/disease traits). Inference of whether a functional variant may be shared between the IgG N-glycan trait and gene expression/disease was made based on the HEIDI test: P_HEIDI ≥ 0.05 (likely shared), 0.05 > P_HEIDI ≥ 0.001 (possibly shared) and P_HEIDI < 0.001 (sharing is unlikely).

DEPICT analysis

Gene prioritization and gene set and tissue/cell type enrichment analyses were performed using the DEPICT (25) software version 1.1, release 194 with default parameters (https://data.broadinstitute.org/mpg/depict/). Analysis was conducted for SNPs associated with any IgG N-glycosylation trait group at P < 5.0 × 10⁻⁸ and P < 1.0 × 10⁻⁵. The MHC region was excluded. The significance threshold for DEPICT analysis was set at FDR < 0.05.

Data and code availability

All relevant data are within the manuscript and its Supporting Information files. All computer code used in this research is available at https://github.com/Defrag1236/MV_IGG.

Acknowledgements

The Viking Health Study. We would like to acknowledge the invaluable contributions of the research nurses in Shetland, the administrative team in Edinburgh and the people of Shetland.

The CROATIA-Korcula2 Study. We would like to acknowledge the staff of several institutions in Croatia that supported the fieldwork, including but not limited to The University of Split and Zagreb Medical Schools and Croatian Institute for Public Health.

Conflict of Interest statement. G.L. is a founder and owner of the Genos Glycoscience Research Laboratory, and O.O.Z. is an employee of the said company. Y.S.A. is a cofounder and a co-owner of PolyOmica and PolyKnomics, private organizations providing services, research and development in the field of computational and statistical genomics. The other authors declare that they have no competing interests.

Funding

Russian Science Foundation (grant number 19-15-00115 to A.S.Sh., S.Z.S., E.D.P and Y.S.A.); Ministry of Education and Science of the RF via the Institute of Cytology and Genetics (project 0259-2021-0009/AAAA-A17-117092070032-4 to A.S.Z.); Russian Ministry of Science and Education under the 5-100 Excellence Programme (to Y.A.T.); RCUK Innovation Fellowship from the National Productivity Investment Fund (MR/R026408/1 to L.K.); Croatian National Centre of Research Excellence in Personalized Healthcare grant (#KK.01.1.1.01.0010 to O.O.Z. and G.L.), Medical Research Council (UK) Quinquennial programme (to J.F.W.). The Viking Health Study—Shetland (VIKING) was supported by the MRC Human Genetics Unit quinquennial programme grant ‘QTL in Health and Disease’. DNA extractions and genotyping were performed at the Edinburgh Clinical Research Facility, University of Edinburgh. The CROATIA-Korcula study was funded by grants from the Medical Research Council (UK), and Republic of Croatia Ministry of Science, Education and Sports research grants (108-1080315-0302). The SNP genotyping for the CROATIA-Korcula2 cohort was performed in the core genotyping laboratory of the Clinical Research Facility at the Western General Hospital, Edinburgh, Scotland. The work of C.H. was supported by an MRC University Unit Programme Grant MC_UU_00007/10 (QTL in Health and Disease). The color figure charges were funded by the Russian Science Foundation grant number 19-15-00115.

Author notes

Gordan Lauc, Yurii S. Aulchenko, and Yakov A. Tsepilov contributed equally to this work