An overview of the study design is illustrated in Figure 1 . Briefly, with large-scale summary statistics from independent meta-analyses of genome-wide association study (GWAS), we first performed linkage disequilibrium (LD) score regression to explore genetic correlations of T2D and glycemic traits with pulmonary function. Then, bidirectional two-sample MR analyses with multiple complementary analyses were conducted to unpick the causality/directionality of these associations. In order to select more genetic variants as IVs for T2D and pulmonary function, we utilized different datasets to analyze the two directions. This work was reported in accordance with the STROBE-MR guidance (13). Publicly available datasets generated by published GWAS or consortiums were used for analyses. Therefore, no additional ethical approval would be required.

Summary statistics for T2D were derived from the hitherto largest GWAS released by the DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) consortium, which combined 32 studies from up to 807,553 participants of European ancestry (74,124 cases and 824,006 controls) (14). The GWAS model had been adjusted for age, sex, ancestry, and study-specific covariates if needed. In addition to T2D, we also examined several glycemic traits for secondary analyses. Summary statistics for fasting glucose (N = 58,074), fasting insulin (N = 51,750), HbA1c (N = 123,665), and fasting proinsulin levels (N = 10,701) were available on Meta-Analyses of Glucose and Insulin-related traits Consortium (MAGIC), which meta-analyzed several GWASs to discover genetic regions associated with glycemic traits (15–17). Only individuals of European ancestry without diabetes were included in the GWAS. Adjustment had been made for age, sex, ancestry, and study-specific covariates, whenever appropriate. Of note, the original GWAS for T2D, fasting glucose, and insulin provided both genetic association data with and without body mass index (BMI) adjustment. To avoid collider bias (18), we only analyzed genetic data that had not been adjusted for BMI.

We selected single-nucleotide polymorphisms (SNPs) associated with T2D or glycemic traits at the genome-wide significance (p < 5 × 10⁻⁸) from the original literature. Significant variants within 500 kb of a lead (most significantly associated) variant were considered as independent SNPs. Palindromic SNPs with mirror allele frequency <0.3 were regarded strand-ambiguous; and therefore, we removed them in the present investigation. To make the IVs satisfy the MR assumptions, SNPs exhibiting significant associations with both T2D/glycemic traits and pulmonary function were excluded from all the analyses after screening each selected SNP and its proxy (r² > 0.8) in PhenoScanner (19) and GWAS catalog.

Corresponding summary statistics for spirometry indices, including FEV1, FVC, and FEV1/FVC ratio, were acquired from a non-overlapping meta-analysis of 22 studies supported by SpiroMeta consortium, involving 79,055 participants of European ancestry (20). The original GWASs were adjusted for age, sex, and height.

In reverse direction analysis, we retrieved the summary statistics for “best measures” of pre-bronchodilator FEV1 (N = 345,665), FVC (N = 345,665), and FEV1/FVC < 0.70 (Ncase = 55,907, Ncontrol = 297,408) from the UK Biobank in British White of European ancestry, which were generated by the MRC-IEU consortium and a previous MR study on pulmonary function (21). The GWAS model was adjusted for sex and study-specific covariates, where appropriate.

For IV selection strategy, we first extracted SNPs associated with each spirometry index that reached genome-wide significance threshold (p < 5 × 10⁻⁸). Then, these extracted SNPs were pruned using the PLINK clumping algorithm (r² = 0.001 and clumping distance = 10,000 kb) to ensure that they are mutually independent. We further removed SNPs that showed pleiotropic associations with T2D and/or glycemic traits and are strand-ambiguous.

Corresponding summary statistics for T2D were acquired from another GWAS meta-analysis (without UK Biobank participants) released by the DIAGRAM consortium, including 69,033 subjects of European ancestry (12,171 cases and 56,862 controls) (22). Since the GWAS of glycemic traits did not include the UK Biobank population, the abovementioned GWAS for four examined glycemic traits were still used to ascertain causal effects of pulmonary function on glycemic traits. Detailed information on genotyping, imputation, and quality control of GWAS included in this study has been described in the original publications.

To study the extent to which heritability of T2D/glycemic traits and pulmonary function are shared, we conducted LD score regression analysis based on LD Hub (23). Such a method involves regressing summary-level data from millions of genetic variants across the genome on a measure of each variant’s ability to tag other variants locally (its LD score). In contrast to MR, LD score regression usually focuses on the amounts of genetic correlations and could distinguish population stratification and polygenicity in GWAS. All summary-level data were Z-transformed and quality controlled to meet analytical requirements.

We harmonized the SNP-exposure and SNP-outcome data to rule out strand mismatches and ensure the alignment of SNP allele. The Wald ratios were calculated for each SNP and were combined using inverse variance-weighted (IVW) method to derive overall causal estimates (24, 25). Given that IVW assumes balanced pleiotropy and its estimates may be biased in the presence of pleiotropic IVs, several sensitivity analyses including weighted median, MR-Egger regression, MR robust adjusted profile score (MR-RAPS), and MR pleiotropy residual sum and outlier (MR-PRESSO) were performed to further strengthen the causal inference (26–29). In addition, multivariable MR (30) was carried out to adjust for indirect effects of BMI and smoking by using UK Biobank data generated by the MRC-IEU consortium. These two traits were assumed as major confounders in this study. Horizontal pleiotropy was tested according to MR-Egger intercept (p < 0.05) (31). Cochran’s Q test was applied to quantify the heterogeneity in IVW estimates (p < 0.05). The causal association between two traits was oriented using Steiger’s test (32). We also carried out “leave-one-out” analysis by removing each SNP in turn to evaluate whether the pooled results were driven by a single variant.

