Abstract

Introduction

A large number of genetic loci have been mapped for complex traits of clinical relevance via genome-wide association studies (GWAS). However, our understanding of the biology behind the association is rather limited for most of these loci. In order to move further towards the goal of unraveling the biological mechanisms underpinning complex disease, large-scale omics data have been generated as intermediate phenotypes to help fill in the gap between the genome and organismal level traits. Nonetheless, in contrast to the complex traits that are commonly analyzed, the large number of omics variables often have a strong correlation structure that cannot be neglected¹. Therefore, instead of conventional univariate GWAS, multivariate analysis approaches are needed in GWA analyses to incorporate the ubiquitous, partly genetically regulated, correlation among omics variables.

Many multivariate methods have been developed for genetic analyses^{2, 3} (see Supplementary Note ¹ for details, also review by ref. ⁴). In GWAS context, multivariate omics analysis with explicit modelling of phenotypic covariance faces several challenges. In genetically structured samples, this structure needs to be taken into account. The multivariate methods allowing for covariates can handle population stratification by including principle components of the genomic kinship matrix in the model. However, the effects of small groups of close relatives—if present—are usually not reflected in the leading principal components because such kinship generates weaker LD than large-scale population structure does⁵. In such situations, the method of choice is mixed effect models; however, a fast mixed model analysis for multiple (more than tens of) phenotypes is still difficult. GEMMA⁶ and Limix⁷ were designed to overcome such difficulty, although the computational cost can be high for a large number of phenotypes such as that encountered in omics data. Another difficulty faced in multivariate analysis of omics phenotypes is finding a method for replication, which would allow for both tests of significance and consistency of the model being replicated.

Despite the number of methodological studies, few empirical multivariate GWAS have been published for humans. This may be largely attributed to the statistical complexity and computational difficulty of such analyses as well as outlined above issues addressing population stratification, replication, and interpretation. Among multivariate GWAS conducted to date, most are dealing with very complex phenotypes (e.g., refs ^{8, 9}). These studies usually propose a new method, and investigate real data to demonstrate that one may expect increased mapping power when using these. However, few new loci have been convincingly identified using such methods for complex phenotypes, and no replication was attempted. Therefore—at least for the moment—the practical gain from multivariate GWAS of complex traits seems to be rather limited. The situation appears to be somewhat different for traits more proximal to the genotype. Stephens¹⁰ applied multivariate GWAS to summary level data from the Global Lipids Genetics Consortium (GLGC-2010¹¹, N up to 100,184) and identified 18 new loci on the top of 95 loci identified in single trait analyses. While no replication was provided in this work, importantly, we observe that 11 of the loci identified by multivariate analyses were later confirmed with genome-wide significance in the next, bigger, univariate meta-analysis by GLGC in 2013¹² (N up to 188,577). It appears the yield of new loci is even bigger for metabolomics traits. Inouye et al.¹³ performed multivariate GWAS of 130 NMR metabolites (grouped in 11 sets) in 6600 individuals. The study demonstrated that multivariate analysis doubles the number of loci detected in this sample. Again, no replication was performed, however, among seven novel loci discovered by Inouye et al. via multivariate analysis, three demonstrate genome-wide significant (P<5×10⁻¹¹) association with at least one NMR metabolite investigated recently in a sample of up to 24,925 individuals¹⁴.

These observations strongly suggest that at least some of the novel loci discovered via multivariate GWAS are true positives. The fact that “univariate” replication of novel loci from the lipid and metabolomic studies is not perfect (around 50%), even after substantially increased sample size, suggests that either multivariate analyses are more prone to the appearance of statistical artifacts and/or are subject to different genome-wide significance thresholds, or that the underlying pleiotropic effects captured in multivariate models are hard to capture using standard univariate analyses (see ref. ¹⁰ for examples of such scenarios), and require different approach to replication.

Here, we explore the potential of multivariate methods for studying genetic regulation of glycan variation. Glycans are complex carbohydrates bound to the surface of proteins, whose structure and function they consequently substantially influence. Glycosylation is one of the most abundant post-translational protein modifications^{15, 16}, but knowledge about its biological function was long hindered by glycans’ structural complexity. Unraveling the complex network of genes involved in protein glycosylation can provide not only a better understanding of this fundamental biological process, but might also provide insights into how these molecules could be involved in complex human diseases, and potentially used as biomarkers in prediction of disease susceptibility^17–20.

