Genetic Variants in Genes Related to Lung Function and Interstitial Lung Diseases Are Associated with Worse Outcomes in Severe COVID-19 and Lung Performance in the Post-COVID-19 Condition

TL;DR

This study finds that genetic variants linked to lung function and interstitial lung disease are associated with worse outcomes in severe COVID-19 and lung performance after recovery.

Contribution

The study identifies specific genetic variants associated with severe COVID-19 outcomes and post-COVID-19 lung function.

Findings

TERT rs2736100 and THSD4 rs872471 variants are linked to ARDS severity in COVID-19.

FAM13A rs2609255 is related to post-COVID-19 pulmonary function.

Genetic variants show differential epistasis and correlation in disease outcomes.

Abstract

Genetic variants related to susceptibility to chronic respiratory conditions such as interstitial lung disease (ILD) could share critical pathways in the pathogenesis of COVID-19 and be implicated in COVID-19 outcomes and post-COVID-19. We aimed to identify the participation of genetic variants in lung function and ILD genes in severe COVID-19 outcomes and post-COVID-19 condition. We studied 936 hospitalized patients with COVID-19. The requirement of invasive mechanical ventilation (IMV) and the acute respiratory distress syndrome (ARDS) classification were considered. The mortality was assessed as the in-hospital death. The post-COVID-19 group included 102 patients evaluated for pulmonary function tests four times during the year after discharge. Five variants (FAM13A rs2609255, DSP rs2076295, TOLLIP rs111521887, TERT rs2736100, and THSD4 rs872471) were genotyped using TaqMan assays. A multifactor dimensionality reduction method (MDR) was performed for epistasis estimation. The TERT rs2736100 and THSD4 rs872471 variants were associated with differential risk for ARDS severity (moderate vs. severe, CC + CA, p = 0.044, OR = 0.66, 95% CI = 0.44–0.99; and GG p = 0.034, OR = 2.22, 95% CI = 1.04–4.72, respectively). These variants and FAM13A rs2609255 were also related to pulmonary function post-COVID-19. The MDR analysis showed differential epistasis and correlation of the genetic variants included in this study. The well-known variants in recognized genes related to pulmonary function worsening and interstitial disorders are related to the severity and mortality of COVID-19 and lung performance in the post-COVID-19 condition.

Figures

Full text

1. Introduction

Coronavirus disease 2019 (COVID-19) has created significant challenges worldwide, and our understanding of this disease continues to evolve [1]. Several genome-wide association studies have been performed in multiple ethnic groups, in which genes mainly related to the immune response were associated with the severity of COVID-19; for instance, inflammatory signaling (JAK1), immunometabolism (SLC2A5 and AK5) [2], host antiviral defense mechanisms and mediators of inflammatory organ damage (IFNAR2 and TYK2) [3], interferon signaling (IL10RB and PLSCR1), leucocyte differentiation (BCL11A), and blood-type antigen secretor status (FUT2) [4]. In addition, the first genome-wide significant association for the post-COVID-19 condition was for the FOXP4 locus, previously associated with COVID-19 severity, lung function, and cancers. Remarkably, the study has supported the role of pulmonary dysfunction in COVID-19 severity and in the development of long COVID or post-COVID-19 condition [5].

Single-Nucleotide Variants (SNVs) related to susceptibility to pulmonary function and respiratory diseases [6,7] could share critical pathways involved in the susceptibility and pathogenesis of COVID-19. These biological routes, including those involving the FAM13A, DSP, TOLLIP, TERT [8], and THSD4 genes, have not been entirely studied in COVID-19 and post-COVID-19 condition. These genes have been previously related to lung function and susceptibility to chronic and respiratory lung diseases, such as chronic obstructive pulmonary diseases, lung cancer, interstitial lung disease, idiopathic pulmonary fibrosis, and mycobacterium tuberculosis [8,9,10,11,12].