The variance (R ²) in each trait explained by the SNPs was calculated using a previously reported formula (33). The strength of IVs was measured based on F-statistic. We estimated the a priori statistical power using an online calculator (https://sb452.shinyapps.io/power/).

Statistical analyses were performed using RStudio version 3.6.2 with “TwoSampleMR” and “MRPRESSO” packages. We used the Bonferroni correction (corrected p: 0.05/30 groups for bidirectional analyses = 0.002) to penalize multiple testing. A two-sided p-value <0.002 was deemed as statistically significant.

Our analyses showed significant genetic correlations of two spirometry indices, FEV1 and FVC with T2D (FEV1: r_G = −0.19; SE = 0.03; p = 7e−12; FVC: r_G = −0.22; SE = 0.03; p = 2e−14) and fasting insulin levels (FEV1: r_G = −0.35; SE = 0.07; p = 3e−7; FVC: r_G = −0.38; SE = 0.07; p = 3e−8) after adjusting for multiple testing ( Figure 2 ). There was also suggestive evidence supporting genetic correlations of FEV1 and FVC with fasting glucose and HbA1c. Replication analysis using alternative summary statistics yielded a similar pattern of the primary results.

Through a series of rigorous exclusion criteria, we finally included 192 SNPs for T2D, seven for fasting glucose, six for fasting insulin, 43 for HbA1c, and eight for fasting proinsulin. The F-statistics for all traits of interest were much greater than 10, which implied that our MR analyses are unlikely to be susceptible to weak instrument bias. The characteristics of the SNPs used as IVs are shown in Supplementary Tables 1 – 8 . The information on a priori power calculation is presented in Supplementary Table 9 .

We found a significant association between genetic predisposition to higher risk of T2D and reduced FEV1 (multivariable beta = −0.060; 95% CI = −0.100, −0.021; p = 0.002) and FVC (multivariable beta = −0.086; 95% CI = −0.125, −0.047; p < 0.001) and increased FEV1/FVC ratio (multivariable beta = 0.043; 95% CI = 0.008, 0.077; p = 0.015) in both univariable and multivariable models ( Figure 3 ). The results were non-significant when we applied weighted median and MR-Egger methods, which produced causal estimates with wider CIs. However, the direction of point estimates from all established MR models is basically concordant.

Our results also showed some evidence that genetically predicted fasting glucose levels were negatively associated with FEV1/FVC ratio (univariable beta = −0.270; 95% CI = −0.488, −0.053; p = 0.015) ( Supplementary Tables 10 – 13 ). Similar results were observed in other MR models, with p-values less than 0.05 but not exceeding 0.002 after Bonferroni correction. This causal association was not diminished in multivariable MR accounting for BMI and smoking (multivariable beta = −0.273; 95% CI = −0.463, −0.080; p = 0.015).

For some traits such as HbA1c, there is some evidence to report the presence of heterogeneity among the SNPs. However, MR-Egger regression demonstrated that our results are unlikely to be affected by pleiotropic effects. Steiger’s test confirmed that the direction of causal association is correct. In the “leave-one-out” analysis, no single SNP was found to significantly affect the results.

A total of 205 SNPs for FEV1, 264 for FVC, and 106 for FEV1/FVC < 0.70 were identified as valid IVs, which explained 3.38%, 4.76%, and 1.77% of variance, respectively. The F-statistics for FEV1, FVC, and FEV1/FVC < 0.70 are 59, 65, and 59, respectively ( Supplementary Tables 6 – 8 ). Our study had sufficient power to detect very weak associations of pulmonary function with T2D and glycemic traits ( Supplementary Table 9 ).

Genetically predicted higher FEV1 (univariable OR = 0.77; 95% CI = 0.63, 0.94; p = 0.011) and FVC (univariable OR = 0.82; 95% CI = 0.68, 0.99; p = 0.038), but not FEV1/FVC ratio < 0.70 (OR=1.20; 95% CI = 0.72, 2.01; p = 0.480), was associated with lower risk of developing T2D based on multiple MR methods ( Figure 4 ). Multivariable MR yielded a weaker but still suggestively significant association between genetically predicted FEV1 and the T2D risk after conditioning for BMI and smoking (multivariable OR = 0.80; 95% CI = 0.65, 0.98; p = 0.032). However, the causal effect of FVC on T2D was less significant in multivariable analyses (multivariable OR = 0.86; 95% CI = 0.71, 1.05; p = 0.146) ( Figure 4 ).

We found suggestive evidence for association of genetically predicted FEV1 (multivariable beta = 0.126; 95% CI = 0.030, 0.222; p = 0.010) and FVC (beta = −0.095; 95% CI = 0.019, 0.171; p = 0.015) with fasting proinsulin levels in most MR models. In addition, there are suggestive associations of genetically predicted FVC (beta = −0.041; 95% CI = −0.077, −0.006; p = 0.025) and FVC/FEV1 < 0.70 (beta = 0.419; 95% CI = 0.096, 0.742; p = 0.013) with fasting insulin in a single MR model ( Supplementary Tables 14 – 16 ).

MR-Egger intercept showed some evidence of possible pleiotropy in the analysis of FEV1/FVC < 0.70 and fasting insulin (p = 0.016). The direction of the associations was detected to be corrected by Steiger’s test. “Leave-one-out” analysis suggested that no single SNP has substantial impact on the overall results.