Previous GWAS have analyzed either the N-glycans released from all plasma proteins^{21, 22} or focused on N-glycosylation of a single protein—Immunoglobulin G (IgG), which is first isolated from other plasma proteins, followed by quantification of enzymatically-released glycans²³. These studies uncovered six loci for total plasma and nine loci for IgG glycosylation^21–23; with only one locus (FUT8) overlapping between the two. The majority of plasma proteins are synthesized in the liver and pancreas²⁴ while immunoglobulins are synthesized specifically in cells of the immune system²⁵. The lack of overlap between results of these genetic studies suggests different mechanisms of biological control of glycosylation in these two tissues.

Here, we apply multivariate methods to the IgG N-glycosylation traits and empirically test whether novel loci can be convincingly replicated. To utilize data from populations with high kinship, we formulate and implement a multivariate GWAS workflow based on combination of a linear-mixed-model-based phenotypic transformation, MANOVA, and multiple regression.

Results

Joint analysis of IgG glycans identifies five novel loci

A single protein (IgG) was isolated from the plasma of 1960 individuals from the population of the Orkney Islands in Northern Scotland (the ORCADES cohort²⁶). The N-linked glycans were assayed using ultra performance liquid chromatography (UPLC), resulting in 23 quantitative measurements.

The heart of our analysis procedure is, essentially, MANOVA statistics (see Methods for details). However, in order to avoid the effect of confounding by population genetic structure and kinship, prior to MANOVA analysis each phenotype undergoes linear-mixed-model-based GRAMMAR+ transformation we have previously developed²⁷ as an improvement of GRAMMAR procedure²⁸ for rapid association analysis in pedigrees and samples form genetically isolated populations.

We performed multivariate GWAS combining all the 23 IgG N-glycosylation phenotypes using MANOVA. The same procedure was also applied to eight subsets of the traits based on different chemical and structural properties of glycans, namely galactosylation, monogalactosylation, digalactosylation, sialylation, monosialylation, disialylation, fucosylation and bisecting GlcNAc (N-acetyl glucosamine). The definition of subgroups is given in Supplementary Table ². The distribution of the observed test statistic was contrasted to that expected under the null hypothesis, similar to the genomic control²⁹ (Supplementary Fig. ¹). Given that we have analyzed nine groups of traits, although they are strongly overlapping, we decided to take a conservative approach and have considered nominal MANOVA P-values<(5×10⁻⁸)/9=5.6×10⁻⁹, as genome-wide significant.

In total, nine GWA scans were performed. Our multivariate analysis replicated five out of nine previously established²³ IgG N-glycosylation loci (ST6GAL1, B4GALT1, FUT8, SMARCB1-DERL3 and SYNGR1-TAB1-MGAT3; Fig. 1, Supplementary Table ³) and detected five new loci (IGH, ELL2, HLA-B-C, AZI1, FUT6-FUT3, named by genes according to functional candidacy; Table 1, Supplementary Table ⁶, Supplementary Data ¹). Among the five novel loci, four were not detectable using the conventional univariate analysis. The IGH locus could be identified using univariate GWAS on GP9, however, it was not mapped in the previous univariate GWAS with even larger discovery sample²³ (N=2247).

Open in a separate window

Fig. 1

Manhattan plots of multivariate GWAS for IgG N-glycosylation phenotypes in the ORCADES discovery population. The known and novel loci are labeled in black and red, respectively. 23 IgG N-glycosylation phenotypes were analyzed together and also in eight different functional subgroups, including sialylation and galactosylation. a Analysis of 23 N-glycosylation traits; b Analysis of eight sialylation phenotypes; c Analysis of 17 galactosylation phenotypes. The horizontal dashed lines represent the genome-wide significant P-value threshold of 5×10⁻⁸/9=5.6×10⁻⁹ and the genome-wide suggestive significant threshold of 5×10⁻⁸

Table 1

Novel loci detected via multivariate GWAS for IgG N-glycosylation phenotypes

Candidate genes	Phenotypes grouping	Top variant	EA	Discovery: ORCADES		Replication I: KORCULA + VIS				Replication II: TWINSUK
				β S (s.e.)	P	P MANOVA	β S (s.e.)	P S	r ρ	P MANOVA	β S (s.e.)	P S	r ρ