FAM13A has been linked with an increased risk of COPD and is associated with β-catenin signaling. In SARS-CoV-2 infection, Wnt/β-catenin signaling affects immune cell function, suppresses the proliferation and self-renewal of alveolar macrophages, and promotes the production of pro-inflammatory mediators [13]. Moreover, FAM13A is involved in the TGF-β1-induced fibrotic response in the airway epithelium [14], which is recognized as the dominant chronic immune response in COVID-19 [15]. Meanwhile, desmoplakins (DSP) are critical components of desmosomes that maintain cardiac and respiratory epithelial integrity [16]. Mainly, the rs2076295 locus is directly responsible for differential RNA expression of DSP in primary epithelial cells [11], and the low levels of DSP in the lung cause increased expression of epithelial–mesenchymal transition and extracellular matrix genes promoting pulmonary fibrosis [17], one of the major concerns arising after the recovery from COVID-19 [18].

TOLLIP (toll-interacting protein) is essential to immune cell activation, cell survival, pathogen defense, and other biological processes related to inflammation and innate immune response, which are crucial in virus infections like COVID-19 [19]. Indeed, its deficiency results in ACE2 stabilization and elevated SARS-CoV-2 infection [20]. On the other hand, THSD4 (thrombospondin type I domain-containing 4) has been positively correlated with inflammatory markers in patients with COVID-19 [21].

Finally, the TERT gene responsible for producing human telomerase reverse transcriptase (hTERT) is closely related to the viral RNA-dependent RNA polymerase (RdRP) [22] and linked to inadequate suppression of the innate immune response, increased risk of developing a cytokine storm, lung injury, and severe COVID-19 [23]. The rs2736100 variant has been related to the telomeres’ length [24], and it has been described that short telomeres limit the alveolar epithelial regenerative potential [25] and could be involved in the prognosis of patients recovered from COVID-19.

We hypothesize that SNVs in these lung function candidate genes are implicated in COVID-19 pathogenesis. This study aims to identify the association of FAM13A, DSP, TOLLIP, TERT, and THSD4 variants with severe COVID-19 and post-COVID-19 condition in a Mexican mestizo population.

2. Results

2.1. Clinical and Demographical Characteristics

We included 936 subjects with COVID-19 diagnosis. The patients were classified according to survival (survivors and non-survivors), IMV requirement (IMV and non-IMV), and ARDS severity (mild, moderate, and severe). These groups’ clinical and demographic data are presented in Table 1 and Table S1. The survival comparison presented higher differences than the other classifications. Older subjects with comorbidities (systemic arterial hypertension (SAH), cardiopathies, and preexisting chronic respiratory diseases (CRD, i.e., chronic obstructive pulmonary disease, obstructive sleep apnea syndrome, and/or asthma)) and lower levels of PaO_2_/FiO_2_ at hospital admission were more frequently in the non-survivors group, and most of the patients in this group required IMV during their hospital stay.

Meanwhile, patients requiring IMV were older, had a more extended hospital stay, and presented lower PaO_2_/FiO_2_ at hospital admission and higher mortality than patients in the non-IMV group. The frequencies of males, co-morbidities, and smokers were similar in both groups, although all were slightly higher in the IMV group (Supplementary Table S1).

For patients classified according to ARDS severity, differences were observed in the PaO_2_/FiO_2_ at hospital admission and outcomes related to the severity of the disease, such as days of hospital stay, mortality, and IMV requirement. The age and co-morbidities were not statistically significant between the mild, moderate, and severe groups (Table 1).

2.2. Genetic Association of FAM13A, TERT, DSP, TOLLIP, and THSD4 Variants with the Survival and Severity of Patients with COVID-19

