Parenchymal fibrosis, a late feature of severe COVID-19 disease, shares characteristics with idiopathic pulmonary fibrosis (IPF). The impact of rs35705950-T, a functional polymorphism upstream of the MUC5B gene and an established risk factor for IPF, on COVID-19 outcomes in an ancestrally diverse population is unclear.
Rs35705950-T was associated with fewer COVID-19 hospitalizations in a trans-ancestry metaanalysis conducted in the Million Veteran Program (MVP, ntotal=511,965; OR=0.89 [0.82-0.97]) and in joint meta-analysis with the Host Genetics Initiative (ntotal=506,174; OR=0.90 [0.86-0.95]). In subgroup analyses of MVP participants of European ancestry, rs35705950-T was associated with fewer post-COVID-19 pneumonia events, with evidence supportive of a protective dose-response relationship for each copy of the rs35705950-T allele. These results support a potentially protective effect of rs35705950-T against COVID-19 hospitalizations and post-infection pneumonia events.
A respiratory disease caused by a novel coronavirus, now known as SARS-CoV-2, was first reported towards the end of 2019. Despite massive public health measures and vaccination initiatives, the coronavirus disease-2019 (COVID-19) pandemic remains a major global health threat. By February 2022, the COVID-19 pandemic had caused more than 424 million confirmed infections resulting in more than 5.8 million deaths worldwide (1).
Parenchymal fibrosis is a late complication of severe respiratory infections due to COVID-19 (2)(3)(4). Among chronic lung diseases, idiopathic pulmonary fibrosis (IPF) (5), a disorder characterized by progressive pulmonary scarring which is associated with a median survival of 2-3 years in the absence of lung transplantation (6), shares several risk factors with those for severe COVID-19 disease, including advanced age (7), cardiovascular disease, diabetes, and a history of smoking (5). Thus, common pathological processes may be shared between the fibrotic response towards COVID-19 infection and those underlying IPF.
IPF likely develops from a multifaceted interaction between genetic and environmental factors, age-related mechanisms, and epigenetic profibrotic reprogramming (8,9). One of the most robust genetic risk factors identified for IPF susceptibility is rs35705950-T, a common G to T transversion located approximately 3 kb upstream of the mucin 5B, oligomeric mucus/gelforming (MUC5B) gene (10,11). Laboratory evidence supports that rs35705950-T is: 1) a functional variant located within an enhancer subject to epigenetic programming and 2) contributes to pathologic mis-expression in IPF (12). Notably, while rs3705950-T has been robustly associated with increased susceptibility, the same allele has also been associated with decreased mortality in IPF (13), though whether this paradoxical effect is attributable to pleiotropy or index event bias remains controversial (14).
Given the high minor allele frequency (MAF) of rs35705950-T (~11% among individuals of European ancestry) and possible shared pathophysiological pathways between IPF and severe COVID-19 disease, we examined the association between rs35705950-T and the clinical outcomes of COVID-19 infection in the Million Veteran Program (MVP), a multi-ancestry cohort of over 650,000 U.S. Veterans with detailed electronic health record (EHR) and genotyping data (15). Following our primary analysis in the MVP, we validated our results with a comparable analysis conducted in the Host Genetics Initiative (HGI), a global collaboration of over 160 genetic studies assembled to facilitate rapid discovery and dissemination of COVID-19 related science (16).
Data from the MVP, a multi-ancestry genetic biobank sponsored by the United States Veterans Affairs (VA), were analyzed (15). All protocols were approved by the VA Central Institutional Review Board and all participants provided written informed consent. Genotyping was performed using a custom Thermo Fisher Axiom genotyping platform (MVP 1.0) which included direct genotyping of rs35705950-T. Ancestry was defined using harmonized ancestry, race, and ethnicity (HARE) derived from self-report and genetic ancestry data (17). Individuals from three major ancestry groups: European (EUR), African (AFR), Hispanic (HIS) were included.
Demographic and pre-existing comorbidity data were collected from questionnaires and the VA EHR; "pre-COVID" data was from the time of enrollment into the MVP to September 30, 2019.
The cohort demographics and a description of the clinical conditions for all genotyped MVP participants and COVID-19 positive MVP participants (18) that were evaluated in this study are shown in Table 1. (22), rheumatoid arthritis with interstitial lung disease (RA-ILD) (23) and idiopathic pulmonary fibrosis (IPF)diagnosis (24). These variables were derived from the data within two years prior to the index date of diagnosis (Table E1).
Post-index pneumonia was defined from ICD codes within 60 days after the index date.
Association with COVID-19 pneumonia events (pneumonia60d) were performed among patients who received COVID-19 PCR testing at VA sites (Table E2). The ICD codes used to define pneumonia events within 60 days COVID-19 infection (pneumonia60d) are presented in Table E3.
The analysis recommendations from the COVID-19 HGI were utilized to test for associations between the rs35705950-T allele and COVID-19 outcomes. First, we conducted analyses within MVP by each ancestry group using plink2a. Inverse-variance weighted meta-analyses were then performed with summary statistics from Release 5 (01/18/2021) of the HGI using GWAMA (18) (Additional details in supplemental methods).
Firth logistic regression as implemented in the R (v3.6.1) package "brglm2" (version 0.7.1) was used in the sensitivity analyses as it provides a bias-reduced estimate in the setting of small sample sizes and is most powerful for analyses of genetic mutations. This was relevant since adjustment for multiple variables can lead to small sample sizes within covariate categories.
Associations between COVID-19 outcomes and rs35705950-T were performed separately by ancestry, with adjustment for age, age 2 , sex, BMI, CCI, smoking history, asbestosis, RA-ILD and IPF (Table E1), and first 20 principal components. The meta-analyses were performed using random-effects models in "metafor" (version 2.4-0).
Interactions between COVID-19 infection status and rs35705950-T on the outcome of COVID-19 pneumonia within 60 days (pneumonia60d) were assessed using a multiplicative interaction term followed by stratified analyses by COVID-19 infection status. The Firth logistic model for interaction included the independent variables of interest: COVID-19 infection status, rs35705950-T, and their multiplicative interaction term, adjusted for pre-index pneumonia (yes/no, within 2 years pre-index), age, age 2 , sex, BMI, CCI, smoking history, asbestosis, RA-ILD and IPF diagnosis, and first 20 ancestry-specific principal components (PC1-20). The SNP rs35705950-T was modeled as a continuous variable (additive genetic model with values 0,1,2) so that the interaction odds ratio (OR) is equal for heterozygous vs homozygous WT and homozygous mutation vs heterozygous. The additive interaction of COVID-19 and rs35705950-T (for every increase in one allele, i.e., 0 to 1, or 1 to 2 copies) was also assessed with the estimate of excess risk due to interaction (RERI) and the 95% confidence interval estimated with the Delta Method.
Associations between rs35705950-T and pre-existing comorbid conditions and laboratory values were examined using clinical data prior to the COVID-19 pandemic (Sept 2019).
Individuals with ≥2 Phecodes(25) were defined as cases. Phecodes with <200 cases within each ancestry group were excluded, resulting in 1618 (EUR, 1289 (AFR), 994 (HIS) Phecodes. LabWAS was conducted for 69 clinical tests; for individuals with repeated measures, the median of the individuals' EHR record was used. Logistic/Firth regression and linear regression were used for Phecodes and laboratory measurements, respectively. A Bonferroni-adjusted p-value threshold of 1.2 x 10 -05 (0.05/3901) accounted for all the models tested across three ancestries for significance. Analyses were performed using PLINK2(26) (Additional details in Supplemental Methods).
Our study included 19,168 COVID-19 positive patients from three major ancestry groups (EUR, AFR, HIS). The minor allele frequency for rs35705950-T was 9% among the entire MVP and 8% among COVID-19 positive individuals (Table 1); among these 14.4% (n=2,758) and 0.7% (n=138) were carriers of 1-copy and 2-copies of MUC5B rs35705950-T, respectively (Table 1, Figure E1).
ancestry meta-analysis within MVP (MVP-ALL), and joint meta-analysis with HGI (META) are shown in Figure 1 and Table 2. In the joint meta-analysis, the most significant association was between rs35705950-T and fewer hospitalization events compared with population controls (OR = 0.90 [0.86-0.93], p=8.99 x 10 -5 , Figure 1 and Table 2). The MUC5B rs35705950-T allele was not associated with reduced COVID-19 positivity in trans-ancestry meta-analysis within MVP (OR=0.98 [0.95-1.01], p=0.06) but was nominally significant (p<0.05) in the joint metaanalysis with HGI (OR=0.97 [0.95-1]; p=0.02). The MUC5B rs35705950-T allele was not associated with severe COVID-19 disease (critically ill).
To further evaluate the robustness of our findings from the primary analyses, we performed sensitivity analyses among MVP participants from 3 ancestry groups (EUR, AFR, HIS) who had COVID-19 testing performed within the VA (n=136,164). Associations between rs35705950-T and COVID-19 positivity and outcomes with adjustment for additional covariates, including BMI, smoking status, CCI, asbestosis, RA-ILD and IPF diagnosis, which were not available in HGI, were performed in MVP participants. We observed similar effect sizes as the primary analyses (Table E4), with the association between rs35705950-T allele and reduced risk of hospitalization, compared with non-hospitalization, due to COVID-19 remaining robust (OR= 0.86 [0.77, 0.95], P=0.004, Table E4). The OR and 95% CI were very similar to the primary analyses (OR=0.88 [0.81, 0.96]. Notably, rs35705950-T was associated with a reduced risk of testing positive for COVID-19 (OR=0.95 [0.90, 0.99], P=0.03) in this MVP sub-population. The MUC5B rs35705950-T allele was not associated with severe outcomes plus mortality (OR = 0,89 [0.67-1.20], p= 0.45) or mortality alone (OR = 0.93 [ 0.74-1.17], p=0.56; Table E4).
Among MVP participants of European ancestry who tested positive for COVID-19 (n=8541), the adjusted odds ratio for post-index pneumonia was 18% less with each additional MUC5B rs35705950-T allele (OR = 0.82 [0.72, 0.93], p=0.001). Among COVID-19 negative MVP participants, the adjusted odds for post-index pneumonia non-significantly increased with each additional MUC5B rs35705950-T allele (OR=1.06 [0.98, 1.15], p=0.13). Increasing copy numbers of the MUC5B rs35705950-T allele were associated with progressive reduction of the post-index pneumonia risk, with adjusted ORs of 10.9 [10.2, 11.7], 8.41 [7.36, 9.61], and 6.47 [4.93, 8.48] in patients with 0, 1, and 2 copies of MUC5B rs35705950-T respectively among COVID-19 positive relative to negative subjects (p < 0.0001, Table 3). This differential effect of an additional MUC5B rs35705950-T allele on post-index pneumonia in COVID-19 positive vs. COVID-19 negative patients was statistically significant (multiplicative scale interaction OR = 0.77 [0.66, 0.89], p=0.0004) in EUR (Table 3). There was also evidence of a negative additive interaction, with the relative excess risk due to interaction (RERI) corresponding to an increase of one copy of rs35705950-T allele estimated as -2.07 (95% CI -3.24, -0.9). This suggests that the relative risk for pneumonia in COVID-19 positive patients is -2.07 less for each increase in one copy of MUC5B rs35705950-T allele than if there were no interaction between COVID-19 and MUC5B rs35705950-T. The number of post-COVID-19 pneumonia events in HIS or AFR was too low to permit further analysis.
To explore clinical conditions and biomarkers associated with the MUC5B rs35705950-T allele which may impact the susceptibility and severity of COVID-19, PheWAS and LabWAS using pre-COVID-era data (through Sept 2019) were performed. The sample sizes for MVP participants included in PheWAS and COVID-19 association studies are shown in Table 1 (Figure E1). The results of the PheWAS are shown in Figure 2 with full summary statistics presented in Table E5.
In PheWAS analyses between rs35705950-T and PheCodes with > 200 cases, significant associations were identified exclusively with pulmonary processes. Increased risk for "other alveolar and parietoalveolar pneumonopathy" (Phecode 504, OR = 2.64 [2.50 -2.78], P = 7.07 x 10 -289 ) and "post-inflammatory pulmonary fibrosis" (Phecode 502, OR = 2.85 [2.65 -3.05], P = 8.90 x 10 -186 ) were significantly associated (P bonferroni < 1.2 x10 -5 ) with rs35705950-T in all three ancestry groups (EUR, AFR, HIS). Among the EUR and AFR groups, rs35705950-T was also associated with increased risk for idiopathic fibrosing alveolitis (Phecode 504.1). Notably, there were no significant associations identified between rs35705950-T and influenza infection (Phecode 481) or bacterial pneumonia (Phecode 480.1) in any of the ancestry groups. The power to detect a difference with these conditions was >95% as there were 4,728 cases of influenza and 10,579 cases of bacterial pneumonia in EUR cohort. There were also no significant associations between rs35705950-T and non-pulmonary conditions potentially relevant to COVID-19 outcomes, including pulmonary embolism/deep vein, systemic inflammatory syndrome / sepsis, or acute renal failure (Figure E2 and Table E5).
The LabWAS of MUC5B rs35705950-T with median values of pre-COVID clinical laboratory tests is shown in Table E6 and Figure E3. Among EUR participants, 10 laboratory tests were significant after Bonferroni adjustment, with the majority belonging to white blood cell counts or fractions, with an increase in monocyte count also significant in the HIS subgroup.
There were no significant associations among the AFR subgroup.
Our study supports that the "T" allele of rs35705950 in MUC5B, which has been associated with an increased risk for the development of IPF, confers a decreased risk for COVID-related hospitalization among MVP participants; among participants of European descent, a decreased incidence of pneumonia following COVID-19 infection was also observed. The protective effect of the rs35705950-T may appear to contradict previous work demonstrating an increased risk for acute respiratory distress syndrome (ARDS) in the pre-COVID era (27) as well as the increased risk of severe COVID-19 disease observed for other well-established causal variants of IPF located in the TERC, DEPTOR, and FAM13A(28) genes. However, our findings are consistent with previous studies conducted in European cohorts (29) and in the HGI (28) which support that pathophysiological changes due to rs35705950-T may result in distinct interactions with SARS-CoV-2 which confer the differential risks observed. The rs35705950 polymorphism is located within an enhancer region of the MUC5B gene, the protein product of which is a major gel-forming mucin in the lung that plays a key role in mucociliary clearance and host defense (30)(31)(32)(33). Consistent with this, mouse knockout models for Muc5b demonstrate increased susceptibility towards bacterial infections of the respiratory tract and persistent inflammation. In contrast to the loss-of-function of knockout models, the "T" allele of rs35705950 results in a gain-of-function and is associated with enhanced expression of the MUC5B transcript in lung tissue from human subjects (24)without clinical IPF. Notably, while excess MUC5B protein is observed in the epithelial cells of respiratory bronchioles and honeycomb cysts (31,32,34) of individuals with IPF, differential expression in lung tissue by rs35705950 genotype has not been consistently observed in patients with advanced disease and may be attributable to the universally elevated MUC5B protein levels in IPF patients regardless of their rs35 705950 allelic configuration (24,32,34).
The functional impact of increased MUC5B expression on clinical outcomes in populations both with and without IPF remains incompletely understood. In a study which antedated the COVID-19 pandemic, a modestly increased risk for the development of ARDS, a major feature of severe SARS-CoV-2 infection, was observed among homozygotes for rs35705950-T who were >50 years old, none of whom had clinical or radiographic evidence of IPF on post-hoc chart review (27). It should be noted, however, that the underlying etiologies for ARDS in the population examined were heterogeneous, with over half of cases attributed to either trauma or non-pneumonia-related sepsis. In contrast to the increased risk for pre-COVIDera ARDS, the rs35705950-T allele was associated with decreased risk for acute respiratory exacerbation events among non-Hispanic white ever-smokers with interstitial lung abnormalities (ILA), defined as non-dependent parenchymal infiltrates which are frequently considered subclinical precursors to clinically-apparent interstitial lung disease (ILD), in the COPDGene study (35). Whether increased mucin production confers protection against viral infections (36), which are believed to contribute to a substantial proportion of acute respiratory exacerbation events in chronic lung disease (37), should be explored in future studies.
Our analyses suggest that the decreasing risk of post-index pneumonia associated with rs35705950-T may be specific to COVID-19. First, to our knowledge, there have been no previously published reports of differential susceptibility or clinical outcomes in non-COVID-19 respiratory viral or bacterial infections by rs35705950-T genotype. Within our own MVP data, there were no associations between pre-pandemic influenza infection or bacterial pneumonia and rs35705950-T allele. Second, although rs35705950-T carrier status or number of allele copies did not impact testing rates for COVID-19 (which were 22.2%, 22.3%, and 22.4% for individuals of European descent with 0, 1, and 2 copies of rs35705950-T, respectively), individuals with the "T" allele demonstrated a trend towards decreasing rates of testing positive in MVP, supporting that rs35705950-T may modulate an individual's susceptibility to infection by SARS-CoV-2. Third, in addition to attenuating the risk for post-index pneumonia exclusively among COVID-19 positive individuals, differential risk by rs35705950-T copy number was observed, supporting a dose-response relationship for this outcome.
We did not observe associations between rs35705950 and severe COVID-19 illness or mortality. This may be due to the low incidence and/or multifactorial causes of these severe outcomes which may include non-respiratory disorders, such as shock and multi-system organ failure. This is supported by the PheWAS analysis, which found nearly exclusive associations between rs35705950-T with pulmonary processes, with no evidence of association (with >80% power to detect differential effects) with other processes potentially relevant to severe COVID-19 disease, including thromboembolic disease, septic shock, and acute renal failure.
The COVID-19 pandemic amplified healthcare disparities due to socioeconomic factors and likely contributed to differences in outcomes (40). However, heterogeneity in the magnitude of disparities between healthcare systems exists, with numerous studies demonstrating reduced levels of disparities and relatively equitable access to care among VA healthcare users (relative to non-VA healthcare users)(41-43); thus, differences in outcomes attributable to socioeconomic factors within the VA may be less than those in the general population. Consistent with this, a study by Trivedi et al (44), reported minimal changes in 30day mortality rates for heart failure and pneumonia hospitalizations following adjustment for a comprehensive panel of socioeconomic factors, including poverty, housing, education, and rurality. Within our cohort, access to care was uniformly high, with >98% of individuals having 1 primary care visit within the preceding 18 months and comparable numbers of patients ≥ with 1 clinical encounter or admission in the preceding year. However, due to the ≥ complexities of socioeconomic factors and their potential impact on and interactions with rs35705950-T, we acknowledge that future studies in this area are needed.
The strengths of our study include a large, ethnically-and geographically-diverse cohort with harmonized prospective outcomes data, directly genotyped rs35705950 data, and the use of robust statistical approaches. In addition to interrogating the relationships between rs35705950-T and clinical outcomes, the availability of rich clinical phenotyping data permitted exploration of the underlying pathobiological mechanisms through PheWAS and LabWAS and support that the differential risks observed by rs35705950-T are likely mediated by pulmonaryspecific processes. Despite these strengths, we acknowledge the following limitations. First, although the MVP is one of the largest and most diverse genomic medicine databases established to date, participants are predominantly male and of European ancestry; this, and the lower minor allele frequency of rs35705950-T among non-European populations may have impacted our power to detect associations between rs35705950-T genotype and clinical outcomes among females and racial/ethnic minorities. Second, due to the timing of the data freeze for our analysis, only short-term outcomes were captured. F uture studies with extended follow-up to explore the potential relationship between rs35705950-T and post-COVID syndromes, including lung fibrosis, are needed. Third, although our analyses were adjusted for critical variables, residual confounding and bias may still exist. Given the association between rs35705950-T and an increased incidence of subclinical ILA and respiratory symptoms (45), individuals harboring the MUC5B variant may exhibit behavior modifications (increased testing or more stringent self-isolation ) which could have introduced biases not be captured or adjusted for in our analyses (28). Fourth, viral subtypes, vaccination status and treatment approaches evolved with time which may have impacted severity outcomes.
Fifth, we acknowledge that the effect estimates associated with rs35705950-T in our study are modest and may limit its application in clinical or predictive algorithms. However, despite these limitations, we assert our results, which identify a potentially protective role of rs35705950-T in
Copy number=2 6.47 (4.93, 8.48) < 0.0001
Multiplicative Scale OR (95% CI) 0.77 (0.66, 0.89) p=0.0004
Additive Scale RERI (95% CI) -2.07 (-3.24, -0.9)
Odds ratios and 95% confidence intervals are estimated from Firth logistic regression adjusting for pre-index pneumonia, age, age 2 , sex, BMI, CCI, smoking, IPF, asbestosis, RA-ILD, and principal components PC1-20, including a multiplicative interaction between additive MUC5B rs35705950-T allele and COVID-19 infection. An estimate of the relative excess risk due to interaction (RERI) and 95% confidence interval (Delta method) corresponding to an additive interaction effect is estimated from the logistic regression model. of European ancestry with no missing data for pre-or post-index pneumonia events after sensitivity analysis.
Table E3 SDR group including the VA National Surveillance Team provided researchers with curated data extracted from electronic health records (EHR) of Veterans, after stripping off identifiable elements, and ensuring data security by requiring restricted access to the data(4). This activity was undertaken centrally by an experienced VA research team in order to ensure consistent applications of dates across the EHR tables, as well as uniform definitions of events preceding and subsequent to COVID-19, in an effort to characterize the trajectory of COVID-19 and other diseases and conditions. VA-wide efforts describing these data are already in public domain (5). Structured data obtained via CPT and ICD9 or ICD10 codes, deposited in the VA Corporate Data Warehouse (CDW),was further enriched with rulebased unstructured events recorded in patient notes via natural language processing (NLP), as was done previously for other projects (6). The purpose of NLP-boosting was to fill gaps in knowledge about the severity of the disease, as well as extract specific dates when procedures (e.g. intubation and extubation) were performed. These curated data were then provisioned in the Million Veteran Program (MVP) study mart for COVID-19 in VINCI ensuring data security behind the VA firewall.
The primary analysis MUC5B SNP rs35705950-T with COVID-19 infection and hospitalization included all the patients who were tested from Feb 2020 through April 2021 using RT-PCR based method (7,8)
For the PheWAS analysis on the MUC5B rs35705950-T allele, we used Phecodes ( 14) and laboratory measures from the clinical data available prior to the onset of COVID-19 infection (Sept 2019).
Individuals with two or more Phecodes were defined as cases, while those without were defined as controls. In our previous study, we found that cases under 200 provide >80% statistical power and low type I error therefore we excluded phecodes with fewer than 200 cases within each ancestry group. It yielded 1618 (EUR), 1289 (AFR), 994 (HIS) Phecodes. We also performed a laboratory-wide association study using 69 clinical labs extracted from the EHR of MVP participants. We used the median of each person's entire lab history. We used logistic regression for Phecodes and linear regression for laboratory measurements in PLINK2 (15). When the logistic regression model failed to converge for binary outcomes, firth regression was used. Regression models were adjusted for sex, age (at the time of data freeze), quadratic term of age, and the first 20 principal components. We applied Bonferroni correction to each ancestry specific analysis to adjust for multiple hypothesis testing, and following threshold were use select significant associations: EUR = 3.09 x 10 -05 (0.05/1618), AFR = 3.8 x 10 -05 (0.05/1289), EUR = 5.03 x 10 -05 (0.05/994).
AJRCCM Articles in Press. Published June 30, 2022 as 10.1164/rccm.202109-2166OCCopyright © 2022 by the American Thoracic Society