ELL2	digalactosylation	rs11135441	T	0.0336 (0.0056)	1.6E−09	4.5E−04	0.0276 (0.0058)	2.1E−06	0.60	6.9E−09	0.0183 (0.0027)	5.5E−12	0.86
	sialylation	rs11135441	T	0.0313 (0.0050)	5.0E−10	1.1E−04	0.0234 (0.0056)	2.7E−05	0.76	3.6E−09	0.0091 (0.0014)	2.7E−10	0.96
HLA-B-C	galactosylation	rs116108880	G	0.0408 (0.0066)	8.3E−10	6.2E−01	0.0188 (0.0053)	4.0E−04	0.39	3.8E−08	0.0160 (0.0031)	2.2E−07	0.74
IGH	N-glycosylation	rs35590487	T	0.0671 (0.0080)	6.3E−17	2.3E−02	0.0218 (0.0074)	3.1E−03	0.10	1.0E−31	0.0534 (0.0040)	5.2E−41	0.76
	monogalactosylation	rs35590487	T	0.0465 (0.0055)	3.6E−17	7.1E−02	0.0178 (0.0040)	9.5E−06	0.90	6.3E−25	0.0358 (0.0032)	1.2E−28	0.86
	galactosylation	rs35590487	T	0.0637 (0.0073)	3.4E−18	8.8E−02	0.0212 (0.0059)	3.3E−04	0.02	4.7E−30	0.0501 (0.0039)	4.1E−38	0.85
	monosialylation	rs35590487	T	0.0325 (0.0046)	9.9E−13	3.5E−01	0.0114 (0.0044)	9.9E−03	0.99	6.3E−24	0.0287 (0.0026)	1.2E−27	0.97
	sialylation	rs35590487	T	0.0379 (0.0054)	1.5E−12	3.5E−01	0.0113 (0.0050)	2.3E−02	0.36	5.9E−23	0.0287 (0.0026)	1.2E−27	0.94
	fucosylation	rs58087925	T	0.0625 (0.0070)	3.8E−19	7.2E−02	0.0242 (0.0059)	4.5E−05	0.18	3.7E−34	0.0523 (0.0038)	8.3E−43	0.93
	bisecting GlcNAc	rs8013055	A	0.0370 (0.0052)	1.0E−12	3.0E−01	0.0077 (0.0058)	1.8E−01	0.70	1.1E−15	0.0261 (0.0029)	1.9E−19	0.96
AZI1	N-glycosylation	rs9319617	C	0.0453 (0.0076)	2.5E−09	9.7E−03	0.0204 (0.0064)	1.4E−03	0.69	7.6E−07	0.0199 (0.0033)	2.2E−09	0.75
	galactosylation	rs9319617	C	0.0422 (0.0073)	2.7E−09	2.6E−03	0.0189 (0.0047)	6.2E−-05	0.79	1.2E−06	0.0176 (0.0029)	8.5E−10	0.78
	fucosylation	rs2659009	A	0.0414 (0.0063)	3.9E−11	5.1E−04	0.0168 (0.0049)	5.8E−04	0.85	6.0E−05	0.0159 (0.0027)	4.0E−09	0.89
FUT6-3	N-glycosylation	rs12019136	A	0.0574 (0.0078)	1.9E−13	3.3E−01	0.0022 (0.0070)	7.5E−01	0.01	3.3E−15	0.0251 (0.0031)	6.3E−16	0.82

Open in a separate window

Nine multivariate GWA scans were performed, including one using all the 23 phenotypes, as well as eight different subgroupings according to type of glycosylation. Replication was performed by (i) MANOVA test in the replication cohorts (P _MANOVA); and (ii) testing the association between the phenotypic score (constructed based on the coefficients estimated in the discovery cohort) and the corresponding genotype dosages (P _S, reported only for replication cohorts). β _S denotes the coefficient of regression of genotype dosage onto the phenotypic score. Consistency of effects was performed by testing the correlation of partial genotype-phenotype correlations (r _ρ) in the discovery and replication cohorts. EA, effect allele. Extended details are given in Supplementary Table ⁶.

Locus-specific pleiotropic model can vary across cohorts

