GWAS summary statistics for T2D (https://www.ebi.ac.uk/gwas/studies/GCST006867) (Xue et al., 2018) and COVID-19 (https://www.ebi.ac.uk/gwas/studies/GCST011073) (COVID-19 Host Genetics Initiative, 2020) were obtained from the GWAS Catalog. The T2D summary statistics consisted of 62,892 case subjects and 596,424 control subjects. The COVID-19 summary statistics consisted of 38,984 case subjects and 1,644,784 control subjects. COVID-19 data were obtained from samples of European ancestry and T2D data were obtained from mixed samples of European (N = 655,666) and South Asian (N = 3,650) ancestry. The T2D data were generated by meta-analysing the Diabetes Genetics Replication and meta-analysis (34,940 cases and 114,981 controls), Genetic Epidemiology Research on Ageing (6,905 cases and 46,983 controls), and the full cohort release of the United Kingdom BioBank databases (21,147 cases and 434,460 controls) by using the software METAL. Willer et al. (2010) There were 5.1 million and 8.9 million genetic variants in the T2D and COVID-19 GWAS summary statistics, respectively. The lift over of single nucleotide polymorphisms (SNPs) and conversion of SNPs into rs IDs were performed using the NCBI Genome Remapping Service and UCSC Table Browser, respectively.

We used conditional quantile-quantile (Q-Q) plots (Andreassen et al., 2013), fold enrichment plots and linkage disequilibrium score regression (Bulik-Sullivan et al., 2015) (LDSC) to evaluate the pleiotropic enrichment and genetic correlations between T2D and COVID-19 GWAS summary statistics. We followed the instructions provided on the GitHub page (https://github.com/bulik/ldsc/) and performed the analysis using Python 3.

To improve the discovery rate of genetic variants correlated with T2D and COVID-19, we computed the cFDR statistics. Andreassen et al. (2013) The cFDR method is based on an empirical Bayesian statistical framework and used the GWAS summary statistics for a trait of interest alongside those for a conditional trait to estimate the posterior probability that an SNP has no association with the primary trait, provided that the p-values for that SNP in both the primary (T2D) and conditional (COVID-19) traits are as small as, or smaller than, the observed p-value. Thus, by re-ranking the test statistics of the main phenotype based on the strength of the connection with the secondary phenotype, this approach increases the likelihood of identification of genetic variations linked with the primary characteristic.

The conjFDR statistic was used to investigate the genetic variations shared by the two phenotypes. The conjFDR statistic is an extension of the cFDR statistic and is defined as the maximum of two mutual cFDRs for a specific SNP. It estimates the posterior probability that an SNP is null for either trait or both, given that the p-values for both phenotypes are as small as, or smaller than, the individual p-values for each trait. For cFDR and conjFDR, we chose a conservative threshold of 0.05 per pairwise comparison. Manhattan plots based on the conjFDR were created to highlight the genomic positions of the common genetic loci.

To minimise possible biases due to complicated linkage disequilibrium (LD) patterns, all analyses were performed after removing SNPs from the extended major histocompatibility complex (MHC) (hg19 position chromosome 6: 25,11,9106–33,85,4733) and the 8p23.1 (hg19 location chromosome 8: 72,42,715–12,48,3982) genomic regions.

We defined the independent significant SNPs according to the Functional Mapping and Annotation (FUMA) protocol (https://fuma.ctglab.nl/). SNPs having a conjFDR <0.05 and at r ² < 0.6 with each other were considered independent significant SNPs and a fraction of the independent significant SNPS in approximate linkage disequilibrium with each other at r ² < 0.1 were considered lead SNPs. In addition, we used the default parameters from FUMA to determine the distinct genomic loci and their borders.

We used SNPnexus (Oscanoa et al., 2020) to annotate the shared SNPs into genes and identify the overrepresented pathways for the genes nearest the identified shared loci between T2D and COVID-19.

PASCAL (Pathway scoring algorithm) (Lamparter et al., 2016) was applied to the T2D and COVID-19 GWAS summary statistics separately. Individual SNPs (p < 0.05) from the summary statistics were first mapped using a 20 kb window around the 5′ and 3’ UTRs. The maximum number of SNPs per gene allowed by PASCAL was 3,000. LD information was retrieved from the 1,000 Genomes European Panel and the significant p-value thresholds for T2D and COVID-19 were 2.26 × 10^–6 (0.05/22,135 genes from the hg19 list) as the entire UCSC list (hg19) used by PASCAL to make calculations is included in this number of genes.

FunCoup v.4.0 (https://funcoup4.scilifelab.se/search/) was employed to expand the list of significant genes between T2D and COVID-19 with their interactors. The FunCoup database combines ten different types of functional couplings among genes to infer functional association networks: protein interaction (PIN), mRNA co-expression (MEX), protein co-expression (PEX), genetic interaction profile similarity (GIN), shared transcription factor binding (TFB), co-miRNA regulation by shared miRNA targeting (MIR), subcellular co-localisation (SCL), domain interactions (DOM), phylogenetic profile similarity (PHP), and quantitative mass spectrometry (QMS). Gene networks for T2D and COVID-19 were constructed with five shared genes between T2D and COVID-19. Expansion parameters for constructing the gene networks included a confidence threshold (0.8), a maximum number of 30 nodes per expansion step, and a query depth of one. The network expansion approach was used to find the strongest interactors for each query gene while ignoring the ties between common neighbours. Furthermore, for each gene network established with the associated p-values, enriched term analyses (Kyoto Encyclopaedia of Genes and Genomes (KEGG), GO biological function, and GO molecular function) were conducted. After adjusting for multiple comparisons using false discovery rate (FDR), gene network representation depicts the most important KEGG pathways based on their q-values. The node sizes represent the gene relevance across the whole network, while each KEGG pathway’s involved nodes are denoted in black.

GENE2FUNC, one of the functions in FUMA (Functional Mapping and Annotation of Genome-Wide Association Studies) (https://fuma.ctglab.nl/), was used to annotate common genes between T2D and COVID-19 and their interactors. Several GENE2FUNC functions were utilised, including a heatmap of gene expression and an enrichment analysis of differentially expressed genes (DEG). Using GTEx v8 (54 tissue types) data, a gene expression heatmap was generated. On the related heatmaps, the average of normalised expression per label (zero means across samples) was presented. TPM (Transcripts Per Million) for GTEx v8 are the expression values. Heatmaps provide normalised expression values (zero mean normalisation of log2 transformed expression), with greater relative expression of a gene indicated by a deep red label and lower relative expression of a gene indicated by a dark blue label.

To detect the pathways shared between T2D and COVID-19, four pathway enrichment analysis methods (GSA-SNP2 (Yoon et al., 2018), MAGENTA (Segrè et al., 2010), PASCAL (Lamparter et al., 2016), and i-GSEA4GWAS (Zhang et al., 2010)) were used to identify enriched pathways in each disease (Figure 1). GSA-SNP2 (Yoon et al., 2018) performs pathway-based analysis by testing the enrichment of associated genes in each pathway using Z-statistics of the random set models to assess pathways and a monotone cubic spline trend to determine SNP counts. MAGENTA (Segrè et al., 2010) performs pathway-based analysis on SNPs within a gene boundary using weighted Kolmogorov-Smirnov statistics that compare the ranks of genes within uniform distributions. PASCAL (Lamparter et al., 2016) calculates gene-based test statistics for all genes and carries forward the gene-based results to conduct a pathway-based test using chi-square statistics (including pathways with 10–200 genes only) that convert the corresponding p-value based on the pathway. The Benjamini–Hochberg procedure was used to account for multiple comparisons in the pathway-based analysis. i-GSEA4GWAS (Zhang et al., 2010) performs pathway analysis by using SNP label permutations to modify GWAS SNP p-values and rectify the genes and gene sets. It then multiplies the proportion ratio factor to the enrichment score to obtain the significant proportion-based enrichment score. Manhattan plots of the gene sets in each pathway were constructed and used to highlight the results of the association test for a given pathway.

Canonical pathways from curated gene sets (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#C2) and ontology gene sets (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#C5) in MSigDB (Subramanian et al., 2005; Liberzon et al., 2015) were used in the pathway enrichment analysis in all four pathway-based association approaches. Significant SNPs (p < 0.05) were mapped to genes if they were located within a range of 20 kb upstream or downstream of genes’ transcription start sites. In order to capture potential regulatory SNPs in a gene’s 5′ and 3′ untranslated regions and to prevent erroneous SNP-to-gene assignments brought on by wider windows, the 20-kb window offered the ideal width (Ghosh et al., 2013). Furthermore, we restricted the downstream analyses to pathways with 10–200 genes only in order to avoid testing over narrow or broad functional pathways (Wang et al., 2007) In addition, significant pathways were determined with a threshold of FDR <0.05.

The data underlying this article are available and downloaded in the GWAS Catalog (Buniello et al., 2019) at https://www.ebi.ac.uk/gwas/studies/GCST006867 (GWAS summary statistics for T2D) (Xue et al., 2018) and https://www.ebi.ac.uk/gwas/studies/GCST011073 (GWAS summary statistics for COVID-19) (COVID-19 Host Genetics Initiative, 2020) and can be accessed with GWAS Catalog study accessions GCST011073 and GCST006867 respectively.

With respect to canonical pathways, we found 84 pathways to be the most significant (FDR <0.05) among the gene sets with a KEGG antigen processing and presentation pathway (p-value = 1.78 × 10^–13, q-value = 3.45 × 10^–10) using the GSA-SNP2 method (Supplementary Table S1A). In addition, among 198 pathways in total, the most significant pathway was the cyclin D associated events in the G1 pathway (Reactome) (p-value < 0.001, q-value = 0) with i-GSEA4GWAS (Supplementary Table S1B). Among 39 pathways, the most significant one was the maturity onset diabetes of the young (KEGG) (p-value = 7.78 × 10^–8, q-value < 0.01) with PASCAL (Supplementary Table S1C). We also found the most significant pathway was the maturity onset diabetes of the young in the KEGG database (p-value = 5 × 10^–6, q-value = 0) with MAGENTA (Supplementary Table S1D).

Using the gene ontology database, the most significant gene set in the Gene Ontology Biological Process category was the insulin secretion pathway out of 130 pathways in total (p-value = 1.61 × 10^–13, q-value = 8.02 × 10^–10) with GSA-SNP2 (Supplementary Table S1A). In addition, the most significant gene set in the Gene Ontology Biological Process category was the interferon-gamma-mediated signalling pathway out of 383 pathways in total (p-value < 0.001, q-value = 0) with i-GSEA4GWAS (Supplementary Table S1A). We found the most significant gene set in the Gene Ontology Biological Process category was the insulin secretion pathway among 112 pathways in total (p-value = 0, q-value = 0) with PASCAL (Supplementary Table S1C). We did not find any significant pathways among the gene sets with MAGENTA.

In the canonical pathway database, we did not identify any significant pathways among the gene sets when using GSA-SNP2. However, we found 21 significant pathways among the gene sets, with Reactome chemokine receptors binding a chemokines pathway (p-value < 0.001, q-value = 0) with i-GSEA4GWAS (Supplementary Table S2A). Among the gene sets, 38 pathways were found to be significantly (FDR <0.05) associated with COVID-19 and the Reactome linked glycosylation pathway was the most significant (p-value = 0, q-value = 0) with PASCAL (Supplementary Table S2B). We did not find any significant pathways among the gene sets with MAGENTA.

With respect to gene ontology pathways, we did not find any significant pathways among the gene sets with GSA-SNP2. However, 23 pathways were significant among the gene sets, with the Gene Ontology Molecular Function G-protein coupled chemoattractant receptor activity pathway being the most significant (p-value < 0.001, q-value < 0.001), with i-GSEA4GWAS (Supplementary Table S2A). We found 151 pathways had a significant association with COVID-19 in the Gene Ontology Biological Process category, where the synapse assembly pathway was the most significant (p-value = 0, q-value = 0), with PASCAL (Supplementary Table S2A). We did not find any significant pathways among the gene sets when using MAGENTA.

To determine the common pathways between T2D and COVID-19, we compared the significant pathways that were shared between them. We found four pathways (Supplementary Table S3A) in common between T2D and COVID-19 using iGSEA4GWAS: the chemokine binding, G-protein coupled chemoattractant receptor activity pathways (CCR2 and CCR3), the TFAP2 family pathway (TFAP2B), and the ventricular cardiac muscle cells differentiation pathway (RARB and PROX1) in the canonical pathway and gene ontology database. The chemokine binding pathway plays a role in recruiting immune cells to infection sites, and G-protein coupled chemoattractant receptor functions in mediating leukocytes’ chemotaxis and promoting innate and adaptive host immune responses. Furthermore, we found 15 pathways (Supplementary Table S3B) shared between T2D and COVID-19 using PASCAL, which were associated with various organs and biological processes, such as the heart, axons and calcium channels, in the canonical pathway and gene ontology database. However, no overlapping pathways were found between the four pathway-based analysis software. Five genes (CCR2, CCR3, TFAP2B, RARB and PROX1) (Figure 2) were used to construct the gene expression heatmaps with GTEx v8 (representing 54 tissues) to investigate their expression in all tissue types.

Genome-wide LD score regression analyses showed significant positive genetic correlations between T2D and COVID-19 (r for genetic = 0.15, p-value = 0.01). We observed an enrichment of associations with COVID-19 across different levels of association with T2D (Figure 3A), indicating a small polygenic overlap between COVID-19 and T2D. We also constructed reverse conditional Q-Q plots (Figure 3B) for T2D conditional upon different levels of association with COVID-19. At a threshold of cFDR <0.05, we identified 471 loci associated with T2D conditional upon COVID-19 (Figure 4A). The reverse cFDR analysis revealed 17 loci associated with COVID-19 conditional upon T2D (Figure 4B).

Based on a threshold of conjFDR <0.05, we identified two loci shared between COVID-19 and T2D: ABO (rs505922, intronic) and NUS1 (rs3924604, intronic) (Figure 5). By comparing the effect directions of the shared independent loci (conjFDR <0.05), both these independent loci were showing consistent direction. One (rs505922) was showing positive effect while another one (rs3924604) was showing negative effect. These two genes were identified through two distinct SNPs using SNPnexus. (Oscanoa et al., 2020). We found eight pathways (Supplementary Table S4) to be significantly overrepresented among the genes nearest the identified loci shared between COVID-19 and T2D, with the defective DHDDs causing retinitis pigmentosa 59 pathway (p-value = 1.88 × 10^–4) being the most significant.

PASCAL analysis revealed an association of 394 genes in T2D (Supplementary Table S5A) and an association of 58 genes in COVID-19 (Supplementary Table S5B). After comparing the significant genes in T2D and COVID-19, five genes (Supplementary Table S5C) were found to be shared between T2D and COVID-19: PTPRD, CSMD1, MAGI1, ASIC2, DAB1.

FunCoup has detected several interactors for T2D and COVID associated genes identified by PASCAL. Gene network for shared genes between T2D and COVID-19 was constructed (Figure 6) after including 35 genes (30 subnetwork genes plus 5 query genes) (Supplementary Table S6A) and considering 92 links between them. Enrichment analysis for KEGG and GO terms (Supplementary Table S6B) has shown association of different biological processes, including endocytosis (q-value = 4.58 × 10^–3), central nervous system development (q-value = 3.64 × 10^–9), and protein domain specific binding (q-value = 5 × 10^–9) etc.

Gene expression heatmap based on GTEx V8 RNA-Seq data for T2D and COVID-19 associated genes togethers with FunCoup interactors (5 genes +30 subnetwork query genes) (Figure 7) was constructed. Specifically, six genes (ATP2A2, APP, ATN1, CTNNB1, ACTN4 and HSPA8) display increased expression levels in all available tissues compared with other genes included in the analysis. Moreover, most of the genes display the trend showing low or moderate relative expression levels on brain tissues. Differential expression gene analysis (DEG) (Figure 8) showed that all brain tissues were highly upregulated while breast mammary tissue, ovary tissue and adipose tissue were highly downregulated contrast to other tissues.