{'fraction_non_alphanumeric': 0.03706616851833645, 'fraction_numerical': 0.025218759217382754, 'mean_word_length': 5.213500305436774, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 1, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of Potential Pleiotropic Genes for Immune and Skeletal Diseases Using Multivariate MetaCCA Analysis

Background Immune and skeletal systems physiologically and pathologically interact with each other. Immune and skeletal diseases may share potential pleiotropic genetics factors, but the shared specific genes are largely unknown. Objective This study aimed to investigate the overlapping genetic factors between multiple diseases (including rheumatoid arthritis (RA), psoriasis, osteoporosis, osteoarthritis, sarcopenia, and fracture). Methods The canonical correlation analysis (metaCCA) approach was used to identify the shared genes for six diseases by integrating genome-wide association study (GWAS)-derived summary statistics. The versatile Gene-based Association Study (VEGAS2) method was further applied to refine and validate the putative pleiotropic genes identified by metaCCA. Results About 157 (p<8.19E-6), 319 (p<3.90E-6), and 77 (p<9.72E-6) potential pleiotropic genes were identified shared by two immune diseases, four skeletal diseases, and all of the six diseases, respectively. The top three significant putative pleiotropic genes shared by both immune and skeletal diseases, including HLA-B, TSBP1, and TSBP1-AS1 (p<E-300), were located in the major histocompatibility complex (MHC) region. Nineteen of 77 putative pleiotropic genes identified by metaCCA analysis were associated with at least one disease in the VEGAS2 analysis. Specifically, the majority (18) of these 19 putative validated pleiotropic genes were associated with RA. Conclusion The metaCCA method identified some pleiotropic genes shared by the immune and skeletal diseases. These findings help to improve our understanding of the shared genetic mechanisms and signaling pathways underlying immune and skeletal diseases.


INTRODUCTION
Bone formation and bone resorption are dynamic processes occurring in bone at all times, and its imbalance contributes to bone metabolic disorders, such as osteoporosis (OP) [1]. The immune system is one of the major systems in determining the balance of bone turnover [1]. The concept of osteoimmunology has recently been proposed to summarize novel insights into functional interdependence between immune and skeletal systems at the anatomical, vascular, cellular, and molecular levels [2]. It is well established that bone and adaptive immune systems appear at the same stage of vertebrate evolution [3]. This phenomenon suggested that the immune system is required for bone during its evolution and vice versa.
Immune and skeletal systems physiologically and pathologically interact with each other. The important role of the immune system in multiple bone disease pathologies, such as OP, osteoarthritis (OA), and rheumatoid arthritis (RA), is now well established. Osteoporosis is a disorder of decreased bone mass, microarchitectural deterioration, and fragility fractures (F) [4]. OP is a prevalent inflammatory bone loss condition in elderly people. RA is a typical chronic inflammatory joint disease characterized by persistent inflammation in multiple joints, progressive bone erosions, and loss of function. The formation of bone erosions and structural damage in inflamed joints in RA were partially attributed to osteoclasts, which were generally regarded as the dominant source of both receptor activators of NF-κB (RANK) ligand (RANKL) and osteoprotegerin (OPG) [5]. The RANKL to OPG ratio is a key determinant of osteoclast differentiation and bone resorption. Psoriasis (PsO) was associated with higher rates of pathological fractures, particularly of the vertebrae, pelvis, femur, tibia, and fibula [6,7]. Previous studies also reported that skeletal muscle was a potent regulator of immune system function [8]. The immune system has been reported as a contributor to sarcopenia (SP), which is a syndrome of low muscle mass with low muscle strength and/or low physical performance.
Immune and skeletal systems are strongly interlinked through a quantity of shared regulatory factors, including cytokines, chemokines, transcription factors, and receptors [3]. Although previous findings have broadened our understanding of the interplay between these two systems, the shared molecular mechanism and signaling pathways between immune and skeletal diseases are still much less explored. Pleiotropy is a phenomenon in one gene that affects multiple phenotypes, leading to genetic correlation among phenotypes [9]. Identifying potential pleiotropic genes among phenotypes can help investigate the overlapping etiology of multiple diseases. The genome-wide association study (G-WAS)-derived summary statistic data can be utilized to identify the shared genes for immune and skeletal diseases with canonical correlation analysis (metaCCA) framework, which was able to test the association between all SNPs located at the same gene and multiple phenotypes [10]. This method incorporates the GWAS summary statistics, linkage disequilibrium (LD) structure, and phenotypic correlations between phenotypes to enhance the statistical power for identifying novel genetic associations [11]. MetaCCA extends the statistical technique of canonical correlation analysis to the setting where original individual-level records are not available and employs a covariance shrinkage algorithm to achieve robustness [10].
This study adopted the metaCCA method to assess the overlapping genetic correlation among two immune diseases (RA and PsO), four bone skeletal disorders (OP, OA, F, and SP), and all of these six diseases using the GWAS summary statistics data, respectively. The versatile Gene-based Association Study (VEGAS2) approach was further applied to validate the putative pleiotropic genes identified by metaCCA. Our findings may provide clues for understanding the shared pathogenesis of immune and skeletal systems and may serve as a reference for further genetic research or drug development.

GWAS Datasets
The summary statistics of six GWAS datasets were downloaded from GWAS Catalog (available at https://www.ebi.ac.uk/gwas/). The RA GWAS dataset comprising of ~10 million SNPs is derived from a GWAS meta-analysis with 14,361 cases and 43,923 controls of European subjects [12]. The dataset for OA GWAS included ~16 million SNPs and was downloaded from UK Biobank involving 32,970 European subjects [13]. The PsO GWAS summary data were obtained from a meta-analysis with 10,588 cases and 22,806 controls of European subjects [14]. The FNK-B-MD dataset contained the association results for approximately 10 million SNPs and were obtained from a GWAS meta-analysis with 32,965 European subjects published by the Genetic Factors for Osteoporosis (GEFOS) Consortium [15]. The summary statistics for fracture (F) included ~13 million SNPs and were based on a GWAS from a European cohort consisting of 53,184 cases and 373,611 controls [16]. Many fracture sites were included but not limited to ankle, leg, hip, spine, and arm. The SP-related traits (appendicular lean body mass) GWAS dataset was estimated in 28,330 European subjects from a subset of 15 cohorts, containing summary statistics of more than 2.3 million SNPs [17]. Basic information for the GWAS studies was summarized in Table  S1.

Data Preparation
The SNPs in the six summary statistics were pruned before multivariate analysis using the genotype data from the 1000 Genomes project as a reference [18]. Briefly, six datasets were combined to select overlapping SNP; then, an LD-based algorithm was adopted to prune SNPs with high pairwise correlations. Default values of the PLINK 1.9 software (50, 5, and 0.2) were set as parameters when calculating LD values (r 2 ) between each SNP pair. The LD was calculated for windows that contained 50 SNPs. The SNP with a lower frequency of the minor allele was excluded for each pair with r 2 > 0.2. The calculation window was then shifted forward by 5 SNPs. Then the above process was repeated until each SNP pair was in low LD.
The SNPs were annotated to their corresponding genes using the Genome Reference Consortium Human genome build 37 (GRCh37) (available at https://www.gencodegenes.org/human/releases.html.), which contained 32,801 genes. An SNP was mapped to a certain gene when it was located between the transcription starting and ending point of that gene [11].
The regression coefficient β and corresponding standard error (SE) of each SNP originatingin GWAS summary statistics of the six different phenotypes were obtained. The regression coefficient β was normalized according to: (1) 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 of each trait.

MetaCCA Analysis
MetaCCA analysis was performed to identify the putative pleiotropic genes among the six phenotypes. Detailed principles and formulas of this approach were described elsewhere [10]. Conducting metaCCA analysis requires a crosscovariance matrix between all genotypic SNPs (Ʃ XX ), a phenotypic correlation structure between phenotypes (Ʃ YY ), and the GWAS summary statistics (Ʃ XY ). The full covariance matrix (Ʃ) can be calculated according to: (2) The full covariance matrix (Ʃ) was applied to obtain the final genotype-phenotype association results, which contained a significant test for each gene in the identified canonical correlation. The p-value was adjusted using the Bonferroni approach. An adjusted p-value < 0.05 was defined as the cutoff value for significant correlation. The metaCCA test was based on the script of the R package "metaCCA".

Versatile Gene-based Association Study (VEGAS) Approach
Gene-based association analysis was performed using the Versatile Gene-based Association Study-2 (VEGAS-2) method to refine and validate the pleiotropic genes identified by metaCCA. VEGAS-2 is a gene-based approach that assesses the correlation between a phenotype and multiple SNPs within a gene and accounts for linkage disequilibrium between those SNPs to estimate the association between each SNP and each phenotype individually [19]. This method has shown higher sensitivity and lower false-positive rates compared to other gene-based approaches. Gene with a correlation p-value < 1.56E-5 (=0.05/3,209) was significantly correlated with the diseases in the VEGA-2 analysis.

Functional Annotation
For the pleiotropic genes validated by VEGAS-2, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEEG) pathway analysis using online tool DAVID (available at http://david.abcc.ncifcrf.gov/). The GO analysis contained biological processes, cellular components, and molecular function. Multiple testing corrected by FDR with a significant threshold at 0.05 was employed to assess the degree of enrichment in each GO term or KEGG pathway.
Protein-protein interaction (PPI) analysis for the genes of interest was performed using the online tool STRING (available at https://string-db.org/).

Potential Pleiotropic Genes for Six Diseases
To investigate the overlapping genetic genes shared by the six diseases, we applied metaCCA analysis using GWAS summary statistics data. After gene annotation and SNP pruning, 17,284 SNPs located in 5,145 genes were used in metaCCA analysis. The genes with correlation p-values < 9.72E-6 (= 0.05/5,145) were significantly correlated with the six diseases. About 77 potential pleiotropic genes were significantly associated with multiple diseases (p < 9.72E-6) ( Table 1). The canonical correlation r between potential pleiotropic genes and diseases ranged from 0.0363 to 0.3567. Specifically, among the 77 genes we identified, 57 were from chromosome 6.

Potential Pleiotropic Genes for Two Immune Diseases
MetaCCA analysis found the pleiotropic effects shared by two immune diseases (RA and PsO). After gene annotation and SNP pruning, 19,290 SNPs located in 6,103 genes were used in metaCCA analysis. The genes with correlation p-values < 8.19E-6 (=0.05/6,103) were significantly correlated with the two autoimmune diseases. A total of 157 genes with significant threshold (p < 8.19E-6) were identified as putative pleiotropic genes ( Table S2). The canonical correlation r between potential pleiotropic genes and diseases ranged from 0.0302 to 0.4769.

Potential Pleiotropic Genes for Four Skeletal Diseases
MetaCCA analysis was performed to identify potential pleiotropic genes for four skeletal traits (OP, OA, SP, and F). After gene annotation and SNP pruning, 67,362 SNPs located in 12,834 genes remained for metaCCA analysis. The genes with a significant threshold (p-values < 3.90E-6 (=0.05/12,833)) were identified as potential pleiotropic genes for the four skeletal diseases. We identified 319 putative pleiotropic genes that were significantly associated with multiple diseases (p < 3.90E-6) ( Table S3). The canonical correlation r between potential pleiotropic gene and phenotype ranged from 0.0438 to 0.9455. These pleiotropic genes were located dispersedly at all human chromosomes.
KEGG pathway analysis showed that the 77 pleiotropic genes were mainly enriched in pathways related to autoimmune/auto-inflammatory diseases and infection ( Table 4).

DISCUSSION
The present study performed multivariable analytical me-taCCA by combining six available independent GWAS summary statistics datasets to identify the genes shared by immune and skeletal diseases. To refine and confirm the genes identified in metaCCA analysis, the gene-based VEGAS-2 analysis was conducted for these six diseases, respectively. We identified 77 potential pleiotropic genes in the metaCCA analysis. Among these putative pleiotropic genes, 19 were associated with at least one disease in the VEGAS-2 analysis. In particular, 14 of these 19 potential pleiotropic genes were previously reported to be associated with more than one of these six diseases. By trimming the multiple testing burden at the initial filtering step, the metaCCA approach allows for the detection of more phenotype-related loci than what would be identified by VEGAS-2 alone. The identification of putative pleiotropic genes and the associated biological pathways may provide a better understanding of the shared genetic factors involved in the development of immune and skeletal diseases.
Among the 19 confirmed pleiotropic genes, 18 genes (PLCL2,DAP,COG6,ARID5B,IRF5,SLC44A4,CCR6,UQCC2,MAGI3,DCLRE1B,SYNGAP1,AFF3,RSBN1,ANKRD55,TAP2,C2) were associated with RA in the VEGAS-2 analysis. Notably, 13 of these genes were involved in the pathogenesis of RA in previous studies. Several genes, such as PLCL2, COG6, RIF5, ARID5B, and CCR6, were reported to be the susceptibility loci for RA [20][21][22]. Chemokine (C-C motif) receptor 6 gene (CCR6) encodes an important protein that is expressed in memory T-cells and immature dendritic cells and plays a critical role in B-cell maturation and differentiation [23]. A variant of the CCR6 rs3093024 variant was significantly associated with RA-risk both in a Pakistani and Chinese population [24]. Univariate SNP-multivariate phenotype analysis in the MetaCCA method indicated that rs3093026 was significantly associated with the six diseases (p=5.91E-08) ( Table  S5). Interestingly, the SNPs rs3093024 and rs3093026 were in strong LD (r 2 =0.7). This is probably because the real associated SNP was removed during the process of data preparation using an LD-based pruning method, which was performed to remove one SNP of pairs with an R 2 value greater than 0.2. Removing the real associated SNPs in the LDbased pruning method may also prevent the pleiotropic genes that are truly associated with multiple diseases from being detected. We noticed that cytokines such as RANKL and TNF-α, which play an important role both in immune and skeletal systems, were insignificantly associated with these six diseases in the metaCCA analysis. This may also result from the absence of SNPs in the gene region. As expected, there was no SNP mapped to the gene regions of RAN-KL and TNF-α.
As a chemokine receptor, CCR6 is not only the susceptibility gene for RA, through binding to a specific chemokine, CCL20, to exert its function on the maturation, differentiation, and migration of immune cells, but also plays an important role in the pathogenesis of OA, PsO and OP [23]. The CCL20-CCR6 axis is upregulated in the lesional skin of human psoriasis, and psoriasiform dermatitis can be suppressed by blocking this axis using a specific antibody [25]. In addition, CCL20 chemokine could induce osteoclast differentiation. Previous studies demonstrated that the expression of RANK on CD14 + cells was positively correlated with that of CCR6. Monocytes expressing both RANK and CCR6 differentiate into osteoclasts [26]. The above findings may help drive force for understanding the osteoimmune axis.
Of the 77 putative pleiotropic genes identified by the me-taCCA approach, several (HLA-B, HLA-DMA, HLA-DMB, HLA-DOA, HLA-DOB, HLA-DPA2, HLA-DPA3, HLA-DQA1, HLA-F, C2, TAP2) were HLA genes, which played a critical role in the processing and presentation of antigens [27]. It has been well established that the HLA region was strongly associated with autoimmune disease (AID). For instance, the association of components of the HLA class II -encoded HLA-DRB1-DQA1-DQB1 haplotype has been determined with several AIDs such as RA. Molecules encoded by this region play a key role in exogenous antigen presentation to CD4+ Th cells, indicating the importance of this pathway in AID initiation and progression. However, the association between the HLA genes and skeletal diseases has been poorly studied. A previous study aimed to investigate the relationships between polymorphisms of the HLA-B gene and postmenopausal osteoporosis showed that the frequency of the HLA-B* 3501 allele was significantly higher in postmenopausal osteoporosis patients than in the control group in a Chinese Han population and likely an important risk factor for postmenopausal osteoporosis [28]. Our results hinted that HLA genes' role in skeletal diseases might not be adequately explored.
GO analyses suggested that the identified putative pleiotropic genes for two immune diseases were significantly enriched in biological processes related to "positive regulation of T cell proliferation" (GO: 0042102, FDR=2.51E-2), "immune response" (GO: 0006955, FDR=2.51E-2), "interferon-gamma-mediated signaling pathway" (GO: 0060333, FDR=2.51E-2) and "antigen processing and presentation of peptide or polysaccharide antigen via MHC class II" (GO: 0002504, FDR=3.33E-2). These processes were reported to be associated with inflammatory autoimmune disorders. It suggested that the result of metaC-CA should be reliable and reasonable. The metaCCA method is an effective technique to systematically and comprehensively identify pleiotropic genes associated with multiple complex diseases. By leveraging large GWAS summary statistics of six diseases, the metaCCA approach could increase the sample size and statistical power of the study compared to the univariate GWAS analysis based on a cross-sectional population [10]. However, the proportions of variability could not be assessed due to no such data being available. Additionally, biological experimental verification of these pleiotropic genes identified in this study is required.

CONCLUSION
In summary, by applying metaCCA and VEGAS2 methods, we identified several putative pleiotropic genes associated with two immune diseases (e.g., RA and PsO) four skeletal disorders (e.g., OP, OA, F, and SP), and all of these six diseases by using GWAS summary statistics data, respectively. These findings may provide a novel molecular basis for understanding the interaction of the pathogenesis between immune diseases and skeletal diseases.

ETHICS APPROVAL AND CONSENT TO PARTICI-PATE
Not applicable.

HUMAN AND ANIMAL RIGHTS
No animals/humans were used for studies that are the basis of this research.

CONSENT FOR PUBLICATION
Not applicable.

AVAILABILITY OF DATA AND MATERIALS
Not applicable.

FUNDING
The study was supported by Natural Science Foundation of China (81872681, 31401079 and 81473046), the Science and Technology Project of Suzhou (SS202050), and a Project of the Priority Academic Program Development of Jiangsu Higher Education Institutions.