Geometrically, the procedure for estimating the phenotype score projects the SNP dosage vector onto the hyperplane defined by the phenotype vectors in the discovery cohort (Supplementary Fig. ¹⁰). Based on the same hyperplane of phenotypes, an exact replication would yield a projected SNP dosage vector in the replication cohort which has a correlation of 1 with that of the discovery cohort, i.e., in Supplementary Fig. ¹⁰, θ should be nearly zero. Having a significant positive genetic effect on the phenotype score in the replication cohort only replicates the fact that |θ|<90°. Namely, the “score replication” above replicates the genetic effect on the same phenotype score but does not guarantee that all the partial correlations between the SNP dosage and the phenotypes are consistent. Therefore, for each newly discovered locus, we tested the correlation between the set of these partial correlations (ρ) in the discovery cohort to that in each replication cohort. In Table 1 and Supplementary Table ⁶, r _ρ and P _ρ represent the estimates and significance of such correlations.

The partial correlations between each genotype dosage and the phenotypes can be viewed—up to a constant—as the partial coefficients from multiple regression of the dosage on the phenotypes (shown in Supplementary Data ¹). The r _ρ values are rather high (all>0.74) when comparing ORCADES and TWINSUK, whereas when contrasting ORCADES and the Croatians, for HLA, FUT6 and some phenotype groups of IGH, such correlation estimates appear to be close to zero. As an example, for the IgG galactosylation phenotypes, we visualize the partial correlation contrast for the IGH and AZI1 loci (Fig. 2). As we standardized the partial correlations in each cohort to z-scores, the regression slope in each panel represents r _ρ. We can see that the effects of the AZI1 locus correlate well across replication cohorts, whereas for the IGH locus, the correlation is low for ORCADES vs. Croatian cohorts. From Table 1 and Supplementary Table ⁶, one can see this is not an exclusive feature of the galactosylation group, but can be observed for other trait groups as well. Because differences in experimental/sample collection procedures between the cohorts would not be locus-specific, they are unlikely to be the cause. We may speculate that the observation of different multivariate association patterns across populations has a genetic explanation, such as different LD structures between British and Croatian populations and/or presence of specific environmental factors modulating the action of the loci in question.

Open in a separate window

Fig. 2

Comparison of the estimated genotype–phenotype partial correlations between discovery and replication cohorts. The IGH and AZI1 loci for the 17 traits in the IgG galactosylation group are displayed. The partial correlations were standardized to z-scores so that the regression slope and confidence intervals (red curves) represent the amount of correlation in the effect sizes between each named cohort. The gray bars represent standard errors

Connection to immune-related tissues and disease

We used DEPICT software³⁰ to perform gene prioritization, gene set and tissue enrichment analyses. Eighteen analyses (GWAS of 9 trait sets, using either genome-wide significant SNPs or SNPs with MANOVA P<1×10⁻⁵) were run. Results are reported in Supplementary Data ². The gene sets demonstrating most significant enrichment included regulation of protein kinase activity and Endoplasmic Reticulum-nucleus signaling pathway. The top enriched tissues/cell types included B-lymphocytes, plasma cells, antibody-producing cells, plasma (see Fig. 3 and Supplemetary Fig. ¹¹–¹² for trait group-specific results).

Open in a separate window

Fig. 3

DEPICT enrichment analysis results for IgG N-glycosylation loci. (a) Tissues. (b) Cell types. Genes in the associated loci are highly expressed in cells and structures of the hemic and immune systems (with emphasis on antibody-producing cells), and to a lesser degree in the skeleton and glands of the digestive system. See Supplementary Data ² for details