The alleles and genotype frequencies were not statistically significantly different in the survival comparison. However, a marginal p-value was observed for the recessive model in the DSP rs2076295 variant. Then, the G allele of this variant was found to be associated with non-survival risk in the linear regression analysis adjusting for age, SAH, cardiopathies, CRD, PaO_2_/FiO_2_ at admission, and the IMV requirement (p = 0.047), even when age, PaO_2_/FiO_2_, and IMV requirement were only considered as co-variables (p = 0.039). Additionally, the TG genotype was found to have a low risk of death when the last three variables were included as co-variables in the analysis (p = 0.011, OR = 0.58, 95% CI = 0.39–0.88), as well as for the TT + TG genotypes (p = 0.021, OR = 0.63, 95% CI = 0.43–0.93). The remaining variants were not statistically significant when the logistic regression was performed (Table 2).

The frequencies were not statistically different in the IMV and non-IMV comparisons, even in the logistic regression model, which included age and PaO_2_/FiO_2_ at admission as co-variables. However, trending values were observed for the THSD4 rs872471 variants (Supplementary Table S2).

In the analysis with the ARDS severity, the A allele of the TERT rs2736100 variant was observed as a low-risk factor when compared to the mild vs. moderate groups, and C was a risk allele when the moderate group was compared to the severe group. Likewise, in the genetic association models, the conjunction of CA + CC genotypes was observed as a low-risk factor for the moderate severity of ARDS. In contrast, the GG genotype of the THSD4 rs872471 variant was more frequent in the severe ARDS group compared to the moderate ARDS group. Nevertheless, none of these variants remained associated when the logistic regression model was performed, adjusting for PaO_2_/FiO_2_ at admission, steroid administration, and the days since symptom onset (Table 3).

2.3. Genetic Association Study in the Post-COVID-19 Condition

One hundred and two patients were followed after COVID-19 discharge and included in the post-COVID-19 group. The patients were both male and female, had a median age of 56 years, were overweight, and around one-third of the group smoked and had co-morbidities such as T2DM and/or SAH. Most of the patients required IMV during their hospitalization (Table 2). More than half of the group (63.7%) were classified as having moderate ARDS during their hospitalization due to COVID-19, 13.7% as mild, and 22.5% as having severe ARDS.

Table 4 shows the values of the pulmonary function tests at each follow-up. The median value for each test was similar, and no statistical differences were observed between the four assessments.

Genetic variants were also determined in the post-COVID-19 group. The allele and genotype frequencies are shown in Supplementary Table S3. Overall, the frequencies were similar to those reported in the previous tables (Table 3, Tables S2 and S3) for the different COVID-19 groups.

We wondered if some of the studied variants were associated with susceptibility to the post-COVID-19 condition. Therefore, we compared the frequencies of the post-COVID-19 group with those of the survivors without post-COVID-19, but we did not find any significant differences.

We also evaluated whether there were differences in the values of the pulmonary function tests according to the genotypes of the studied variants. We performed a generalized linear model for each variable in the pulmonary function test at the four follow-ups (Table 4) as the dependent variable, and each studied genetic variant served as the independent variable (Table 1), adjusting for age and sex as covariates due to their known influence on pulmonary capacity [26,27].

The influence of genetic variants on pulmonary function test scores was observed at different times. In the 4th follow-up, the DL_CO_ (%) of TT subjects was different from that of the GG carriers of the FAM13A rs2609255, and the FEV1 (%) varied according to the TERT rs2736100 genotypes when the pulmonary function tests were performed for the first time at hospital discharge (Supplementary Figure S1). Moreover, the influence of the DSP rs2076295 genotypes was observed in the FEV_1_/FVC values of the 3rd follow-up and the FEV_1_, FVC, and ∆SpO_2_ (%) after 6MWT after a year of the hospitalization discharge (4th follow-up, Supplementary Figure S1). Finally, the herein-associated THSD4 rs872471 variant was related to the meters reported by the 6MWT in the first two follow-ups and the FEV_1_ and the DL_CO_ values (both in percentages) of the last follow-up (Supplementary Figure S2). Although significant differences were observed in the mentioned cases, we acknowledge the underpower of the analyses due to the small sample sizes observed mainly for the homozygous minor allele in some genetic variants. The remaining analyses of the different pulmonary function tests and genetic variants were not statistically significantly different.

