Identification of 67 Pleiotropic Genes Associated With Seven Autoimmune/Autoinflammatory Diseases Using Multivariate Statistical Analysis

Although genome-wide association studies (GWAS) have a dramatic impact on susceptibility locus discovery, this univariate approach has limitations in detecting complex genotype-phenotype correlations. Multivariate analysis is essential to identify shared genetic risk factors acting through common biological mechanisms of autoimmune/autoinflammatory diseases. In this study, GWAS summary statistics, including 41,274 single nucleotide polymorphisms (SNPs) located in 11,516 gene regions, were analyzed to identify shared variants of seven autoimmune/autoinflammatory diseases using the metaCCA method. Gene-based association analysis was used to refine the pleiotropic genes. In addition, GO term enrichment analysis and protein-protein interaction network analysis were applied to explore the potential biological functions of the identified genes. A total of 4,962 SNPs (P < 1.21 × 10⁻⁶) and 1,044 pleotropic genes (P < 4.34 × 10⁻⁶) were identified by metaCCA analysis. By screening the results of gene-based P-values, we identified the existence of 27 confirmed pleiotropic genes and highlighted 40 novel pleiotropic genes that achieved statistical significance in the metaCCA analysis and were also associated with at least one autoimmune/autoinflammatory in the VEGAS2 analysis. Using the metaCCA method, we identified novel variants associated with complex diseases incorporating different GWAS datasets. Our analysis may provide insights for the development of common therapeutic approaches for autoimmune/autoinflammatory diseases based on the pleiotropic genes and common mechanisms identified.

Introduction

Autoimmune/autoinflammatory diseases are chronic conditions initiated by loss of immunological tolerance to self-antigens (1). In Europe and North America, these conditions occur with an estimated incidence of ~90 cases per 100,000 person-years and a prevalence of between 7.6 and 9.4% (2, 3). The chronic nature of such diseases has a significant impact in terms of the utilization of medical care, direct and indirect economic costs and quality of life. In addition, extensive clinical and epidemiologic observations have shown that autoimmune/autoinflammatory diseases are characterized by familial clustering of multiple diseases, epidemiological co-occurrence, overlapped autoantibody level, and the efficacy of therapies across diseases. These observations provide evidence that different autoimmune/autoinflammatory diseases share a substantial portion of their pathobiology and genetic predisposition underlies disease susceptibility (4–6).

The genetic effect of a single nucleotide polymorphism (SNP) or gene on two or more phenotypic traits can be described as pleiotropic and the outcome is genetically correlated. In general, this concept concerns across-trait architecture (7). To date, 186 statistically significant susceptibility loci have been identified in genome-wide association studies (GWAS) of autoimmune/autoinflammatory diseases (2, 8). These were confirmed in subsequent studies, with more than half of the susceptibility genes found to be shared by at least two distinct autoimmune/autoinflammatory diseases (9–12). For example, PTPN22 c.1858C>T (rs2476601) was identified as a susceptibility gene in independent GWAS across multiple autoimmune/autoinflammatory diseases (13). In addition, there is evidence that loci associated with predisposition to one disease can have effects on the risk of a second disease, although the risk alleles for the two diseases may not be the same (9). However, the evidence for specific shared risk variants is modest, and consequently the genetic mechanisms underlying the patterns of disease aggregation remain unclear. It is, therefore, important to identify shared genetic risk factors acting through common biological mechanisms and to assess the overlapping pathophysiological relationships among autoimmune/autoinflammatory diseases using effective analytical approaches.

GWAS is a standard univariate approach used to investigate and identify potentially causal or risk-conferring genetic variants associated with common human diseases at the level of individual measurements (14, 15). GWAS, especially those with large sample size, and meta-analysis of multiple studies, have a dramatic impact on susceptibility locus discovery and in addition, highlight and extend the previously observed commonality among autoimmune/autoinflammatory diseases. However, with the identification of millions of SNPs and a growing number of phenotypes, this univariate approach has been used to detect complex genotype-phenotype correlations with limited success. Furthermore, studies of statistical methods have confirmed that multivariate analysis provides higher statistical power for the detection of unexplained heritability due to the consideration of correlations not only among multiple SNPs, but also among different traits or diseases (16). In previous studies, bivariate analysis has been used to investigate genetic risk factors associated with complex traits, and multivariate analysis of this aspect of complex diseases is rare (17, 18). Therefore, multivariate analysis of the publicly available GWAS summary statistics in particular is highly essential and relevant to the identification of pleiotropic genes.

Cichonska et al. (19) recently employed meta-analysis using canonical correlation analysis (metaCCA) to allow multivariate representation of both genotypic and phenotypic variables based on the published univariate GWAS summary statistics. This new approach has been applied to identify potential pleiotropic genes associated with lipid-related measures, psychiatric disorders, and chronic diseases (19–22). In this study, the genetic pleiotropy-informed metaCCA method was used to identify shared variants and pleiotropic effects in seven autoimmune/autoinflammatory diseases: celiac disease (CEL), inflammatory bowel disease (IBD), which includes Crohn's disease (CRO) and ulcerative colitis (UC), multiple sclerosis (MS), primary biliary cirrhosis (PBC), rheumatoid arthritis (RA), systemic lupus erythematosus (SLE) and type 1 diabetes (T1D). In addition, gene-based association analysis was used to refine pleiotropic genes. GO term enrichment analysis and protein-protein interaction network analysis were applied to explore the potential biological function of the identified genes.