To link to established knowledge and to provide insights regarding potential underlying functions of the IgG glycosylation loci, we searched association databases using PhenoScanner (http://www.phenoscanner.medschl.cam.ac.uk/) for our reported top variants (8 and 6 for the new and known loci, respectively) and obtained 2011 association records. Filtered based on FDR<5%, we identified associations between SNPs detected in this study and 17 complex diseases and disease-related traits (Fig. 4, Supplementary Table ⁵).

Open in a separate window

Fig. 4

Pleiotropic network of the novel (orange) and known (blue) loci discovered using the multivariate method. The IgG N-glycosylation phenotypes are shown in light green and complex traits/disease in red. The associations with complex traits were filtered based on an FDR cutoff of 5%. CAD, coronary heart disease; BMI, body mass index; WHR adj BMI, waist–hip ratio adjusted for BMI

Discussion

We start with discussion of the methodological implications, and continue to discuss biological aspects of our findings.

Our analyses are of methodological interest over and above the new biological insight they have provided.The fact that we discover and replicate a number of new loci using a multivariate approach is consistent with results previously demonstrated by Inouye et al.¹³ for metabolomics, and, although to a lesser extent, by Stephens¹⁰ for classical lipids. As we have demonstrated, the P-value adjustment method TATES did not show the same power as our multivariate analysis. In principle, TATES produces quite similar results to the univariate GWAS, because it does not jointly model the correlations among the phenotypes and genotypes, but instead adjusts the univariate GWAS results accounting for the correlations of univariate analyses statistics. However, phenotypic correlation and correlation between statistics are not necessarily consistent with each other, especially at a single genetic variant. As such, TATES is not an ideal method for multivariate analysis, but rather one for multiple testing correction.

To examine the discovery power of different multi-trait methods, we contrasted the power of our efficient MANOVA statistic with GEMMA (as a representitive for mixed model based methods) and TATES (as a representative for P-value correction methods) under different scenarios, where we see that, for unrelated individuals, the MANOVA statistic performs no worse than the other methods (Supplementary Fig. ⁷). Although GEMMA should perform better than our two-step approach, the difference is only expected to become substantial when the sample includes a large group of highly related individuals, which is not the case with our discovery cohort. Because the variance-covariance structure of the phenotypes is modeled explicitly, we found that MANOVA and GEMMA methods are no worse than univariate or TATES in all scenarios investigated, and in some scenarios they substantially outperform TATES (Supplementary Fig. ⁷). This highlights the importance of joint modeling of multiple omic phenotypes, which is able to identify loci that were missed in univariate GWAS.

Intuitively, the power of multivariate analysis should depend on the amount and pattern of correlations among the phenotypes, and between the phenotypes and genotypes. MANOVA in general has good power for intermediately correlated phenotypes, as it approximates Fisher’s method for independent phenotypes, while nearly identical phenotypes are not expected to contribute additional information. Beyond mathematical correlation, biological relevance of the phenotypes should be considered in practice, for example to link genetic variants to (sub)groups of omics measurements reflecting specific biological processes, such as the glycosylation groups in this study.

Certainly, increasing the number of tested phenotype groups would lead to more multiplicity of the statistical testing procedure. Here, we focus on nine particular groups and Bonferroni correct the number of GWA scans, in order to be conservative in reporting discoveries and follow the conventional rules in such genetic research. However, we would also emphasize that multiplicity is not simply a problem that reduces power, but also offers an opportunity to discover more.

Another perspective of looking at power increase via multivariate analysis is that the power is gained by reducing noise in the measurements when combining multiple correlated phenotypes. For example, in our replication procedure, the effect at a single locus is interpreted by expressing the genotype dosage as a linear combination of multiple phenotypes, which is similar to our common knowledge that summing up multiple repeated measurements of a trait would lead to the reduction of relative residual variance and therefore improve statistical power³¹. Given this, multivariate analysis can be applied, not only to omics data, but also to any correlated phenotypes in general, where power may be gained due to this statistical property.

Methods for meta-analysis are critical in current statistical genomics: they allow for pooling of results of analyses from different studies thus achieving large sample sizes, statistical power and replication. As MANOVA was applied for the multivariate GWA analyses, there was no single effect size estimate for multiple phenotypes available. Meta-analysis based on the estimated effects size would therefore not be straightforward. In a large meta-analysis, we can combine the P-values with a weight for each cohort depending on its sample size. At the same time, pooled estimates of the coefficients for constructing the phenotype score could be obtained, though the meta-analysis based on such multiple coefficients will be less powerful. Meta-analysis in a replication procedure can be easily conducted based on the one degree-of-freedom test on the phenotype score estimated in the discovery cohort. However, our results suggest that pleiotropic models may differ between cohorts, which calls for extra caution (e.g., use of correlation between partial regression coefficients) when performing multivariate replication.

Our multi-phenotype method identified five novel loci associated with the human IgG N-glycome; all of them were convincingly replicated. These loci include such positional candidates as the immunoglobulin heavy locus (IGH), IgH transcription elongation factor ELL2, human leukocyte antigen (HLA-B-C), and a locus including fucosyltransferases 6 and 3 (FUT6-FUT3); additionally, we have identified a locus on chromosome 17 flanked by AZI1 and TMEM105 (Table 1, Supplementary Figs. ²–⁶).

The novel chromosome 19 association includes the FUT3-5-6 gene cluster (Supplementary Fig. ⁶). This locus was recently shown to be associated with glycosylation of plasma proteins²¹; based on our results, this is the second known locus shared in common between total plasma protein and IgG glycosylation. Products of this gene cluster are three fucosyltransferases, enzymes that catalyze transfer of fucose from the donor guanosine-diphosphate fucose to the acceptor molecules. The main products of these genes are the Lewis x and Lewis a structures that determine Lewis blood groups. However, according to current knowledge, these enzymes catalyze fucosylation of antennary GlcNAc³², resulting in glycan structures that are not found on IgG, and therefore the mechanism through which this locus affects IgG glycosylation is not obvious.

The newly established locus on chromosome 6 includes the human leukocyte antigen class I HLA-C and HLA-B genes. The HLA super-locus on chromosome 6 is a gene-rich region that has been associated with more than a hundred, mostly autoimmune, diseases. There are at least 132 protein-coding genes, as well as the classical HLA genes. Many products of these genes are molecules involved in either innate immunity or the adaptive immune response—such as the classical HLA which encode the antigen-presentation apparatus³³.

Several potential candidates are found at the novel chromosome 17 locus. We note that this locus was suggestively associated in the previous univariate analysis of Lauc et al.²³, which included our discovery cohort ORCADES as one of the five cohorts analyzed. The peak of association is flanked by the AZI1 and TMEM105 genes, and includes several other genes, among which SLC38A10 was prioritized by DEPICT. Potentially of greater interest might be IKZF3, which encodes a transcription factor that interacts with IKZF1. The IKZF1 locus was associated with IgG glycosylation in the previous GWAS and suggestively rediscovered in this study. These two transcription factors are involved in regulation of differentiation and proliferation of B lymphocytes³⁴, the cells where immunoglobulins are synthesized.

For the IGH locus (Supplementary Fig. ⁴), no gene was prioritized by DEPICT. The locus contains genes encoding the heavy chains of immunoglobulins, also including immunoglobulin G (IGHG genes). Immunoglobulin G consists of two biologically different regions—the antigen binding fragment (Fab) and crystallizable fragment (Fc). While the Fab region is responsible for binding with antigens, the Fc region is responsible for binding with effector molecules and cells³⁵, guiding the immune response. Both regions can be glycosylated, with the majority of glycans coming from the Fc region and leading the immune response³⁶. While a biological relationship between IGHG and IgG is obvious, it is not immediately clear what mechanism could link variation in IGHG region and IgG glycosylation. Similarly, it is interesting that in the novel associated interval on chromosome 5, the ELL2 gene is prioritized by DEPICT. ELL2 encodes the RNA polymerase II transcription elongation factor, which plays a role in immunoglobulin secretion. ELL2 regulates exon skipping of IGH and is necessary for processing mRNA transcribed from IGH ³⁷.

In a previous GWAS, Lauc et al.²³ used single-trait analyses to detect nine loci, of which four contained genes encoding glycosyltransferases with obvious links to IgG N-glycosylation. Here, we find an additional five loci, of which only one contains a gene directly involved in protein glycosylation. Our results show that genetic control of IgG glycosylation is a complex process involving multiple biological pathways. Another interesting observation is that we observe clear biological links between some of the positional candidate genes, for example IGHG and ELL2, IKZF1 and IKZF3. With future larger studies, we should be able to further illuminate the complexity of the genetic control of glycosylation.

Methods

Electronic supplementary material

Acknowledgements

Author contributions

Initiated and coordinated the study: X.S. and Y.S.A. Developed statistical methods: X.S. and Y.S.A. Analyzed the data: X.S. Investigated and interpreted discoveries: X.S., L.K., S.S. and D.W. Contributed to replication analysis: M.M. Contributed to simulation studies: Z.N. Contributed discovery data: J.F.W. Contributed replication data: I.R., O.P., C.H. and T.D.S. Contributed IgG glycomics data: I.T.-A., M.P.-B., G.L. Contributed to writing: X.S., L.K., M.M., T.D.S., J.F.W., G.L. and Y.S.A.

Notes

Footnotes

Contributor Information

Xia Shen, Email: [email protected].

Yurii S. Aulchenko, Email: [email protected].

References