2.4. Multifactor Dimensionality Reduction (MDR) Analyses

We looked for the influence of gene–gene interactions on clinical outcomes. The sample sizes were balanced using the under-sampling method. The best MDR model was selected when the values of testing accuracy were >0.55 and maximum cross-validation consistency was obtained [28]. The TOLLIP rs111521887 variant was not included in the analyses due to its low minor allele frequency.

The best MDR model for the mortality study had a testing accuracy of 0.57, 10/10 of cross-validation consistency, and included the three SNVs DSP rs2076295, FAM13A rs2609255, and THSD4 rs872471. The distribution of cases and controls according to the different combinations is shown in Figure 1a. Some combinations could be observed as high-mortality-risk genotypes. The interaction graph showed a synergistic interaction between the DSP and FAM13A and between the FAM13A and the THSD4 variants. Antagonism effects were observed for the interaction of THSD4 with DSP and TERT. In contrast, no association was observed for TERT with FAM13A and DSP loci (Figure 1b).

Meanwhile, when the IMV requirement was evaluated, the best model had a testing accuracy of 0.56 (10/10 cross-validation) and included in the model the four assessed variants (Figure 2a). The interaction graph mainly revealed synergistic interactions (epistasis) (Figure 2b), with the strongest synergism for the TERT and FAM13A loci. Similarly, a powerful synergistic interaction was found between TERT and FAM13A when the severity of ARDS was compared in the analysis (moderate as controls and severe as cases; Figure 3a), but antagonism effects were observed for DSP with TERT and FAM13A loci. The best model for ARDS severity included the combination of the TERT and FAM13A variants and presented a testing accuracy of 0.66 (10/10 cross-validation) (Figure 3b).

3. Discussion

In our study investigating 936 hospitalized patients with COVID-19, we identified an association between a decreased risk of death and the TG heterozygous genotype and the recessive model (TT + TG genotypes) in DSP rs2076295. Recent research on COVID-19 highlights a relationship of the DSP gene with the lncRNA PIRAT (PU.1-induced regulator of alarmin transcription) that, during SARS-CoV-2 infection, drives the expression of alarmin S100A8 and S100A9, which together promote the differentiation of monocytes to macrophages and contribute to myeloid imbalances during severe COVID-19 [31]. In idiopathic pulmonary fibrosis (IPF), the G allele has been identified as the risk variant and associated with a decreased DSP RNA expression, resulting in a dysfunction in cell migration and the remodeling process in the lungs. In this sense, the protective effect observed with the common allele could be related to maintaining airway epithelial integrity [11,17]. However, the trade-off observed for IPF and COPD [17] can also be observed herein, since in the post-COVID-19 group, the highest pulmonary test scores were found for the patients with the GG genotype, suggesting differential activity of desmoplakins according to the timely disease process (acute or chronic) and the patient’s physiological status (i.e., age, comorbidities, and inflammatory status).

Meanwhile, THSD4 has been studied multiple times in COPD [32]. Although little has been explored for COVID-19, we found that in the genetic association models, the GG genotype of the THSD4 variant rs872471 was more frequent in the severe ARDS group. This relationship is likely seen in the interaction that the THSD4 gene has with the metabolism of the family of extracellular calcium-binding thrombospondins involved in inflammation, apoptosis, coagulation, wound healing, and the interaction with extracellular components, receptors, and cytokines that could aggravate the disease [32]. COVID-19 patients had significantly higher levels of thrombospondins than healthy controls and were positively correlated with inflammatory markers, including ESR, CRP, PCT, ferritin, and biochemical parameters [21]. Moreover, the THSD4 gene has been related to two chemokine genes (CCL18 and CXCL12) and a TNF family receptor that participates in respiratory tract mucosal immunity and influences lung function, which is affected in moderate and severe cases of COVID-19 [33], as we found that it was related to the 6MWT, the % FEV_1_, and DL_CO_. Precisely, thrombospondin-1 has been correlated with % FEV_1_ in patients with COPD and is involved in activating the TGF-β pathway [34].