Materials and Methods

GWAS Datasets

All the GWAS summary statistics of seven autoimmune/autoinflammatory diseases investigated in the present study were downloaded from ImmunoBase (website: https://www.immunobase.org/), which is a web-based resource focused on the genetics and genomics of immunologically-related human diseases. The CEL data were obtained from a second-generation GWAS of 4,533 cases and 10,750 control subjects including 523,402 SNPs (23). The association summary statistics of IBD, including 9,735,446 imputed SNPs, were obtained from a meta-analysis with a total sample size of 59,957 subjects (24). The MS dataset consisted of 464,357 genotyped or imputed SNPs from a collaborative GWAS involving 9,772 cases and 17,376 controls of European descent collected by 23 research teams from 15 different countries (25). The PBC dataset including 1,134,141 SNPs were obtained from a meta-analysis (2,764 cases and 10,475 controls) and an independent cohort study (3,716 cases and 4,261 controls) (26). The RA dataset was also obtained from a GWAS meta-analysis of 5,539 autoantibody-positive RA cases and 20,169 controls of European descent, followed by replication in an independent set of 6,768 RA cases and 8,806 controls, which included a total of 8,254,863 SNPs (27). The SLE dataset comprised 7,219 cases and 15,991 controls of European ancestry, yielding a total of 7,915,250 SNPs from a new GWAS, meta-analysis of published GWAS and a replication study (28). The T1D dataset consisting of 8,781,607 SNPs was extracted from a Mendelian randomization analysis with 5,913 T1D cases and 8,828 reference samples (29). All the samples in the present study came from Northern and Western European ancestry (CEU) population. The summary statistics have been undertaken genomic control in the individual study or meta-analysis. Furthermore, the ImmunoBase was searched using a global cutoff for minor allele frequency (MAF) <99% for all datasets. Further details of the inclusion criteria and quality control methods used in the different GWAS studies are described in the original publications (23–29). It should be noted that the summary statistics of GWAS or meta-analysis were required to include P-values, regression coefficients and standard error.

Data Preparation

When dealing with the various datasets, we first combined the seven summary statistics for the common SNPs studied in all datasets. The result of 324,031 overlapping SNPs was then selected to for multivariate analysis. Second, gene annotation was completed according to the 1,000 Genome datasets (website:https://www.cog-genomics.org/static/bin/plink/glist-hg19) using PLINK1.9. Third, the linkage disequilibrium (LD)-based SNP pruning method was used to remove SNPs with large pairwise correlations. This SNP pruning method was performed using 50 SNPs as a window where the LD was calculated between each pair of SNPs. The MAF was also used as a criterion in the SNP pruning method to remove the SNPs with smaller MAFs for pairs with R² > 0.2. Following this initial removal of SNPs in high LD, each sliding window of five SNPs forward and the process was repeated until there were no pairs of SNPs with MAFs > 0.2 (30). All datasets were pruned using the HapMap 3 CEU genotypes as a reference panel (website: http://www.sanger.ac.uk/resources/downloads/human/hapmap3.html). Following this pruning procedure, 41,274 SNPs located in 11,516 gene regions remained and were included in the metaCCA analysis. Finally, the regression coefficient beta was normalized before conducting the metaCCA analysis in instance when the individual-level dataset genotype and phenotype matrices were not standardized. Standardization was achieved subsequently according to the following equation:

$\begin{array}{ll}{\beta }_{gp}^{STANDR}=\frac{1}{\sqrt{n}S{E}_{\text{gp}}} \times {\beta }_{\text{gp}}& (1)\end{array}$

where SE_gp is the standard error of β_gp, as given by the original GWAS result, g is the number of genotypic variables, p is the number of phenotypic variables, and n is the sample number for each disease.

MetaCCA Analysis

MetaCCA is an extension of the CCA method, which requires a full covariance matrix (∑), consisting of four covariance matrices, which can be obtained using the following formula:

$\begin{array}{ll}{\displaystyle \sum =}\left(\begin{array}{cc}{\displaystyle \sum ^{\wedge }}XX& {\displaystyle \sum XY}\\ {\displaystyle {\sum }_{XY}^T}& {\displaystyle \sum ^{\wedge }}YY\end{array}\right)& (2)\end{array}$

where ∑XY is a cross-covariance matrix between all genotypic and phenotypic variables, ${\displaystyle \sum ^{\wedge }}XX$ is a genotypic correlation structure between SNPs, ${\displaystyle \sum ^{\wedge }}YY$ is a phenotypic correlation structure between traits, and ${\sum }_{XY}^T$ is the transposition of ∑XY (19). ∑XY is constructed as the normalized regression coefficient βgp:

$\begin{array}{ll}{\displaystyle \sum XY}=\left(\begin{array}{cccc}{\beta }_{11}& {\beta }_{12}& \cdots & {\beta }_{1p}\\ {\beta }_{21}& {\beta }_{22}& \cdots & {\beta }_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ {\beta }_{g1}& {\beta }_{g2}& \cdots & {\beta }_{gp}\end{array}\right)& (3)\end{array}$

In metaCCA, ${\displaystyle \sum ^{\wedge }}XX$ is calculated using a reference database representing the study population, such as the 1,000 Genomes database, or other genotypic data available for the target population. More accurate results will be obtained if ${\displaystyle \sum ^{\wedge }}XX$ is estimated from the target population or a population of the same ethnicity instead of interracial populations (19). In our study, ${\displaystyle \sum ^{\wedge }}XX$ was estimated using the reference SNP dataset of the HapMap 3 CEU.

The phenotypic correlation structure ${\displaystyle \sum ^{\wedge }}YY$ was computed based on ∑XY, with each ${\displaystyle \sum ^{\wedge }}YY$ corresponding to a Pearson correlation coefficient between the vector of β estimates from p phenotypic variables across g genetic variants. It has been demonstrated that the accuracy of estimate increased with the value of g. Thus, ${\displaystyle \sum ^{\wedge }}YY$ values were calculated from summary statistics of 324,031 overlapping SNPs, even if 41,274 of them were used for subsequent analysis.

We next determined whether the full covariance matrix was positive semidefinite (PSD); if not, an iterative procedure was used to shrink the full covariance matrix until ∑ became PSD. In the next analysis, the PSD of the full covariance matrix was entered into the CCA framework to determine the final genotype-phenotype association (19), where the correlation between genotype and phenotype is defined as the canonical correlation, r (31).

In this study, two types of multivariate analysis were evaluated and the conservative Bonferroni correction method was used as the threshold for nominal significance. If the P-value of the canonical correlation r of any SNP was smaller than 1.21 × 10⁻⁶ (=0.05/41,274), it was deemed significantly associated with the seven diseases. Second, multivariate SNP-multivariate phenotype association analysis at the gene level was performed to identify any potential pleiotropic gene. Similarly, genes with a P-value of canonical correlation smaller than 4.34 × 10⁻⁶ (=0.05/11,516) were identified as potential pleiotropic genes associated with multiple autoimmune/autoinflammatory diseases.

Gene-Based Association Analysis

VEGAS2 (Versatile Gene-based Association Study−2) is a method of gene-based association analysis used to calculate the correlation of multiple SNPs in a gene region with one phenotype using original GWAS summary statistics (32). This method has been widely applied in the genetics studies and has been shown to offer higher sensitivity and lower false positive rates compared to other gene-based approaches (33). In the present study, the VEGAS2 method was combined with the metaCCA to refine the genes identified using the metaCCA by computing the gene-based P-value for each specific disease. This analysis was performed using https://vegas2.qimrberghofer.edu.au/, with a threshold of 1.00E-06.

GO Term Enrichment Analysis

In clarifying polygenic associations, it is useful to determine whether or not the implicated genetic variants occur in genes involved in a biological pathway (34). GO term enrichment analysis classifies gene functions based on three main categories: molecular function, cellular component and biological process. We conducted GO term enrichment analysis to determine which GO term were over- or under-represented). All significant genes re-identified by VEGAS2 in our study were annotated and enriched (website: http://amp.pharm.mssm.edu/Enrichr/). An adjusted P < 0.05 in the enrichment analysis was considered to indicate nominal significance (35).

Protein-Protein Interaction Network

Protein-protein interactions (PPIs) are crucial for all biological processes and the networks provide many new insights into protein function (36, 37). In order to detect interactions and associations of all putative pleiotropic genes, PPIs analysis were conducted by searching the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (website: http://string-db.org/), which comprises known and predicted associations from curated databases or high-throughput studies, in addition to other associations derived from text mining, co-expression, and protein homology (38).

Results

Potential Pleiotropic SNPs and Genes Identified by metaCCA Analysis

After gene annotation and SNP pruning, 41,274 SNPs located in 11,516 gene regions were available for the metaCCA analysis. The size of SNP representation among the genes ranged from 1 to 213 SNPs, and the median number of SNPs in each gene was 3.72. For the univariate SNP-multivariate phenotype analysis, 4,962 SNPs reached the Bonferroni corrected threshold (P < 1.21 × 10⁻⁶), and the canonical correlation r between each SNP and phenotype ranged from 0.0372 to 0.6586. The results are presented in a Manhattan plot in Figure 1. For the multivariate SNP-multivariate phenotype analysis, 1,044 genes with a significance threshold of P < 4.34 × 10⁻⁶ were identified as potential pleiotropic genes. The canonical correlation r between genotype and phenotype ranged from 0.0322 to 0.5899.

FIGURE 1

Figure 1. Manhattan plot of –log₁₀(metaCCA) values for univariate SNP-seven autoimmune/autoinflammatory diseases analysis. The red line marks the –log₁₀(metaCCA) value of 5.92 corresponding to P < 1.21 × 10⁻⁶. If the –log₁₀(metaCCA) value of a certain SNP was >5.92, this SNP was identified as a pleiotropic SNP for seven diseases.

Refining the Pleiotropic Genes by Gene-Based Association Analysis

The VEGAS2 algorithm based on original GWAS summary statistics was used to identify association of one gene with a specific disease. In the gene-based association analysis, 19 genes were identified for CEL, 111 for IBD, 16 for MS, 20 for PBC, 19 for RA, 20 for SLE, and 33 significant genes for T1D (P = 1.00E-06).

By screening the results of gene-based analysis P-values, we identified 67 putative pleiotropic genes yielding statistical significance in the metaCCA analysis and found to be associated with at least one disease in the VEGAS2 analysis. In particular, 17 genes were found to be associated with more than one disease in the original GWAS. The results of the metaCCA and VEGAS2 analysis are summarized in Table 1.

TABLE 1

Table 1. The 67 pleiotropic genes identified by the metaCCA and VEGAS2 analysis.

Specifically, 27 of these 67 putative pleiotropic genes had previously been reported to be associated with more than one of these seven diseases. Of these 27 confirmed pleiotropic genes, six genes (ADAD1, CIITA, CLEC16A, IL23R, MAGI3, and PTPN2) were associated with more than one disease in the VEGAS2 analysis of the original GWAS summary statistics. Of the 40 novel putative pleiotropic genes detected, 16 were previously reported to be associated with only one autoimmune/autoinflammatory disease. EFR3B and RBM17 were reported to be associated with T1D only in published studies but were shown to be associated with multiple diseases in the VEGAS2 analysis. The remaining 24 significant genes were implicated as candidate novel pleiotropic genes for these seven diseases. More significantly, nine genes (C1orf141, CALU, CCDC136, FGF2, LOC101927051, MAP4K4, MPZL3, PAPOLG, and SH3BP1) were associated with more than one disease in the VEGAS2 analysis, although they had never been reported to be associated with any autoimmune/autoinflammatory disease. The detailed features of 67 significant pleiotropic genes are shown in Table 2.

TABLE 2

Table 2. The features of 67 significant pleiotropic genes.

GO Term Enrichment Analysis

GO enrichment analysis revealed that the biological functions of these pleiotropic genes were involved mainly in the metabolism of lipids. When the 67 pleiotropic genes associated with autoimmune/autoinflammatory diseases were used as the gene sets for the GO term enrichment analysis, several functional terms were identified as being enriched. For the GO biological process and molecular function, there were 63 and 5 terms were identified to be significantly enriched in the development of autoimmune/autoinflammatory diseases, respectively. Details of the top five significant GO terms are shown in Table 3.

TABLE 3

Table 3. Top five significant GO term enrichment of the 67 pleiotropic genes.

Protein-Protein Interaction Network Analysis in String

The 67 putative pleiotropic genes identified were retrieved from the STRING database. Of these, 63 genes were clearly enriched in two confirmed clusters: JAK2 and MAP3K7 clusters (Figure 2). Three of the novel putative pleiotropic genes detected, FGF2, IL22RA2, and ITGAM, are involved in the JAK2 cluster, which participates in various processes such as cell growth, development, differentiation or histone modifications, and mediates essential signaling events in both innate and adaptive immunity. Three other novel genes, MAP4K4, PRKAA1, and TNIK, were involved in the MAP3K7 cluster, which acts as an essential component of the MAP kinase signal transduction pathway and plays a role in the response to environmental stress and cytokines such as TNF-alpha.

FIGURE 2

Figure 2. Protein-protein interactions between 67 pleiotropic genes associated with seven autoimmune/autoinflammatory diseases.

Discussion

In the present study, a novel analytical approach—metaCCA—was used to explore the common genetic variants associated with autoimmune/autoinflammatory diseases by combining seven available independent GWAS or meta-analysis summary statistics. A total of 67 putative pleiotropic genes were successfully identified after verification using gene-based analysis. In particular, 27 confirmed genes were identified as pleiotropic in previous different types of studies and were validated in the present study, 16 novel pleiotropic genes were previously reported to be associated with one autoimmune/autoinflammatory disease, and 24 candidate novel pleiotropic genes that had never been reported to be associated with any autoimmune/autoinflammatory disease. The improved detection of pleiotropic genes and the associated biological pathways may provide novel insights into the shared genetic factors involved in development of autoimmune/autoinflammatory diseases.

Among the 27 confirmed pleiotropic genes, 6 genes (ADAD1, CIITA, CLEC16A, IL23R, MAGI3, and PTPN2), which play an important role on the pathomechanism of autoimmune/autoinflammatory diseases, were shown to be associated with more than one autoimmune/autoinflammatory disease not only in the literature review, but also in the VEGAS2 analysis using original GWAS summary statistics. For example, common genetic variants of CLEC16A, also known as C-type lectin-like domain family 16A, had been reported to be associated with CEL, IBD, MS, PBC, and T1D (10). As the non-HLA genome-wide significant risk gene, CLEC16A is essential for autophagosomal growth and autophagy processes, which are of major importance for proper immune regulation, including regulation of inflammasome activation (39, 40). Moreover, recent data from murine studies and our PPIs analysis indicated that CLEC16A plays a key role in beta cells functions by regulating mitophagy/autophagy and mitochondrial health (41). PTPN2 is another important and confirmed pleiotropic gene associated with several autoimmune/autoinflammatory diseases we studied (10). The GO term enrichment analysis results suggested that PTPN2 encodes T-cell protein tyrosine phosphatase, acting as a negative regulator of the JAK/STAT signaling pathways downstream of cytokines and playing a prominent role in T-cell activation, signaling and/or effector function, which may represent potential targets for the pharmacotherapy of autoimmune diseases. In addition, Mei et al. (42) also showed that, in the Northeastern Chinese population, PTPN2 polymorphisms are associated with psoriasis, which is another chronic immune-mediated disease with a complex etiology.

Sixteen novel putative pleiotropic genes detected in this study had previously been confirmed to be associated with one form of autoimmune/autoinflammatory disease. Interestingly, EFR3B and RBM17 had been reported to be associated only with T1Din published studies, but were confirmed to be associated with other diseases in the VEGAS2 analysis (13, 43–45). EFR3B is an associated gene located in 2p23, which probably acts as the membrane-anchoring component and is involved in responsiveness to G-protein-coupled receptors. Although Bradfield et al. (13) have confirmed that EFR3B is an associated loci and protein-protein interaction network analysis also provided some protein information, it is still unclear whether this role is direct or indirect. RBM17, which is involved in the regulation of alternative splicing and the utilization of cryptic splice sites, is essential for survival and cell maintenance. Fortunately, genetic and serologic data suggest that the inherited altered genetic constitution located between IL2RA and RBM17 may predispose to a less destructive course of RA (46, 47). Although the 14 remaining genes were identified associated with one form of autoimmune/autoinflammatory disease in the literature review and VEGAS2 analysis, further experimental studies are required to confirm their role as pleiotropic genes associated with autoimmune/autoinflammatory diseases. CUL2, which is associated with response to the hypoxic environment and activation of tumor immunity, has been identified in association with CRO in nine independent case-control series (48). Zhang et al. (49) suggested that human immunodeficiency virus type 1 and simian immunodeficiency virus viral infectivity factor form a CRL5 E3 ubiquitin ligase complex that suppresses virus restriction by host APOBEC3 proteins, and that CUL2 eventually suppresses this pathway and increases the risk of autoimmune/autoinflammatory diseases (50).

Significantly, nine genes (C1orf141, CALU, CCDC136, FGF2, LOC101927051, MAP4K4, MPZL3, PAPOLG, and SH3BP1) were found to be associated with more than one disease in the VEGAS2 analysis, although these genes had never been reported to be associated with any autoimmune/autoinflammatory disease in previous GWAS. MAP4K4 has been enriched in several GO terms including MAP kinase c kinase activity (GO:0008349), which is an important contributor to the risk of developing type 2 diabetes mellitus in a Chinese Han population (51). In addition, Aouadi et al. (52) demonstrated that orally delivered small interfering RNA targeting macrophage MAP4K4 suppresses systemic inflammation, thus implicating this technology as a new strategy to attenuate inflammatory responses in human disease. C1orf141 is another significant candidate novel pleiotropic gene found to be associated with IBD and PBC in our study. It has been recently shown that C1orf141 is a susceptibility variant in psoriasis, a chronic inflammatory hyperproliferative cutaneous disease (53). Further studies are required to confirm our findings and provide a detailed description of each candidate novel gene and the associated pathomechanism.

Systematic and comprehensive searches for pleiotropic genes and their effects are essential for an understanding of the mechanisms underlying the development of autoimmune/autoinflammatory diseases (4). Compared to the univariate GWAS analysis based on a cross-sectional population, our study was cost-effective and reliable, not only due to the increased sample size achieved by integrating the summary statistics of seven large GWAS. This approach also increased the statistical power of the study and provided a wealth of information by the simultaneous analysis of multiple related autoimmune/autoinflammatory diseases. However, because of a lack of detailed original individual measures, we were unable to determine whether the effects of pleiotropic genes on risk of these diseases are direct or indirect. Alternative approaches and experimental studies are required to validate these novel genes identified in this study.

In summary, we have provided convincing evidence of the existence of 27 confirmed pleiotropic genes and highlighted 40 novel pleiotropic genes associated with autoimmune/autoinflammatory diseases by performing a systematic multivariate analysis of the open GWAS data using metaCCA. Furthermore, we have illustrated potential biological functions of these pleiotropic genes and our results contribute to a better understanding of common genetic mechanisms, and eventually the development of improved diagnosis, prognosis and targeted therapies.

Data Availability Statement

All datasets generated for this study are included in the article.

Author Contributions

XJ, YY, and XS conceived and designed the study, analyzed the data, and wrote the manuscript. YY and FX performed the methodology and data curation. XJ, NS, YF, JT, and YL analyzed the data. XS, WW, CS, and HD contributed to the critical revision of the manuscript.

Funding

This study was funded by the National Major Science and Technology Projects of the 13th 5-year plan of China (no. 2018ZX10715009), the National Natural Science Foundation of China (no. 81872597), the National Institutes of Health (no. AR069055, U19 AG055373, R01 MH104680, R01AR059781, and P20GM109036), and the Key Scientific and Technological Research Projects in Henan Province (no. 19A330005).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This manuscript has been released as a preprint at bioRxiv Jia et al. (54). We are grateful for the free online archive and distribution service. The language is modified by MogoEdit (http://www.mogoedit.com/).

References

1. Brooks WH. Involvement of X chromosome short arm in autoimmune diseases: comment on the article by Sharma et al. Arthr Rheumatol. (2018) 70:625–6. doi: 10.1002/art.40411

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Cooper GS, Bynum MLK, Somers EC. Recent insights in the epidemiology of autoimmune diseases: improved prevalence estimates and understanding of clustering of diseases. J Autoimmun. (2009) 33:197–207. doi: 10.1016/j.jaut.2009.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Ludwig RJ, Karen V, Frank L, Ziya K, Katja B, MS. M, et al. Mechanisms of autoantibody-induced pathology. Front Immunol. (2017) 8:603. doi: 10.3389/fimmu.2017.00603

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Cotsapas C, Hafler DA. Immune-mediated disease genetics: the shared basis of pathogenesis. Trends Immunol. (2013) 34:22–6. doi: 10.1016/j.it.2012.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Cotsapas C, Voight BF, Rossin E, Lage K, Neale BM, Wallace C, et al. Pervasive sharing of genetic effects in autoimmune disease. PLoS Genet. (2011) 7:e1002254. doi: 10.1371/journal.pgen.1002254

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Prüßmann J, Prüßmann W, Recke A, Rentzsch K, Ludwig RJ. Co-occurrence of autoantibodies in healthy blood donors. Exp Dermatol. (2014) 23:519–21. doi: 10.1111/exd.12445

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Stearns FW. One hundred years of pleiotropy: a retrospective (vol 186, pg 767, 2010). Genetics. (2011) 187:355. doi: 10.1534/genetics.110.122549

CrossRef Full Text | Google Scholar

8. Armstrong DL, Zidovetzki R, Alarcon-Riquelme ME, Tsaos BP, Criswell LA, Kimberly RP, et al. GWAS identifies novel SLE susceptibility genes and explains the association of the HLA region. Genes Immun. (2014) 15:347–54. doi: 10.1038/gene.2014.23

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Sirota M, Schaub MA, Batzoglou S, Robinson WH, Butte AJ. Autoimmune disease classification by inverse association with SNP alleles. PLoS Genet. (2009) 5:e1000792. doi: 10.1371/journal.pgen.1000792

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Li YR, Li J, Zhao SD, Bradfield JP, Mentch FD, Maggadottir SM, et al. Meta-analysis of shared genetic architecture across ten pediatric autoimmune diseases. Nat Med. (2015) 21:1018–27. doi: 10.1038/nm.3933

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Winkler TW, Day FR, Croteau-Chonka DC, Wood AR, Locke AE, Maegi R, et al. Quality control and conduct of genome-wide association meta-analyses. Nat Protoc. (2014) 9:1192–212. doi: 10.1038/nprot.2014.071

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Guerini FR, Bolognesi E, Manca S, Sotgiu S, Zanzottera M, Agliardi C, et al. Family-based transmission analysis of HLA genetic markers in Sardinian children with autistic spectrum disorders. Hum Immunol. (2009) 70:184–90. doi: 10.1016/j.humimm.2008.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Bradfield JP, Qu H-Q, Wang K, Zhang H, Sleiman PM, Kim CE, et al. A genome-wide meta-analysis of six type 1 diabetes cohorts identifies multiple associated Loci. PLoS Genet. (2011) 7:e1002293. doi: 10.1371/journal.pgen.1002293

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Tang CS, Ferreira MAR. A gene-based test of association using canonical correlation analysis. Bioinformatics. (2012) 28:845–50. doi: 10.1093/bioinformatics/bts051

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Kettunen J, Tukiainen T, Sarin A-P, Ortega-Alonso A, Tikkanen E, Lyytikainen L-P, et al. Genome-wide association study identifies multiple loci influencing human serum metabolite levels. Nat Genet. (2012) 44:269–76. doi: 10.1038/ng.1073

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Inouye M, Ripatti S, Kettunen J, Lyytikainen L-P, Oksala N, Laurila P-P, et al. Novel Loci for metabolic networks and multi-tissue expression studies reveal genes for atherosclerosis. PLoS Genet. (2012) 8:e1002907. doi: 10.1371/journal.pgen.1002907

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Evangelou E, Ioannidis JPA. Meta-analysis methods for genome-wide association studies and beyond. Nat Rev Genet. (2013) 14:379–89. doi: 10.1038/nrg3472

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Chung D, Yang C, Li C, Gelernter J, Zhao H. GPA: A statistical approach to prioritizing GWAS results by integrating pleiotropy and annotation. PLoS Genet. (2014) 10:e1004787. doi: 10.1371/journal.pgen.1004787

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Cichonska A, Rousu J, Marttinen P, Kangas AJ, Soininen P, Lehtimaki T, et al. metaCCA: summary statistics-based multivariate meta-analysis of genome-wide association studies using canonical correlation analysis. Bioinformatics. (2016) 32:1981–9. doi: 10.1093/bioinformatics/btw052

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Jia X, Yang Y, Chen Y, Cheng Z, Du Y, Xia Z, et al. Multivariate analysis of genome-wide data to identify potential pleiotropic genes for five major psychiatric disorders using MetaCCA. J Affect Disord. (2019) 242:234–43. doi: 10.1016/j.jad.2018.07.046

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Jia X, Yang Y, Chen Y, Xia Z, Zhang W, Feng Y, et al. Multivariate analysis of genome-wide data to identify potential pleiotropic genes for type 2 diabetes, obesity and coronary artery disease using MetaCCA. Int J Cardiol. (2019) 283:144–50. doi: 10.1016/j.ijcard.2018.10.102.

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Chen Y-C, Xu C, Zhang J-G, Zeng C-P, Wang X-F, Zhou R, et al. Multivariate analysis of genomics data to identify potential pleiotropic genes for type 2 diabetes, obesity and dyslipidemia using Meta-CCA and gene-based approach. PLoS ONE. (2018) 13:e0201173. doi: 10.1371/journal.pone.0201173

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Dubois PCA, Trynka G, Franke L, Hunt KA, Romanos J, Curtotti A, et al. Multiple common variants for celiac disease influencing immune gene expression. Nat Genet. (2010) 42:295–302. doi: 10.1038/ng.543

PubMed Abstract | CrossRef Full Text | Google Scholar

24. de lange KM, Moutsianas L, Lee JC, Lamb CA, Luo Y, Kennedy NA, et al. Genome-wide association study implicates immune activation of multiple integrin genes in inflammatory bowel disease. Nat Genet. (2017) 49:256–61. doi: 10.1038/ng.3760

CrossRef Full Text | Google Scholar

25. Sawcer S, Hellenthal G, Pirinen M, Spencer CCA, Patsopoulos NA, Moutsianas L, et al. Genetic risk and a primary role for cell-mediated immune mechanisms in multiple sclerosis. Nature. (2011) 476:214–9. doi: 10.1038/nature10251

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Cordell HJ, Han Y, Mells GF, Li Y, Hirschfield GM, Greene CS, et al. International genome-wide meta-analysis identifies new primary biliary cirrhosis risk loci and targetable pathogenic pathways. Nat Commun. (2015) 6:8019. doi: 10.1038/ncomms9019

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Stahl EA, Raychaudhuri S, Remmers EF, Xie G, Eyre S, Thomson BP, et al. Genome-wide association study meta-analysis identifies seven new rheumatoid arthritis risk loci. Nat Genet. (2010) 42:508–14. doi: 10.1038/ng.582

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bentham J, Morris DL, Graham DSC, Pinder CL, Tombleson P, Behrens TW, et al. Genetic association analyses implicate aberrant regulation of innate and adaptive immunity genes in the pathogenesis of systemic lupus erythematosus. Nat Genet. (2015) 47:1457–64. doi: 10.1038/ng.3434

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Censin JC, Nowak C, Cooper N, Bergsten P, Todd JA, Fall T. Childhood adiposity and risk of type 1 diabetes: a Mendelian randomization study. PLoS Med. (2017) 14: e1002362. doi: 10.1371/journal.pmed.1002362

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhang Q, Wu K-H, He J-Y, Zeng Y, Greenbaum J, Xia X, et al. Novel common variants associated with obesity and type 2 diabetes detected using a cFDR method. Sci Rep. (2017) 7:16397. doi: 10.1038/s41598-017-16722-6

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Seoane JA, Campbell C, Day INM, Casas JP, Gaunt TR. Canonical correlation analysis for gene-based pleiotropy discovery. PLoS Comput Biol. (2014) 10:e1003876. doi: 10.1371/journal.pcbi.1003876

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Mishra A, Macgregor S. VEGAS2: software for more flexible gene-based testing. Twin Res Hum Genet. (2015) 18:86–91. doi: 10.1017/thg.2014.79

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Wojcik GL, Kao WHL, Duggal P. Relative performance of gene- and pathway-level methods as secondary analyses for genome-wide association studies. BMC Genet. (2015) 16:34. doi: 10.1186/s12863-015-0191-2

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Lv W-Q, Zhang X, Zhang Q, He J-Y, Liu H-M, Xia X, et al. Novel common variants associated with body mass index and coronary artery disease detected using a pleiotropic cFDR method. J Mol Cell Cardiol. (2017) 112:1–7. doi: 10.1016/j.yjmcc.2017.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. (2013) 14:128. doi: 10.1186/1471-2105-14-128

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhang G, Zhang W. Protein-protein interaction network analysis of insecticide resistance molecular mechanism in Drosophila melanogaster. Arch Insect Biochem Physiol. (2019) 100:e21523. doi: 10.1002/arch.21523

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Fierabracci A, Milillo A, Locatelli F, Fruci D. The putative role of endoplasmic reticulum aminopeptidases in autoimmunity: Insights from genomic-wide association studies. Autoimmun Rev. (2012) 12:281–8. doi: 10.1016/j.autrev.2012.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. (2015) 43:D447–52. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Liu JZ, Almarri MA, Gaffney DJ, Mells GF, Jostins L, Cordell HJ, et al. Dense fine-mapping study identifies new susceptibility loci for primary biliary cirrhosis. Nat Genet. (2012) 44:1137–41. doi: 10.1038/ng.2395

CrossRef Full Text | Google Scholar

40. Berge T, Leikfoss IS, Harbo HF. From identification to characterization of the multiple sclerosis susceptibility gene CLEC16A. Int J Mol Sci. (2013) 14:4476–97. doi: 10.3390/ijms14034476

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Li J, Jorgensen SF, Maggadottir SM, Bakay M, Warnatz K, Glessner J, et al. Association of CLEC16A with human common variable immunodeficiency disorder and role in murine B cells. Nat Commun. (2015) 6:6804. doi: 10.1038/ncomms7804

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Mei Q, Liu C, Zhang X, Li Q, Jia X, Wu J, et al. Associations between PTPN2 gene polymorphisms and psoriasis in Northeastern China. Gene. (2019) 681:73–9. doi: 10.1016/j.gene.2018.09.047

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Oldstone MBA, Edelmann KH, Mcgavern DB, Cruite JT, Welch MJ. Molecular anatomy and number of antigen specific CD8 T cells required to cause type 1 diabetes. PLoS Pathog. (2012) 8:1352–62. doi: 10.1371/journal.ppat.1003044

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Honke N, Shaabani N, Zhang DE, Iliakis G, Xu HC, Häussinger D, et al. Usp18 driven enforced viral replication in dendritic cells contributes to break of immunological tolerance in autoimmune diabetes. PLoS Pathog. (2013) 9:e1003650. doi: 10.1371/journal.ppat.1003650

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Sharma A, Liu X, Hadley D, Hagopian W, Chen WM, Onengut-Gumuscu S, et al. Identification of non-HLA genes associated with development of islet autoimmunity and type 1 diabetes in the prospective TEDDY cohort. J Autoimmun. (2018) 89:90–100. doi: 10.1016/j.jaut.2017.12.008

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Knevel R, de Rooy DPC, Zhernakova A, Grondal G, Krabben A, Steinsson K, et al. Association of variants in IL2RA with progression of joint destruction in rheumatoid arthritis. Arthr Rheum. (2013) 65:1684–93. doi: 10.1002/art.37938

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Steer S, Abkevich V, Gutin A, Cordell HJ, Gendall KL, Merriman ME, et al. Genomic DNA pooling for whole-genome association scans in complex disease: empirical demonstration of efficacy in rheumatoid arthritis. Genes Immun. (2007) 8:57–68. doi: 10.1038/sj.gene.6364359

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Rivas MA, Beaudoin M, Gardet A, Stevens C, Sharma Y, Zhang CK, et al. Deep resequencing of GWAS loci identifies independent rare variants associated with inflammatory bowel disease. Nat Genet. (2011) 43:1066–73. doi: 10.1038/ng.952

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Zhang W, Wang H, Li Z, Liu X, Liu G, Harris RS, et al. Cellular requirements for Bovine immunodeficiency virus Vif-mediated inactivation of Bovine APOBEC3 proteins. J Virol. (2014) 88:12528–40. doi: 10.1128/JVI.02072-14

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Brooks JM, Long HM, Tierney RJ, Shannon-Lowe C, Leese AM, Fitzpatrick M, et al. Early T cell recognition of B cells following Epstein-Barr virus infection: identifying potential targets for prophylactic vaccination. PLoS Pathog. (2016) 12:e1005549. doi: 10.1371/journal.ppat.1005549

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Li T-T, Qiao H, Tong H-X, Zhuang T-W, Wang T-T. Association of common genetic variants in mitogen-activated Protein Kinase Kinase Kinase Kinase 4 with Type 2 Diabetes Mellitus in a Chinese Han Population. Chin Med J. (2016) 129:1179–84. doi: 10.4103/0366-6999.181969

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Aouadi M, Tesz GJ, Nicoloro SM, Wang M, Chouinard M, Soto E, et al. Orally delivered siRNA targeting macrophage Map4k4 suppresses systemic inflammation. Nature. (2009) 458:1180–4. doi: 10.1038/nature07774

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Zuo X, Sun L, Yin X, Gao J, Sheng Y, Xu J, et al. Whole-exome SNP array identifies 15 new susceptibility loci for psoriasis. Nat Commun. (2015) 6:6793. doi: 10.1038/ncomms7793

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Jia X, Shi N, Xia Z, Feng Y, Li Y, Tan J, et al. Identification of 67 pleiotropic genes for seven autoimmune diseases using multivariate statistical analysis. bioRxiv [Preprint]. doi: 10.1101/563973

CrossRef Full Text | Google Scholar