DOI: https://doi.org/10.21203/rs.3.rs-2611953/v1
Type 2 diabetes (T2D) is a critical risk factor for peripheral artery disease (PAD). However, the sex differences in genetic basis, causality, and underlying mechanisms of the two diseases are still unclear. Using sex-stratified and ethnic-based GWAS summary, we explored the genetic correlation and causal relationship between T2D and PAD in both ethnicities and sexes by linkage disequilibrium score regression, LAVA and six Mendelian Randomization approaches. We observed stronger genetic correlations between T2D and PAD in females than males in East Asians and Europeans. East Asian females exhibit higher causal effects of T2D on PAD than males. The gene-level analysis found KCNJ11 and ANK1 genes associated with the cross-trait of T2D and PAD in both sexes. Our study provides genetic evidence for the sex difference of genetic correlations and causal relationships between PAD and T2D, indicating the importance of using sex-specific strategies for monitoring PAD in T2D patients.
Type 2 diabetes (T2D) is one of the most prevalent metabolic disorder1 and accounts for about 90% of cases of diabetes worldwide2. One of the most frequent vascular complications of T2D is peripheral artery disease (PAD) 3that is an atherosclerotic occlusive disease occurring within the peripheral limb arteries. Delayed diagnosis and inaccurate evaluation of PAD disease progression may lead to poor prognostic outcomes4,5. Recently, T2D has been reported as a major risk factor for PAD6,7. Therefore, it is urgent to elucidate the genetic mechanisms of the epidemiological association between T2D and PAD.
Researchers have found that women have similar or higher PAD prevalence than men13, despite their lower risk of T2D14–17. Compared with individuals without T2D, the risk of cardiovascular mortality in women with T2D increases about threefold more than in men with T2D 18,19. The underlying reasons are unclear. The shared genetic foundations of the links between T2D and PAD have been observed previously8–11. For example, CDKN2A/B was detected repeatedly in several works8,9. A large-scale multi-ancestry meta-analysis (Ncase=228,499, Ncontorl=1,178,783) discovered 318 novel risk loci of T2D and PAD10. The latest study12 further identified new single nucleotide polymorphisms (SNP) associated with common genetic underpinnings of T2D and PAD in Europeans and East Asians. However, all these studies are unable to interpret the causes of these sex differences of PAD risks in women and men who are suffering from T2D.
Recently, sex-stratified genome-wide association studies (GWAS) are introduced. The method can evaluate and compare the genetic effects of SNPs between two sexes and provide evidence of sexually dimorphic associations20–25. Bernabeu et.al.20 conducted sex-aware analyses about the genetic architecture of 530 complex traits based on 450,000 European individuals in the UK Biobank and discovered > 100 traits with at least one SNP that may account for genetic sex heterogeneity. Since sex-aware single trait GWAS research paradigms have been fully developed, leveraging large population-based genetics studies to explore the sex-specific genetic effect on T2D and PAD is a feasible line of thought.
Importantly, several analysis tools available for GWAS have emerged in recent years. Linkage disequilibrium score regression (LDSC)26 is a powerful method to quantify the associations between traits based on polygenicity (i.e., small genetic effects from many SNPs) on a global scale. LAVA27 is an integrated framework for detecting genomic regions that drive the genetic correlation between traits, and advances the GWAS analysis to the locus-level. Mendelian Randomization (MR)28 is an approach of using genetic instrumental variables to obtain causal inferences from observational data29,30, which can be used to assess the risk factors of diseases. Summary data-based Mendelian randomization analysis (SMR)31 is a derivatization of MR but utilizes expression quantitative trait loci (eQTL) as instrumental variables instead. SMR integrates gene-level expression and methylation with the GWAS summary data to elucidate potential gene-level functional mechanisms. Multi-marker analysis of genomic annotation (MAGMA)32 is another gene-based method that helps identify the genes associated with a particular trait. These methods provide effective means to explore risk factors of disease from large-scale population.
In this study, we collected and generated sex-stratified GWAS summary statistics of T2D and PAD in East Asians and Europeans, respectively. We then applied LDSC and LAVA to estimate the sex differences of genetic correlations between T2D and PAD on both genome- and locus- level. MR was used for identifying potential causal relationships between two traits in two sexes, respectively. Subsequently, MAGMA and SMR were performed to identify the risk genes that are possible to be sex-specific mechanisms and account for the epidemiological association of T2D and PAD. As a result, we successfully revealed the sex differences of genetic correlations between T2D and PAD as well as the heterogeneous putative causal effects in two ethnic groups. And we identified a batch of candidate genes that play a role in the common and different biomechanisms of genetic relationships between two diseases in two sexes. Our study provides genetic insights into the different PAD risks in T2D males and females and suggests the importance of using sex-specific strategies for monitoring PAD status in T2D patients.
First, we leveraged the sex-specific GWAS summary of T2D and PAD to overview the number of genome-wide SNP (in 23 chromosomes) associations with two diseases (Fig. 1). SNPs with GWAS P < 5×10− 8 (T2D) or 1×10− 5 (PAD) were defined as significantly or suggestive significantly associated with the traits (Table S1). In Europeans, 33 and 20 SNPs were identified as associated with T2D males and females, respectively, while 46 and 40 SNPs were indicated as associated with PAD males and females, respectively. In East Asians, the AGEN-T2D (BBJ-T2D) GWAS has uncovered 64 (68) and 66 (34) independent SNPs associated with T2D in males and females, respectively, while BBJ-PAD has revealed 19 and 4 SNPs associated with PAD in males and females, respectively. These results provided a genome-wide scan of genetic heterogeneity induced by sex and ethnicity.
To confirm whether sex or ethnicity difference leads to dissimilarity in genetic correlations between T2D and PAD quantificationally, we then performed stratified single trait LDSC analysis and estimated the heritability (h2) of each disease (Fig. 2a, Table 1). In East Asians and Europeans, the heritability of two traits in males is close to that in females (Table S2). Mild genomic inflation was observed in the TAD datasets (all λGC > 1.1), while nearly no inflation was found in the PAD datasets (all λGC < 1.1) (Table 1).
Genetic correlations (rg) between T2D and PAD were computed using cross-trait LDSC, and significant genetic correlations were observed in all groups (the rg ranges from 0.27 ± 0.07 to 0.82 ± 0.41) (Fig. 2b, Table 1). In East Asians, the genetic correlations between T2D and PAD in females, [rg=0.60 ± 0.33 (using AGEN-T2D and BBJ-PAD GWAS), P = 6.60×10− 2 and rg=0.82 ± 0.41, P = 4.36×10− 2 (using BBJ-T2D and BBJ-PAD GWAS) ], are higher than in males, [rg=0.27 ± 0.07, P = 9.98×10− 5(using AGEN-T2D and BBJ-PAD GWAS) and rg=0.33 ± 0.07, P = 9.61×10− 6 (using BBJ-T2D and BBJ-PAD GWAS) ] (Table S2). In Europeans, we found slightly higher genetic correlations between T2D and PAD in females (rg=0.76 ± 0.26, P = 4.00×10− 3) than in males (rg=0.50 ± 0.08, P = 2.71×10− 10) (Table S2).
Local genetic correlation analysis was performed by LAVA that helps to detect genetic regions associated with T2D or PAD. LAVA divided the genome into 2,495 LD-independent regions for Europeans, among which more than a half showed significant univariate signals (local h2) in T2D (male:1,286/51.5%; female:1,288/51.6%) and PAD (male:1,406/56.4%; female:1,542/61.8%) (Fig. 2d, Table S4). Within these genomic regions, 843 (33.8%) and 943 (37.8%) are detected as simultaneously significant in T2D and PAD, respectively. We then estimated the local genetic correlations (rg ) between T2D and PAD via the bivariate tests in LAVA, and found only 30 and 7 genomic regions in European males and females reached significance after Bonferroni correction (Pajust <0.05/ [number of univariate genetic signals]). Among them, only the region (chr22:26,166,935 − 27,192,923) are shared by two sexes (Fig. 2d, Table S6).
The genomes of East Asian males and females from BBJ cohort are divided into 1,445 LD-independent regions by LAVA. Among them, 676(46.8%) and 657(45.5%) have significant (P < 0.05) local h2 in T2D, and 741(51.3%) and 876(60.6%) in PAD for females and males, respectively, which include 439 and 496 regions shared by T2D and PAD and were used for bivariate tests in LAVA (Table S3). The bivariate tests indicated 17 and 13 regions showing significant genetic correlations (rg) between T2D and PAD in East Asian males and East Asian females. None of these regions was simultaneously significant in both sexes (Fig. 2c, Table S5).
Moreover, the genetic correlations of T2D and PAD in 12 genomic regions (e.g., Chr13:71,978,413 − 73,934,089) showed reversed directions in two sexes, which is responsible for the significant differences of genetic correlation between sexes (Table 2). These results indicate that the genetic correlations between T2D and PAD in genomic regions are different in two sexes.
Since the ubiquitous genetic correlations in T2D and PAD were observed, we wonder if there are sex differences in the putative causal relationship between two traits. We applied six MR models to further assess the sex-specific causal relationship between two traits (Fig. 3, Table S3). The causal effects were quantified with odds ratios (ORs) in a liability scale. Results showed that East Asians presented evident unidirectional causal relationships in two sexes. As shown in Fig. 3a and Fig. 3b, T2D is a putative causal effect on PAD while PAD is not a significant effect on T2D in both males (T2D on PAD, OR = 1.21 ~ 1.31/1.14 ~ 1.23, P = 3.67×10− 11~5.98×10− 3/3.64×10− 7~3.26×10− 2; PAD on T2D, OR = 1.01 ~ 1.09/0.95 ~ 1.02, P = 5.81×10− 3~1/4.98×10− 1~9.71×10− 1) and females (T2D on PAD, OR = 1.22 ~ 1.52/1.22 ~ 1.48, P = 2.08×10− 6~4.06×10− 1/5.47×10− 5~5.09×10− 2; PAD on T2D, OR = 0.99 ~ 1.08/1.00 ~ 1.02, P = 5.81×10− 3~8.80×10− 1/2.41×10− 1~7.70×10− 1) (Table S7). Thus, males and females with T2D have approximately 1.18 ~ 1.25 and 1.35 ~ 1.43 times (i.e., the average OR of six MR estimates in the BBJ and AGEN cohort) risk of suffering from PAD compared to individuals without T2D, respectively. Further paired t-test on liability-scale causal effects (\({{\beta }_{\text{M}\text{R}}}_{\text{l}\text{i}\text{a}\text{b}\text{i}\text{l}\text{i}\text{t}\text{y}}\)) presented significant (P < 0.05) gaps between two sexes in the putative causal effects of T2D on PAD (in three of six MR models) in BBJ cohort (Table S8). In Europeans, we found bidirectional causal relationships in males (T2D to PAD, OR = 1.40 ~ 1.51, P = 3.52×10− 7~7.07×10− 2; PAD to T2D, OR = 1.05 ~ 1.08, P = 2.91×10− 9~1.73×10− 1), but no significant causality in females (T2D on PAD, OR = 0.89 ~ 1.21, P = 4.63×10− 2~8.77×10− 1; PAD on T2D, OR = 1.01 ~ 1.04, P = 1.13×10− 2~6.50×10− 1) (Fig. 3c, Table S7). These results validated sex and ethnicity as non-negligible variables of causal effects between T2D and PAD.
MR models might be confused by the instrumental variables that do not directly affect the outcome by exposure factors (i.e., correlated or uncorrelated horizontal pleiotropy). Thus, we used MR-Egger, CAUSE and MR-PRESSO to assess whether there is horizontal pleiotropy in the MR estimates (Table S9). In all groups, MR-Egger intercept is very close to zero and P is > 0.05, suggesting low impacts on causality induced by uncorrelated pleiotropy. According to the global tests of MR-PRESSO, our raw MR showed satisfying power and no outlier was filtered out. Besides, the CAUSE fitting model (causal model vs the sharing model) presented a significant (i.e., less than 0.05) ELPD P in the T2D-to-PAD direction (except for European females). These results indicate the robust accuracy of our MR analyses.
Using MTAG, we performed a meta-analysis of T2D GWAS and PAD GWAS and revealed 36 novel SNPs in total as the shared genetic variants associated with two traits (Table S10). In East Asians, we found 23 novel SNPs (in ABO and EML2) in males and 1 SNP (rs9350269, located on CDKAL1) in females associated with the cross-trait of T2D and PAD. In Europeans, seven SNPs (rs73632607, rs73632608, rs6640462, rs5934661, rs73631502, rs17255125, and rs73631501) in TBL1X and four SNPs (rs10757277, rs10811656, rs10757278 and rs10757279) in DMRTA1 were identified associated with the cross-trait genetics of T2D and PAD in males and females, respectively. These SNPs represent the shared genetic mechanisms between T2D and PAD.
To identify sex-specific risk genes associated with T2D and PAD, MAGMA was applied to both single trait GWAS and cross-trait GWAS. We identified numerous significant genes (i.e., MAGMA Padjust<0.05) in each sex-stratified GWAS (Table 3). The chi-square test was applied to test the number of overlapped genes between two sexes, aiming to find whether gender affects the disease mechanisms. We observed that the number of overlapped genes between two sexes is not larger than expected (Table S11). Only the number of common genes of two sexes associated with single trait T2D is significantly more than expected by chance (P = 1.79×10− 4 and 1.26×10− 2 for UKB-T2D and BBJ-T2D, respectively). Thus, the number of overlapped genes associated with the cross-trait T2D and PAD between two sexes are close to random, which may suggest the specific mechanisms for the co-occurrence of two diseases in males and females.
We additionally performed SMR analysis for genes showing significance (Padjust < 0.05) in MAGMA analysis. A gene was considered as significantly associated with the traits if the corrected P of SMR < 0.05/ [number of MAGMA significant genes]. The results are shown in Fig. 4 and Table S12. For East Asians, PABPC4, STK17B and PSMC3IP are associated with the single trait of female T2D, while TLE1, HMGCR and MYRF are associated with the single trait of male T2D. HP is associated with the single trait of male PAD, while no gene was found associated with female PAD. Cross-trait of T2D and PAD identified three novel genes, PIM3, RSRC1, and RPL14 in East Asian males. For Europeans, FARSA, HBQ1, KCNJ11 are associated with the single trait of female T2D, while CHMP4B, KIF11, MLX, ZBTB46, PSMC3IP and SLC39A10 are associated with the single trait of male T2D. BTN3A2, PRSS16 and SLC41A3 are associated with the single trait of the female PAD, while IREB2, PSMA4 and PSRC1 are associated with the single trait of male PAD. Cross-trait of T2D and PAD reveals JMY and NUDT5 as functional genes associated with cross-trait of T2D and PAD in European males. KCNJ11 was a cross-ethic and cross-sex gene discovered as associated with cross-trait of T2D and PAD (except European males). ANK1 was an ethic-specific but cross-sex gene associated with the cross-trait of T2D and PAD in East Asians. These genes may play important roles in the development of PAD in T2D patients cross ethnicities or sexes.
In this study, we utilized large-scale GWAS summary statistics of Europeans and East Asians to explore the sex difference of genetic correlation and causal relationship between T2D and PAD. We revealed the potential risk SNPs and functional genes contributing to the shared genetic etiology between T2D and PAD in two ethnic groups, and further provided evidence supporting the sex difference (or similarity) in T2D-PAD association.
To our knowledge, this is the first study to investigate the sex-differential genetic structure and causality between T2D and PAD. Observational studies identified significant phenotypic correlations between T2D and PAD 33,34, and indicated patients suffering from T2D with higher risks to developing into PAD. A subsequent study further investigated the shared genetic architecture between T2D and PAD in Europeans and East Asians12. Previous work has opened avenues to better understand the genetic heterogeneity underlying T2D and PAD. Meanwhile, we noticed that stratifying the analysis by sex in scientific research is becoming a trend35. Many studies21,23−25 have revealed sex differences in heritability and genetic correlations in many traits. Thus, our efforts to explore the sex differences between T2D and PAD are necessary for filling the field blank.
By conducting LDSC analysis, we observed in both ethnicities that females have higher genetic correlations (rg) between T2D and PAD than males. The estimated heritability h2 of T2D of PAD is similar between sexes. When we used LAVA to analysis the local genetic correlations (rg) between T2D and PAD, we found only few genomic regions showing significant rg in both sexes, indicating the sex-specific genetic correlations in two sexes. In Europeans and East Asians, we further found that males have less numbers of heritability-significant genomic regions than the females in both T2D and PAD while more genomic regions in males showing significant genetic correlations between T2D and PAD in males than in females. However, when we calculated the average h2 (P < 0.05 in single trait) and rg (Padjust< 0.05 in both traits) in two sexes, we found the h2 is lower in females than in males, but the rg is higher (FigS1). These findings may consistent our observations that females in T2D and PAD have higher genetic correlations than males in the global genetic correlation analysis.
When we applied multiple MR methods to investigate the putative causal effect of T2D on PAD in two sexes. Highly concordant results from six models confirmed that T2D has a putative causal effect on PAD in East Asian males, East Asian females and European males, while no evident effect was found in European females. More importantly, we found a stronger putative causal effect of T2D on PAD in females than in males among East Asians. However, in Europeans, we found a stronger putative causal effect of T2D on PAD in males than in females. These results embodied the complexity of causality across different populations and genders. Since the sample size of the current study was limited, especially when considering PAD patients, a more accurate conclusion of this observation is anticipated in the future.
This study confirmed that the well-known ethnic- and sex- differences in T2D and PAD prevalence36–38 are associated with genetic factors. For example, the T2D prevalence in men is slightly higher than that of women in Western European or Asian descent14–17, the PAD prevalence is similar or higher in women than in men13, and the PAD prevalence is higher in non-white women than white women39–41. Otherwise, the risk factors and clinical manifestations of T2D and PAD also differ by sex42–46. All these previous studies are consistent with our conclusion that T2D in females has higher putative causal effects on PAD than in males in East Asians. Importantly, this study also uncovered the similar and different mechanisms of occurring T2D and PAD in two sexes. We conducted MAGMA and SMR analysis to identify a batch of sex-specific and ethnicity-specific functional genes associated with single trait or cross trait of T2D or PAD. The high heterogeneity of risk genes in two sexes indicates the importance of sex-aware research. Notably, KCNJ11 is a cross-ethnicity genes relevant to T2D-PAD causality in two sexes. KCNJ11 encodes the pancreatic β-cell KATP channel and is a potential susceptibility gene for T2D in various populations47–49. ANK1 is a risk gene of T2D in East Asian males and East Asian females according to our analysis. ANK1 is a known T2D susceptibility gene which alters DNA-protein complex binding and confers impaired insulin release50,51. Here, we confirm the association of this gene in both sexes. These results remind us that KCNJ11 and ANK1 and related pathways are potential drug targets for PAD in T2D patients. Regarding to the GO enrichment analysis (FigS2), the enriched functions of candidate T2D-associated genes for both sexes in Europeans and East Asians are highly similar, suggesting common pathogenesis for T2D. Importantly, we noticed that HP (Chr16:72,088,403 − 72,094,954) showed functional association with PAD by two gene-based methods in BBJ males, which is located within the genomic region showing significant genetic heritability and genetic correlation between T2D and PAD (Chr16:69,541,877 − 73,543,801) in East Asians according to LAVA analysis (Table S3, Table S5). Thus, the heritable hot spot region/gene turned out to be functionally important for the T2D-PAD relationship. HP has previously been reported an additional independent cardiovascular disease risk factor in diabetic patients52,53, our results corroborated this point from a new perspective.
There are many limitations in our analysis inevitably. First, SNPs in the MHC region were excluded to avoid the effects caused by the complicated LD pattern within this region, which may cause underestimations of the shared genetic basis between T2D and PAD. Second, using current European- or East Asian-based PAD GWAS, we could not obtain enough (i.e., at least 10) instrumental SNPs when using 5×10− 8 as the P threshold for MR analysis. A loosen threshold was set (with P < 1×10− 5) instead, which may violate the assumption of MR assumptions. However, our results of multiple MR models and sensitivity analysis validated a slight influence, suggesting our findings are credible. Third, the lack of eQTL summary data from East Asians led to the potential insufficiency of our SMR results and limited the power to find trans-ethnic functional genes related to T2D and PAD.
In conclusion, our study provides the first evidence for the sex difference of genetic association between T2D and PAD as well as the heterogeneous putative causal effects in two ethnic groups. We found a batch of sex-specific SNPs, genes, and regions which might involve in the pathogenesis of T2D and PAD. Our study further advances the previous non-sex-stratified work and improved our understanding of the pathogenesis of PAD in T2D patients. It will not only provide a theoretical basis for a precise therapeutic design for PAD in T2D patients but also promote equity in allocating research and medical resources to both men and women.
The East Asian sex-stratified GWAS summary statistics datasets of T2D and PAD were obtained from the BioBank Japan (BBJ) and the Asian Genetic Epidemiology Network (AGEN). BBJ (http://jenger.riken.jp/en/)54 is a database that collects disease samples and GWAS data from > 200,000 registered Japanese individuals (male: female = 53.10%: 46.90%; average baseline age at 62.70 and 61.5 for men and women). The sex-stratified BBJ-T2D GWAS55 (Ncase = 40,250 [male = 25,705, female = 14,545], Ncontrol = 170,615 [male = 82,774, female = 87841]) and BBJ-PAD GWAS56 (Ncase = 3,593 [male = 2,784, female = 809], Ncontrol = 208,860 [male = 106,563, female = 102,297]) were both derived from a linear mixed model via SAIGE57, adjusted by age and top five genetic principal components. In addition, we downloaded sex-stratified GWAS summary statistics of T2D from the AGEN (https://blog.nus.edu.sg/agen/summary-statistics/t2d-2020/). The AGEN-T2D GWAS1 contains 22 East Asian population-based cohorts (Ncase = 55,397 [male = 28,027, female = 27,370], Ncontrol = 162,425 [male = 89,312, female = 135,055]) and was generated by a fixed-effect inverse-variance-weighted meta-analysis using METAL58. Considering the sample overlap between AGEN-T2D and BBJ-PAD (~ 14%) is less than that between BBJ-T2D and BBJ-PAD, in cross-trait analysis (e.g., LDSC, LAVA, and MR), we used the analysis between AGEN-T2D and BBJ-PAD as the main result and the ones between BBJ-T2D and BBJ-PAD as a sensitivity analysis.
The European sex-stratified GWAS summary statistics datasets of T2D and PAD were generated with our in-house UK Biobank (UKB) data (application number: 65805). UKB involves > 500,000 individuals of European ancestry from the UKB cohort (male: female = 45.60%: 54.40%; average baseline age at 56.74 and 56.35 for men and women). Quality control was conducted by removing: (1) The sample with non-British ethnicity and a missing rate > 0.05. (2) SNPs with call rate < 0.05, minor allele frequency (MAF) < 0.01, imputation INFO score < 0.8, Hardy-Wein-berg test P < 1× 10− 6. Finally, a set of ~ 425,000 individuals and 8.8M SNPs were included in this study. Cases of a certain disease were determined with the code of the International Classification of Diseases 10th (ICD-10) (UKB Field ID: 41270). An individual with an ICD-10 code of non-insulin-dependent diabetes mellitus (E11) or diseases of arteries, arterioles and capillaries (I70-I79) was identified as a T2D or a PAD patient1, and the corresponding control is the individual without T2D or PAD. Then, we used BOLT-LMM57 to generated sex-stratified UKB-T2D GWAS (Ncase = 22,813 [male = 14,059, female = 8,754], Ncontrol = 402,585 [male = 181,145, female = 221,440]) and UKB-PAD GWAS (Ncase = 6,294 [male = 3,712, female = 2,582], Ncontrol = 419,104 [male = 191,492, female = 227,612]) after adjusting age and cryptic relatedness. We converted the SNP effects (\(\beta\)) and their associated standard errors (\(SE\)) to a quantitative scale (\(\beta {\prime }\) and \(SE{\prime }\)) using the approximate approach \(\beta {\prime }\left(\text{o}\text{r} SE{\prime }\right)= \frac{\beta \left(\text{o}\text{r} SE\right)}{\mu (1-\mu )}\), where µ is the proportion of cases59.
We also generated the sex-stratified reference genome (including 23 chromosomes) of East Asians and Europeans. The East Asian reference genome was generated with 491 Chinese males and 831 Chinese females. The European reference genome was generated with 10,000 British males and 10,000 British females randomly sampled from the whole British population. All reference genome were hg19-based, and SNPs with MAF < 0.01 or located within the major histocompatibility complex (MHC) region (chromosome 6: 28,477,797 − 33,448,354)60 was excluded.
We used LDSC to investigate the genetic correlation between T2D and PAD in males and females across the genome. The pre-computed LD scores of two populations were drawn from 481 East Asians and 489 Europeans from the 1000 Genomes and are available at the LDSC website (https://alkesgroup.broadinstitute.org/LDSCORE/). Then, we produced the sex-stratified LD score in the X-chromosome through the in-house reference. Before applying LDSC analysis, the SNPs to input were filtered with HapMap3 data31 (for autosomes) and the reference genome (for X-chromosome) for the purpose of quality control. The SNPs which were strand-ambiguous (i.e., the A/T and G/C SNPs) were excluded.
Single trait LDSC was conducted to estimate the sex-specific liability-scale heritability (h2) of T2D and PAD on global levels. The population prevalence of T2D (male: 8.6%61, female: 6.8%61) and PAD (male: 6.3%62, female: 7.0%62) in East Asians were set according to previous studies conducted by Wu et al.61 and Wang et al.62. Similarly, the corresponding prevalence of T2D(male: 3.8%63, female: 3.1%63) and PAD (male: 6.4%64, female: 5.1%64) in Europeans were collected based upon existing research reported by The Information Centre63 and Kroger et al.64.
Then, we implemented sex-stratified cross-trait LDSC to calculate the genetic correlation (rg) between T2D and PAD with unconstrained intercept26,65. For genetic correlation (rg), a P < 0.05/2 = 0.025 was considered significant. To quantify the sex difference in rg, the Z statistic was constructed as follows: \(Z=\frac{{Z}_{male}^{{\prime }}-{Z}_{female}^{{\prime }}}{\sqrt{s}}\), where \({Z}^{{\prime }}\) is converted from rg using the Fisher's Z-transformation66 (\({Z}^{{\prime }}=0.5\text{l}\text{n}\left(\frac{1+{r}_{g}}{1-{r}_{g}}\right)\)), \(s\) is defined as the standard error for the difference66(\(s=\sqrt{\frac{1}{{n}_{male}-3}+\frac{1}{{n}_{female}-3}}\)), and \(n\) refers to the sample size of rg66 (\(n=\frac{1-{r}_{g}}{{SE}_{rg}^{2}}+2\)). Differences were deemed significant if the P of the Z statistics was less than 0.05. The Z statistic of h2 was also established.
Local Analysis of [co]Variant Annotation (LAVA)27 was performed to estimate genetic correlations on particular genomic regions. It accounts for correlated SNPs due to LD by converting marginal SNP effects into their combined effects, based on an external LD reference. The inputted GWAS data and reference genome data are the same as LDSC analysis. For locus definition, we download the Europeans’ example locus file (based on 1000g data phase 3, build GRCh37/hg19) provided by LAVA’s support data(https://github.com/josefin-werme/LAVA/tree/main/support_data)67, which includes 2,495 semi-independent blocks of ~ 1 Mb. And we used approximately LD-independent breakpoint dataset of East Asians provided by Berisa et al68, which partitions the genome into 1445 blocks of ~ 1.8 Mb (based on 1000g data phase 1, build GRCh37/hg19). Using LAVA, we detected the univariate local genetic signal within each trait (i.e., the local h2). Loci that have significantly estimated local h2 (defined with a default threshold of P < 0.05 in univariate test) in both T2D and PAD were then used for bivariate tests (two-sided). P of the local rg were corrected (i.e., Padjust = 0.05/ [number of bivariate tests]) and the threshold was set to 0.05.
Z statistic was constructed to test the sex difference of local rg on particular region using the same method mentioned above. The lower and upper bound of the 95% CI of rg were transferred into the standard error of rg using the following equation: \({SE}_{rg}=\frac{{{r}_{g}}_{upper}-{{r}_{g}}_{lower}}{2*1.96}\). The threshold of significant sex difference was defined by P < 0.05.
For each sex, we performed Mendelian Randomization (MR)28 methods to infer the causal relationship between T2D and PAD. MR28 is a technique that can be used to determine the causal effects of a risk factor (i.e., exposure) on a trait (i.e., outcome). Six models were used in our analysis to lower the weakness or intrinsic biases in a single method, including Inverse-variance-weighted (IVW)69, MR-Egger70, generalized summary-data-based Mendelian Randomization (GSMR)71, weighed median72, weighted mode73, and the causal analysis using summary effect estimates (CAUSE)74.
Of these, five methods (IVW, MR-Egger, GSMR, weighted median, and weight mode) selected instrumental SNPs survived from LD clump with GWAS P < 5×10− 8 or 1×10− 5 when T2D or PAD served as “exposure” since no significant SNP remains in PAD GWAS at the P threshold of 5×10− 8. The LD clump was performed with LD \({r}^{2}\)<0.05 within a 1000-kb window. CAUSE selected more instrumental SNPs with GWAS P < 1×10− 3 and with LD \({r}^{2}\) < 0.1. To evaluate pleiotropy, we used MR-Egger intercept (indicates the uncorrelated pleiotropy), CAUSE q (indicates the uncorrelated pleiotropy), CAUSE η (indicates the correlated pleiotropy) and ELPD P (indicates the power of causal model). Small MR-Egger intercept, CAUSE q, CAUSE η and significant ELPD P are ideal. The MR analysis mentioned above can be performed using R packages “TwoSampleMR”, “gsmr” and “cause”. If significant pleiotropy is detected (P < 0.05), outlier tests would be performed to filter SNPs and check the robustness of the estimates via R package “MR-PRESSO”. A value of Bonferroni-corrected level MR P (i.e., <0.05/12≈4.17×10−3, adjusted by six bi-directional tests) was regarded as a significant causal relationship.
Then, we converted the logit-scale causal effects (i.e., \({\beta }_{MR}\)) estimated by MR models to liability-scale based on the method proposed by Byrne et al75:
,
where \({K}_{\text{x}}\) and \({K}_{\text{y}}\) refer to the population prevalence of exposure and outcome, and \({Z}_{{K}_{\text{x}}}\) and \({Z}_{{K}_{\text{y}}}\)represent the values of the standard normal distribution at the corresponding prevalence. Then, we transformed the liability-scale \(\beta\) to odds ratios (OR). The sex difference of causal effect between two sexes estimated with liability \({\beta }_{MR}\) was evaluated by paired t-test. Significance was set at a P less than 0.05.
We explored the mechanism underlying the genetic relationship between T2D and PAD by investigating the functional gene of single trait T2D, PAD, and their cross-trait meta-analysis. Cross-trait GWAS of T2D and PAD were generated with inverse-variance-weighted meta-analysis via Multi-Trait Analysis of Genome-wide association summary statistics (MTAG)76, for investigating sex-specific risk genes shared by T2D and PAD. MTAG uses GWAS summary statistics from multiple traits and can boost statistical power when the traits are genetically correlated. Novel SNPs were ascertained with criteria of genome-wide significant (MTAG GWAS P\(<\)5\(\times\)10−8) associations with the cross-trait shared architecture of T2D and PAD but were not included in the single trait T2D or PAD.
Novel functional genes were annotated if they exhibited significant associations with the cross-trait shared architecture of T2D and PAD using MAGMA. We used 19,427 protein-coding genes (NCBI 37.3) for annotation, where SNPs were mapped to a gene if they were located within its transcription region. The 1000 Genomes Project European- and East Asian-based LD reference panels were utilized to correct LD structure. The sex-stratified gene-based analysis was performed with the SNP-wise Mean model, using GWAS summary SNP P. The significance thresholds of candidate genes were set at MAGMA Padjust <0.05. The number of significant genes overlapped in two sexes was then used for chi-square test to find whether overlapped genes are more than expected by chance. To investigate the function of candidate genes identified by MAGMA and shared by two sexes, we used R package clusterProfiler77 (https://guangchuangyu.github.io/software/clusterProfiler) to conduct GO analysis (Biological Processes).
Summary data-based Mendelian randomization analysis (SMR)31 was applied to capture sex-differential causal genes from single trait (T2D or PAD) and multi-trait GWAS. SMR uses the top SNPs in cis-eQTL (GWAS P < 5×10− 8) as instrumental variables to detect potential functional genes. Integrating GWAS summary data and eQTL information, it can enhance the power to search for significant gene-trait associations78. In our study, we accessed the public blood-based cis-eQTL summary data79, which contains 19,250 probes and covered > 30,000 Europeans, from the eQTLGen Consortium (https://eqtlgen.org/cis-eqtls.html). Then, SMR analysis was applied to each sex-stratified single trait GWAS, and the cross-trait GWAS of T2D and PAD generated with MTAG76. To correct the unexpected associations confounded by linkage (e.g., causal SNPs are in LD with each other, which alter the gene expression and the trait, respectively), we conducted HEIDI test31 in each SMR estimate. Functional genes were determined with MAGMA-significant genes with SMR FDR < 0.05/ [number of MAGMA-significant genes] and HEIDI P > 0.05 from at least 10 SNPs.
The authors thank the UKB, BBJ and AGEN project for making data available. H.Zhao designed the study. Z.L and H.Zhang conducted analyses, with assistance from H.Zhang, H.Zhao. Z.L, H.Zhang, Y.Y and H.Zhao wrote the manuscript. H.Zhao supervised the study. All authors contributed to the final revision of the paper.
The work was funded by the Natural Science Foundation of China (81801132; HY.Zhao, Sun Yat-sen Memorial Hospital; 61772566, 62041209, and U1611261)
Funding
The work was funded by the Natural Science Foundation of China (81801132, and 81971190; HY.Zhao, Sun Yat-sen Memorial Hospital; 61772566, 62041209, and U1611261)
Competing Interests
All authors state they have no conflict interests.
Author Contributions
The authors thank the UKB, BBJ and AGEN project for making data available.H.Zhao designed the study. Z.L and H.Zhang conducted analyses, with assistance from H.Zhang, H.Zhao. Z.L, H.Zhang, Y.Y and H.Zhao wrote the manuscript. H.Zhao supervised the study. All authors contributed to the final revision of the paper.
Data availability
Summary statistics are publicly available at https://blog.nus.edu.sg/agen/summary-statistics/t2d-2020/ and http://jenger.riken.jp/en/.
Codes and software used in the study are available in the following links: LDSC: https://github.com/bulik/ldsc; LAVA: https://ctg.cncr.nl/software/Lava; TwoSampleMR: https://mrcieu.github.io/TwoSampleMR/; CAUSE: https://jean997.github.io/cause/index.html; GSMR: http://cnsgenomics.com/software/gsmr/; MR-PRESSO: https://github.com/rondolab/MR-PRESSO; MAGMA: https://ctg.cncr.nl/software/magma; SMR: https://cnsgenomics.com/software/smr/; MTAG: https://github.com/JonJala/mtag;
|
Population |
Model (unconstrained intercept) |
Statistics |
T2D |
PAD |
T2D |
PAD |
|
Female |
Male |
|||||
|
European (UKB-T2D and UKB-PAD) |
Single trait LDSC |
Heritability (h2 ± SE) |
0.03 ± 0.00 |
0.00 ± 0.00 |
0.06 ± 0.00 |
0.01 ± 0.00 |
|
Ph2 |
1.50E-33 |
9.56E-02 |
1.48E-40 |
5.58E-06 |
||
|
λGC |
1.15 |
1.00 |
1.20 |
1.05 |
||
|
LDSC intercept |
1.02 |
1.00 |
1.04 |
1.00 |
||
|
Cross-trait LDSC |
Genetic correlation (rg ±SE) |
0.76 ± 0.26 |
0.50 ± 0.08 |
|||
|
Prg |
4.00E-03 |
2.71E-10 |
||||
|
East Asian (AGEN-T2D and BBJ-PAD) |
Single trait LDSC |
Heritability (h2 ± SE) |
0.08 ± 0.01 |
0.00 ± 0.00 |
0.10 ± 0.01 |
0.03 ± 0.01 |
|
Ph2 |
2.76E-23 |
2.95E-01 |
4.01E-21 |
3.16E-09 |
||
|
λGC |
1.22 |
0.99 |
1.18 |
1.06 |
||
|
LDSC intercept |
1.02 |
0.98 |
1.02 |
0.98 |
||
|
Cross-trait LDSC |
Genetic correlation (rg ±SE) |
0.60 ± 0.33 |
0.27 ± 0.07 |
|||
|
Prg |
6.60E-02 |
9.98E-05 |
||||
|
East Asian (BBJ-T2D and BBJ-PAD) |
Single trait LDSC |
Heritability (h2 ± SE) |
0.10 ± 0.01 |
0.00 ± 0.00 |
0.11 ± 0.01 |
0.03 ± 0.01 |
|
Ph2 |
1.64E-28 |
2.95E-01 |
3.38E-25 |
5.65E-10 |
||
|
λGC |
1.16 |
0.99 |
1.17 |
1.06 |
||
|
LDSC intercept |
0.99 |
0.98 |
1.01 |
0.98 |
||
|
Cross-trait LDSC |
Genetic correlation (rg ±SE) |
0.82 ± 0.41 |
0.33 ± 0.07 |
|||
|
Prg |
0.0436 |
9.6131E-06 |
||||
T2D: type 2 diabetes. PAD: peripheral artery disease. h2: heritability. rg: genetic correlation. SE: standard error. Ph2: P for estimated h2.
P rg : P for estimated rg. LDSC: linkage disequilibrium score regression. λGC: Genomic Inflation Factor.
|
Population |
Locus |
Chr |
Start |
Stop |
rg(famale) |
p(famale) |
rg(male) |
p(male) |
Z |
Psexdiff |
|---|---|---|---|---|---|---|---|---|---|---|
|
East Asian (BBJ-T2D and BBJ-PAD) |
368 |
4 |
54658565 |
57458942 |
0.92 |
0.01 |
-0.29 |
0.09 |
-1.93 |
0.03 |
|
661 |
7 |
49211848 |
51674880 |
0.45 |
0.07 |
-0.60 |
0.15 |
-1.72 |
0.04 |
|
|
1101 |
13 |
71978413 |
73934089 |
0.85 |
0.00 |
-0.48 |
0.02 |
-2.56 |
0.01 |
|
|
East Asian (AGEN-T2D and BBJ-PAD) |
232 |
2 |
238502381 |
240328718 |
0.18 |
0.43 |
-0.74 |
0.01 |
-1.87 |
0.03 |
|
269 |
3 |
59704331 |
61958750 |
0.92 |
0.00 |
-0.06 |
0.77 |
-2.26 |
0.01 |
|
|
351 |
4 |
16579695 |
18752293 |
-0.03 |
0.88 |
-0.74 |
0.07 |
-1.74 |
0.04 |
|
|
399 |
4 |
111466213 |
114687995 |
0.60 |
0.00 |
-0.50 |
0.13 |
-1.96 |
0.03 |
|
|
533 |
5 |
174167623 |
176509938 |
0.33 |
0.00 |
-0.56 |
0.08 |
-1.90 |
0.03 |
|
|
1078 |
13 |
29435976 |
30625939 |
0.39 |
0.28 |
-0.67 |
0.10 |
-1.65 |
0.05 |
|
|
1378 |
20 |
6669967 |
8602105 |
0.34 |
0.13 |
-0.62 |
0.04 |
-1.73 |
0.04 |
|
|
European (UKB-T2D and UKB-PAD) |
179 |
1 |
230227248 |
231287764 |
0.94 |
0.00 |
-0.10 |
0.78 |
-1.98 |
0.02 |
|
541 |
3 |
137372142 |
138693846 |
0.87 |
0.00 |
-0.08 |
0.81 |
-1.71 |
0.04 |
|
|
1866 |
13 |
21540734 |
22795188 |
0.85 |
0.00 |
-0.08 |
0.80 |
-1.71 |
0.04 |
T2D: type 2 diabetes; PAD: peripheral artery disease; Chr: chromosome; Start/Stop: start and end position of locus in basepairs; rg(famale)/rg(male): local genetic correlation estimate of given region in females/males; Z and Psexdiff: the Z statistics and the p-value of sex difference tests. The table lists genetic regions that showed significant sex difference on rg.
|
Trait |
Population |
Sex |
N. significant genes (MAGMA) |
N. significant genes (MAGMA and SMR) |
Gene symbol |
||
|
Male-specific |
Female-specific |
In both sexes |
|||||
|
T2D |
European |
Male |
1716 |
1264 |
387 |
5 |
ZBTB46, PSMC3IP, KIF11, SLC39A10, MLX |
|
Female |
3 |
FARSA, HBQ1, KCNJ11 |
|||||
|
East Asian (BBJ) |
Male |
1372 |
1132 |
481 |
3 |
HMGCR, TLE1, MYRF |
|
|
Female |
4 |
ANK1, PSMC3IP, PABPC4, STK17B |
|||||
|
East Asian (AGEN) |
Male |
1406 |
1479 |
499 |
1 |
KCNJ11 |
|
|
Female |
1 |
KCNJ11 |
|||||
|
PAD |
European |
Male |
906 |
926 |
52 |
3 |
IREB2, PSMA4, PSRC1 |
|
Female |
3 |
PRSS16, BTN3A2, SLC41A3 |
|||||
|
East Asian (BBJ) |
Male |
923 |
772 |
65 |
1 |
HP |
|
|
Female |
0 |
- |
|||||
|
Meta |
European |
Male |
1493 |
1265 |
252 |
3 |
NUDT5, JMY, IREB2 |
|
Female |
2 |
FARSA, KIF11 |
|||||
|
East Asian (BBJ) |
Male |
1324 |
1191 |
310 |
3 |
ANK1, TLE1, RPL14 |
|
|
Female |
1 |
ANK1 |
|||||
|
East Asian (AGEN) |
Male |
1348 |
1247 |
315 |
3 |
ANK1, PIM3, RSRC1 |
|
|
Female |
2 |
ANK1, KCNJ11 |
|||||
T2D: type 2 diabetes. PAD: peripheral artery disease. Meta: Meta-analysis of T2D GWAS and PAD GWAS
N. significant genes (MAGMA): The number of significant genes annotated by MAGMA (i.e., P(MAGMA) < 0.05)
N. significant genes (MAGMA and SMR): The number of MAGMA-significant and SMR-significant genes. (i.e., P(MAGMA) < 0.05 and P(SMR) < 0.05/ N. significant genes (MAGMA)).