The A allele of the TERT rs2736100 variant was observed as a low-risk factor compared to the mild vs. moderate group, and C was observed as a risk allele when the moderate group was compared with the severe group. Recent research has pointed out the relationship between the TERT gene and the consistent presence of short telomeres in patients with severe COVID-19 symptoms [35]. As evidence grows that short telomeres increase the risk of severe COVID-19, several molecular mechanisms have been explained to address the connection between short telomeres and severe COVID-19 and suggest that T cell telomere length (TL) affects the adaptive immune response and the innate immune response to SARS-CoV-2 infection. Individuals with long T cell telomeres show a strong T cell response and a strong suppression of the innate immune response accompanied by moderate activity of the innate immune response. In contrast, individuals with short T cell telomeres show weak T cell response and inadequate suppression of the innate immune response accompanied by strong activity of the innate immune response, expressed as a storm of cytokines, lung injury, and severe COVID-19 [35]. Short telomeres were associated with an increased risk of critical illness, defined as admission to the intensive care unit (ICU) or non-ICU death. Finally, lung tissue from patients with very short telomeres shows signs of senescence in structural and immune cells [36].

Finally, SNVs in the FAM13A gene have been associated with an increased risk of COPD and IPF [8]; the most critical part of the FAM13A protein is its N-terminal extension containing the Rho-GAP domain involved in lung endothelial barrier function, which is often deregulated in lung diseases [9]. Patients with COVID-19 exhibit stronger humoral immunity, elevated Rho-GAP and mTOR pathway activities, and higher IFN-I signaling [37]. FAM13A is also associated with β-catenin signaling that is generally activated during injury repair and tissue regeneration [38]. In SARS-CoV-2 infection, Wnt/β-catenin signaling leads to the assembly of an “unconventional” β-catenin–HIF-1α complex that impairs immune cell function, suppresses proliferation and self-renewal of alveolar macrophages, and simultaneously promotes the production of pro-inflammatory mediators in vitro and in vivo [39]. In the post-COVID-19 group, we found that the T allele of rs2609255 variants was associated with lower diffusing capacity of the lungs for carbon monoxide, which agrees with previous studies including patients with IPF [40].

In summary, according to our findings, the involvement of DSP, TERT, and THSD4 during the acute phase of COVID-19 is related to the inflammatory process, cell migration, cytokines release, and the adaptive and innate immune responses. Meanwhile, in the chronic manifestation (post-COVID-19 condition), the main problem is related to the limited regeneration of respiratory cells, as has been previously described [41], characterized by the participation of FAM13A in injury repair and tissue regeneration, DSP in the maintenance of airway epithelial integrity, TERT mediating cell proliferation, and the promotion of fibrotic mechanisms of THSD4 via the TGF-β pathway.

Our study is not exempt from limitations. Firstly, it employed a single-center design and was conducted in a third-level hospital for respiratory diseases, which was focused on severe or critical COVID-19 during the pandemic; thus, patients with mild or asymptomatic disease were not included in this study, nor was a control group, since the diagnostic tests were prioritized for patients and healthcare professionals. Second, the sample collection period was extensive, potentially including different SARS-CoV-2 variants we did not consider. In addition, we did not follow up on mortality during the post-COVID-19 stage; therefore, the results are focused on survivor patients who have been discharged after one year since COVID-19 and in which there could also be follow-up attrition rates. However, our results provide some perspectives on the influence of SNVs in genes related to ILDs, respiratory traits, and severe COVID-19 outcomes.

4. Materials and Methods

4.1. Subjects Included

Patients were recruited in the Instituto Nacional de Enfermedades Respiratorias Ismael Cosio Villegas (INER) at Mexico City, Mexico, a third-level referral center for respiratory disease, where severe cases of COVID-19 were attended. We evaluated two groups of patients with a confirmed diagnosis of COVID-19 through an RT-PCR SARS-CoV-2 test (Supplementary Figure S3); the first included those hospitalized in the INER, consecutively enrolled from August 2020 to December 2022, and who were >18 years old. Patients were excluded from the analyses when their information (outcome, clinical, and demographic variables) was unavailable. This group was further classified according to three clinical outcomes: (1) non-survival, evaluated as in-hospital deceased; (2) requirement of invasive mechanical ventilation (IMV); and (3) severity of the acute respiratory distress syndrome (ARDS) according to PaO_2_/FiO_2_ upon hospital admission (mild > 200, moderate 100–200, and severe < 100). Only <1.0% of the patients included in this study were administered at least one dose of the vaccine; therefore, this information was not considered in the analyses.

The second group comprised patients with a post-COVID-19 condition followed up for 12 months (four follow-ups every three months) after discharge from severe COVID-19 in the INER. All of the patients in this group were diagnosed with severe COVID-19 when hospitalized (April 2020 to August 2021) and required respiratory support through IMV or high-flow nasal cannula oxygen therapy. Fifty-seven subjects from this group were also included in the survivors group. Therefore, they were only included in the post-COVID-19 group for the survivors vs. post-COVID-19 study. They showed pulmonary dysfunction determined by a decrease in their forced vital capacity or desaturation in the 6 min walking test (6MWT) and/or interstitial thickening in the computed tomography, as previously reported [42,43]. The tests performed in each follow-up were spirometry, DLCO, and 6MWT.

4.2. Genotyping

Genomic DNA was isolated from blood cells using standard techniques. The FAM13A rs2609255, TERT rs2736100, DSP rs2076295, TOLLIP rs111521887, and THSD4 rs872471 variants were genotyped using TaqMan ^TM^ assays. These intron SNVs are found in genes selected according to their relevance in pulmonary diseases and lung function. These variants affect regulatory motifs and the expression evaluated from eQTL studies (Table 5) according to HaploReg v4.2 [44].

4.3. Statistical Analysis

Categorical data are presented as absolute numbers (frequency in percentage). Continuous data are described as medians (interquartile range, Q_1_–Q_3_). The normality distribution of the data was assessed with the Shapiro–Wilk and/or Kolmogorov–Smirnov tests, as appropriate. χ2 test and Mann–Whitney U or Kruskal–Wallis tests were used to compare variables between groups. The genetic association study was assessed in PLINK v1.07 [45], and a logistic regression model was performed to adjust for those co-variables that had statistical differences between the comparison groups (see Table 1). The evaluation of the pulmonary function test scores according to genotype and co-variable adjustment was evaluated using a generalized linear model. The statistical tests were conducted in RStudio v4.1.2 [46]. A p-value < 0.05 was set as statistically significant in all cases.

In addition, we evaluated gene–gene interactions using the non-parametric method multifactor dimensionality reduction (MDR) v3.0.2 [29,30]. Three analyses were performed for each evaluated phenotype (mortality, IMV requirement, and moderate vs. severe ARDS). The comparison groups were unbalanced; therefore, we performed under-sampling, which involves randomly removing members of the over-represented group from the dataset until it is balanced. Subjects with missing values were not included. All procedures were performed using ten-fold cross-validation. The selection criteria for a final candidate model were based on the best cross-validation consistency, training accuracy, and testing accuracy (>0.55). A permutation test determined the statistical significance of the training and testing accuracy for the best candidate model.

5. Conclusions

The well-known variants in recognized genes related to pulmonary function worsening (THSD4 rs872471) and interstitial disorders (FAM13A rs2609255, DSP rs2076295, TERT rs2736100) are related to the severity and mortality of COVID-19 and lung performance in the post-COVID-19 condition.