Pharmacogenomic Study of SARS-CoV-2 Treatments: Identifying Polymorphisms Associated with Treatment Response in COVID-19 Patients

TL;DR

This study identifies genetic factors that influence how well patients respond to treatments for COVID-19, aiming to improve personalized medicine.

Contribution

The study identifies 16 genes and key biological pathways linked to treatment response in COVID-19 patients.

Findings

Significant polymorphisms in 16 genes were found to affect response to corticoid and immunomodulator therapies.

Enrichment analyses highlighted innate immune and drug metabolism pathways as key in treatment outcomes.

The study provides a framework for precision medicine in managing future coronavirus outbreaks.

Abstract

Background/Objectives: The COVID-19 pandemic resulted in 675 million cases and 6.9 million deaths by 2022. Despite substantial declines in case fatalities following widespread vaccination campaigns, the threat of future coronavirus outbreaks remains a concern. Current treatments for COVID-19 have been repurposed from existing therapies for other infectious and non-infectious diseases. Emerging evidence suggests a role for genetic factors in both susceptibility to SARS-CoV-2 infection and response to treatment. However, comprehensive studies correlating clinical outcomes with genetic variants are lacking. The main aim of our study is the identification of host genetic biomarkers that predict the clinical outcome of COVID-19 pharmacological treatments. Methods: In this study, we present findings from GWAS and candidate gene and pathway enrichment analyses leveraging diverse patient samples from the Spanish Coalition to Unlock Research of Host Genetics on COVID-19 (SCOURGE), representing patients treated with immunomodulators (n = 849), corticoids (n = 2202), and the combined cohort of both treatments (n = 2487) who developed different outcomes. We assessed various phenotypes as indicators of treatment response, including survival at 90 days, admission to the intensive care unit (ICU), radiological affectation, and type of ventilation. Results: We identified significant polymorphisms in 16 genes from the GWAS and candidate gene studies (TLR1, TLR6, TLR10, CYP2C19, ACE2, UGT1A1, IL-1α, ZMAT3, TLR4, MIR924HG, IFNG-AS1, ABCG1, RBFOX1, ABCB11, TLR5, and ANK3) that may modulate the response to corticoid and immunomodulator therapies in COVID-19 patients. Enrichment analyses revealed overrepresentation of genes involved in the innate immune system, drug ADME, viral infection, and the programmed cell death pathways associated with the response phenotypes. Conclusions: Our study provides an initial framework for understanding the genetic determinants of treatment response in COVID-19 patients, offering insights that could inform precision medicine approaches for future epidemics.

Figures

Full text

1. Introduction

Coronaviruses are a subfamily of RNA virus that cause a variety of respiratory diseases, including common cold and severe acute respiratory syndrome (SARS). The Beta (β) class of coronaviruses that comprises SARS-CoV-1, MERS-CoV, and SARS-CoV-2 infect the lower respiratory tract, causing a wide variety of symptoms, including cough, anosmia, fever, and headache, among others [1], that can degenerate into pneumonia and may affect the function of other organs, including the liver, heart, kidney, and brain [2]. SARS-CoV-2, responsible for the COVID-19 pandemic, caused 777 million reported cases and 7 million confirmed deaths worldwide, according to recent WHO data. The number of COVID-19 cases and deaths per year have decreased drastically due to vaccination and immunisation of the population. However, newer COVID-19 variants have significantly increased infections, and future coronavirus epidemics are highly likely [3,4]. The high mutation rates of coronaviruses; environmental factors (e.g., temperature, humidity, radiation, and pollution); and, particularly, host susceptibility to infection (e.g., genetic predisposition) may contribute to infection susceptibility [5,6].

To date, there is no specific treatment for the immune reaction produced by SARS-CoV-2 infections [2]. The rapid mutation of coronaviruses may lead to the inefficacy of treatments in the long term, especially those targeting replicative mechanisms [7]. The current COVID-19 treatment options are based on therapies used for previous coronavirus outbreaks (i.e., SARS-CoV-1 and MERS) and include anti-inflammatory compounds such as glucocorticoids, antivirals, antibiotic, antiparasitic, and/or anti-inflammatory compounds previously used for other infectious and non-infectious diseases [8]. About 10–20% of patients do not respond to treatment and develop an aberrant immune system response that results in a cytokine storm, causing severe symptomatology and even death [9]. Monoclonal antibodies such as the immunomodulator Tocilizumab are used for the treatment of severely affected COVID-19 patients suffering from cytokine storms, with varied success [10,11]. Furthermore, one of the leading causes of mortality is the adverse reactions (ADRs) induced by COVID-19 treatments. Tocilizumab and corticoids may cause immunosuppression and increased risk of infection when administered at high doses [12]. A recent study showed a superior incidence (4.75-fold) of ADRs in COVID-19 patients compared to non-COVID patients, with Tocilizumab associated with the higher rate of ADRs [13]. The reasons behind treatment failure and adverse reactions are unclear. Identifying the risk factors for treatment failure may help to better select treatment and/or doses in future epidemics.

Pharmacogenetic studies have revealed genetic variants that may influence the clinical outcome by altering drug metabolic rates and/or drug targets. Genetic variants in hepatic cytochrome P450 (CYP) metabolic enzymes, transporter proteins, targeted cytokines, and virus entry proteins may play an important role in the host response to SARS-CoV-2 infections. For instance, about 85% of the current medications are metabolised by CYP enzymes, which are known to harbour genetic variants affecting their metabolic rates. CYPs functional variants have been reported to influence the response to treatments such as Lopinavir, Ritonavir, and Hydroxychloroquine in several studies [14,15]. CYPs variants may also play an important role in the bioavailability of most corticoids, currently used as the first line of treatment for SARS-CoV-2 infections. The enzymes CYP3A4 and CYP2D6 are the main metabolic pathways of several corticoid compounds, including Dexamethasone, Prednisone, Hydrocortisone, and Fludrocortisone [16]. However, the influence of known CYP functional variants on the response to corticoid treatments has been determined mainly in asthma but not in COVID-19 patients [17]. Solute carrier organic anion transporter family member 1B1 (SLCO1B1), ATP-binding cassette subfamily C member 1 (ABCB1), and ATP-binding cassette subfamily B members 1 and 2 (ABCC1 and ABCC2) transporter proteins also play an important role in the pharmacokinetics of corticoid treatments. Genetic polymorphisms in SLCO1B1, ABCB1, ABCC1, and ABCC2 have been associated with the response to the antiretroviral Lopinavir and antibiotics like Azithromycin [18,19]. However, the influence of these genetic variants on the response to treatment of COVID-19 infections has not been investigated to date. Proteins linked to host susceptibility to infections may also be related to the treatment response and sustained inflammation [20]. In our own studies, we have shown that variants in genes encoding diverse immunoregulatory interleukins (IL4, IL6, and IL10) are associated with susceptibility to invasive pneumococcal disease [21]. Genetic variants in IL6, IL10, and C-reactive protein (CRP) genes have also been associated with the severity and community-acquired pneumonia [22]. Immunomodulator treatments directed to inhibit the cytokine storm and, specifically, the IL6 pathway in severe patients may be affected by these genetic polymorphisms. Several studies have described associations between genetic polymorphisms within the IL-6 receptor (IL6R) gene and Tocilizumab response [23]. Understanding the host pharmacogenetic profile can provide useful information to help fight the virus infection and reduce mortality. However, further evidence is required before using pharmacogenetic information for the personalisation of COVID-19 treatments.

In summary, there are currently few predictors of the clinical response to COVID-19 pharmacological treatments. The main aim of our study is the identification of host genetic biomarkers that predict the clinical outcomes of COVID-19 pharmacological treatments. This information will help to personalise the treatment of COVID-19 and other coronavirus infections, improve its efficacy, and reduce patient morbidity and mortality.

2. Materials and Methods

2.1. Patients and Sample Selection

Samples from patients diagnosed with COVID-19 were sourced from the Spanish Consortium for COVID-19 Research (SCOURGE; https://github.com/CIBERER/Scourge-COVID19, accessed on 1 April 2022). The samples were obtained between March and December 2020 in 34 medical centres across 25 Spanish cities. A previous publication contains a detailed description of the study cohort [6]. Briefly, COVID-19 infection was diagnosed through PCR testing or clinical assessment. Sample collection and data management were carried out by biobanks associated with the participating centres following informed consent. The whole project was approved by the Galician Ethical Committee (ref. 2020/197) on 10 April 2020. Additional approval was obtained by Ethics and Scientific Committees of the participating centres. All samples and data were processed using standardised procedures, with data management facilitated through REDCap electronic data capture tools hosted at the Centro de Investigación Biomédica en Red (CIBER). For this study, we selected exclusively the patients treated with corticoids and immunomodulators that are known to modulate the cytokine storm, the critical factor that leads to deterioration and mortality in COVID-19 disease.

2.1.1. Clinical Sample

Patients who had received corticoid (n = 2202; 1442 males and 760 females, mean age = 66 years, SD = ±15) or immunomodulator treatments (n = 849; 633 males and 216 females, mean age = 63 years, SD = ±12) were included in the study. The combined sample of patients treated with immunomodulators and/or corticoids group consisted of n = 2487 (1651 males and 836 females, mean age = 66 years, SD = ±14) (see Table 1).

2.1.2. Response Phenotype

Treatment response was assessed using the following data available in all participating centres: survival at 90 days, admission to the intensive care unit (ICU), radiological affectation, and type of ventilation. More detailed or specific response information (i.e., levels of ferritin; interleukins such as IL-6, IL-1β, and IL-10; C-reactive protein (CRP); D-dimer; lactate dehydrogenase (LDH); Troponin; or Prothrombin time (PT)) was not available for most participants.

Survival at 90 days (yes/no) refers to whether the patients survived up to 90 days after hospital discharge. Patients who did not survive at 90 days included patients who passed away in hospital and patients who were discharged but the sequelae of the coronavirus were too severe and passed away soon after. Patients who died more than 90 days after discharge were considered to have died due to causes other than the coronavirus infection. Admission to the intensive care unit (ICU) (yes/no) was commonly used for severe patients with chances of survival during what is known as the first wave of COVID-19 and may have been influenced by the age and severity of the patients. The radiological affectation (yes/no) refers to the changes or abnormalities that can be observed in radiographic images (X-ray or CT scan findings) in the lungs of patients infected with SARS-CoV-2. These changes may include opacities (areas whiter than normal, indicating fluid accumulation or inflammation in lung tissues); infiltrates (presence of fluid in the lung spaces, which may indicate inflammation or damage); and consolidations (accumulation of fluid, inflammatory cells, or scar tissue) that are indicative of the presence and severity of COVID-19 infection in the lung tissue. In our study, we established the radiological affectation as a binary variable to distinguish between patients who received treatment that severely reduced the inflammation and did not develop radiological affectation and patients that did not improve and developed lung affectation. Radiological affectation values have been employed in predicting disease states, particularly amidst the unprecedented circumstances of the pandemic. While the specific dates of the radiological assessments were not documented, it is presumed that patients underwent treatment upon admission to the hospital, followed by subsequent imaging evaluations [24]. The type of ventilation is a categorical variable (1 = no ventilation; 2 = conventional oxygen therapy; 3 = nasal cannulas and IMV) directly related to the state of the disease. If patients fail to respond adequately to medication, disease progression ensues, resulting in an increased requirement for oxygen therapy. As the disease advances, oxygen therapy may need to be escalated to more invasive measures to prevent hypoxia, a severe and potentially fatal condition characterised by viral pneumonia and marked decreases in blood oxygen levels. This progression can lead to acute respiratory distress syndrome (ARDS), organ failure, and, ultimately, death.

Most patients survived more than 90 days after being discharged, with an overall survival rate of 83.8%. Interestingly, the survival rate was slightly higher among patients treated with immunomodulators (85.8%) compared to those treated with corticoids (82.5%). Among patients receiving immunomodulators, approximately half required admission to the ICU, whereas only 23% of patients receiving corticoids were admitted. Additionally, half of the patients in the immunomodulator cohort required invasive ventilation, while 50% of patients in the corticoid and combination therapy cohorts required non-invasive ventilation. Regarding the radiological findings, only 5% of patients in the corticoid and combination therapy cohorts showed no radiological abnormalities, whereas, in the IMM cohort, only 1.37% of patients were unaffected radiologically.

2.2. SNP Genotyping

The samples were genotyped using the Axiom Spain Biobank Array (Thermo Fisher Scientific, Waltham, MA, USA) in accordance with the manufacturer’s instructions at the Santiago de Compostela Node of the National Genotyping Centre (CeGen-ISCIII; http://archivo.xenomica.org, accessed on 1 April 2022). This array interrogates 757,836 polymorphisms, including rare variants from exonic regions specifically chosen from the Spanish population’s genetic profile. Quality control (QC) of the GWAS results was performed by the Santiago de Compostela Node of the National Genotyping Centre, as previously described [6]. Briefly, the QC was carried out using the PLINK1.9 package and the R platform 4.3.1. Variant exclusion criteria: Variants with minor allele frequency (MAF) < 1%, call rate < 98%, and Hardy–Weinberg equilibrium (HWE) [p < 1 × 10^−10^ as recommended [25]] were excluded from the analyses. X chromosome variants were excluded from the GWAS study and analysed separately in the candidate gene analyses. Sample exclusion criteria: Individuals with a call rate < 98% or whose heterozygosity rate deviated >5 standard deviation (SD) from the mean heterozygosity and individuals with an estimated probability < 20% of pertaining to European ancestry were excluded. In related individuals, one individual of each pair of second-degree relatives were excluded (PI_HAT  >  0.25). After QC, 588.117 variants and 1948 patients who had received corticoids and 696 patients with immunomodulator treatments and with available clinical information were considered. The combined sample of patients treated with immunomodulators and/or corticoids consisted of n = 2181.

Variant Imputation

Genetic variants were imputed using TOPMed version r2 reference panel (GRCh38 [26]) on the TOPMed Imputation Server (https://imputation.biodatacatalyst.nhlbi.nih.gov/, accessed on 3 October 2023). The following post-imputation filtering criteria were applied for inclusion: coefficient of determination R-square (Rsq > 0.3), HWE p > 1 × 10^−6^, and MAF > 1%. This dataset encompassed a total of 15,997,581 genetic markers. Bcftools (version 1.18) a software for managing genetic databases, was used for the SNPID annotation.

2.3. Candidate Gene Analysis

We selected polymorphisms within genes implicated in the severity of the disease, metabolism, and targets of compounds used for the treatment of COVID-19 (see Table 2). Specifically, the selected genes were:

• Genes encoding for cytochrome P450 (CYPs) enzymes involved in the metabolism of corticoids: CYP2C8, CYP2C9, CYP2C19, CYP2D6, CYP3A4, and CYP3A5 [9,14,18,23,27,28].

• Transporter genes related to the bioavailability of drugs: ABCB1, ABCB11, ABCC1, ABCC13, ABCC2, ABCC4, and UGT1A1 [9,18,23,27,29,30,31,32].

• Genes directly related to the viral entry into cells of some coronaviruses: ACE, ACE2, and ABO [18,23,29].

• Genes related with the response to immunomodulators and interferon: FCGR3A, IFNAR1, IFNG, IFNGR1, IFNGR2, IFNAR2, IFNLR1, IFNA16, IL-1A, IL1B, IL2, IL4, IL6, IL6R, IL10, IRF7, and TNFA [9,18,23,27,29,33].

• Genes related with the COVID-19 infection-induced cytokine storm: TLR10, TLR8, TLR7, TLR5, TLR4, TLR3, TLR2, TLR1, IL1R1, IL-1A, IL1B, IL2, IL4, IL6, IL6R, IL10 & TNF-α, IFNAR1, IFNAR2, IFITM3, IFNG, IFNGR1, IFNGR2, IFNLR1, IFNA16, and IRF7 [29,33,34,35,36].

2.4. Statistical Analyses

The association of genetic variants with the selected phenotypes was investigated using linear and logistic regression models and the Plink v1.9 package [37]. The models were adjusted for covariates known to be associated with the outcome of the disease (age and sex) [6], as well as the first 10 ancestry-specific principal components (PCs). Significance was determined at p < 5 × 10^−8^ for the GWAS results. Bonferroni corrections for multiple analyses were applied for the candidate gene results, where the significance threshold was 0.05/number of variants included in the analysis for each gene. All the variants with a p-value lower than the threshold established for each gene were considered significant. Analyses were performed separately for each treatment and phenotype.

2.5. Pathway Enrichment Analyses

A gene set enrichment analysis was performed using the WEB-based Gene SeT AnaLysis Toolkit (www.webgestalt.org, accessed on 26 February 2024) to extract Gene Ontology terms (including cellular component, biological process, and molecular function ontologies).

3. Results

The following sections will describe the results by response phenotypes (survival at 90 days, admission in ICU, type of ventilation, and radiological affectation) and type of study (GWAS or candidate gene studies). Table 3 and Table 4 summarise the results from the GWAS analyses and candidate gene studies, respectively. Figure 1, Figure 2 and Figure 3 illustrate the significant results from these analyses.

3.1. GWAS Results

GWAS statistical analyses revealed several genetic loci associated with the type of ventilation in the IMM cohort and the radiological affectation in the COMB and CORT groups at the genome-wide significance level. No statistically significant association was observed when analysing ICU stay and survival at 90 days. Figure 1, Figure 2 and Figure 3 illustrate the Manhattan plots of the most significant results.

3.1.1. Associations with Type of Ventilation

Numerous variants in the gene ANK3 (rs144347645, rs16915354, rs145299149, p rs142740585, rs115701266, rs116165734, rs16915359, rs16915361, rs149847098, and rs144806783) were associated with the type of ventilation in the IMM group (p < 5 × 10^−8^ for all comparisons; see Table 3 and Figure 1).

3.1.2. Associations with Radiological Affectation

Significant associations were found between radiological affectation and variants in the gene RBFOX1 in the CORT (rs551128984, p = 2.01 × 10^−8^) and in the COMB (rs551128984, p = 3 × 10^−9^ and rs72765129, p = 4.75 × 10^−8^) groups. Variants regulating the expression of the gene ZMAT3 were also significantly associated with this phenotype in the CORT group (rs74370746 and rs78451671, p = 2.33 × 10^−8^). A variant in the MIR924HG gene was associated with radiological affectation in the CORT group (rs36036468, p = 4.99 × 10^−8^). Two variants in the ABCG1 gene were associated with radiological affectation in the COMB sample (rs914110892 and rs112302620, p = 1.38 × 10^−8^) (see Table 3 and Figure 2 and Figure 3).

3.1.3. Associations with ICU and Survival at 90 Days

No statistically significant associations were found at the genome-wide level (p < 5 × 10^−8^) for the ICU admission and 90-day survival phenotypes in any of the comparisons performed in the different groups.

3.2. Candidate Genes Results

Single marker analyses of selected variants in the candidate genes revealed several significant associations after correcting for multiple analyses (see Table 4).

3.2.1. Associations with Survival at 90 Days

Several genes involved in the immune response (TLR5, IFNG-AS1, TLR1, TLR6, and TLR10) contained genetic variants associated with survival at 90 days. Associations were found between survival and the TLR5 variants in the IMM cohort (rs55866312 and rs542741410, p < 5.39 × 10^−3^). Associations were also found with the IFNG-AS1 gene variants in the CORT group (rs12306899 and rs12300716, p = 4.05 × 10^−5^) and the combined COMB cohort (rs12300716, rs12306899, rs10878747, rs10878749, rs7301797, rs7306440, rs2870955, rs7133171, rs7137158, and rs11177059, p < 4.74 × 10^−5^). Significant associations were observed with variants located in an overlapping region shared by the genes TLR1 and TLR6 (rs11933455, rs111530790, rs6849400, rs146468588, rs376523214, rs111980996, rs113668069, and rs148035117, p < 3.5 × 10^−4^) in patients treated with corticoids. We identified one significant polymorphism in the TLR10 gene (rs149895872) associated with 90-day survival in the CORT (p = 9.05 × 10^−4^) and in the COMB (p = 1.05 × 10^−3^) groups.

3.2.2. Associations with Admission to the Intensive Care Unit (ICU)

Significant associations were found between admission to the ICU and a variant in the ABCB11 gene in the IMM cohort (rs3770585, p = 2.55 × 10^−4^). Associations were also found with variants in the CYP2C19 (rs12258243 p = 3.19 × 10^−4^) and ACE2 (rs62578917, p = 3.26 × 10^−4^) genes.

3.2.3. Associations with the Type of Ventilation

Statistical analyses revealed a significant association with a variant (rs6742078) in the complex region of UGT1A, a multigenic region that generates nine UGT proteins in the CORT (p = 5.12 × 10^−4^) and COMB (p = 6.73 × 10^−4^) groups. Two associations were found in variants of the TLR4 gene (rs12377632 and rs7868859, p < 6.35 × 10^−4^) in the CORT patients.

3.2.4. Associations with Radiological Affectation

Associations were found between variants in the IL1A gene (rs3783585, rs2071375, and rs697) and radiological affectation in the CORT (p = 7.83 × 10^−5^, p = 5.63 × 10^−4^, and p = 5.63 × 10^−4^, respectively) and COMB (p = 6.406 × 10^−5^, p = 1.128 × 10^−3^, and p = 1.128 × 10^−3^, respectively) groups (Figure 2 and Figure 3).

3.3. Functional Enrichment Analyses Results

A gene set enrichment analysis was performed for each cohort using the top 5000 genes from the GWAS results ranked by the lowest unadjusted p-values using the Gene Ontology databases (Biological Process, Cellular Component, and Molecular Function).

Biological Process (BP): Overrepresentation of genes involved in the regulation of neuron projection development, regulation of transsynaptic signalling, regulation of membrane potential, and dendrite development was detected in the three cohorts. Other pathways that appeared significantly enriched were small GTPase-mediated signal transduction and developmental growth involved in morphogenesis within the IMM cohort, renal system development, and cell–substrate adhesion in the CORT subgroup and muscle system process and sodium ion transport in the COMB group. Genes involved in the regulation of developmental growth appeared enriched in the IMM and CORT cohorts, while the genes involved in cell junction assembly appeared significantly enriched in the CORT and COMB cohorts.

Cellular Component (CC): Genes involved in nervous system and structural components like adherent junctions were enriched in all the cohorts.

Molecular Function (MF): Most of the enriched molecular functions were related to binding domains, such as alcohol, actin, and steroid binding (IMM); calmodulin and phospholipid binding (CORT); phosphoprotein binding (COMB); and scaffold protein and PDZ domain binding, in more than one cohort. There is also gene enrichment in transporter activity, including organic acid transmembrane, metal ion, monoatomic ion, and gated channel activity, present in all the analysed cohorts. Cyclase and nucleoside-triphosphatase regulator activities were found enriched in various cohorts. Enrichment in adhesion and motor activity was observed in the IMM and CORT cohorts, respectively. Finally, enrichment of glutamate receptor activity was observed in all the cohorts.

4. Discussion

Given the severity of SARS-CoV-2 and the unfortunate deaths resulting from inadequate drug responses, identifying predictive factors could tailor treatment in future coronavirus outbreaks [38]. In this study, we investigated genetic predictors of treatment effectiveness in a large group of COVID-19 patients. Several polymorphisms in the genes involved in immune response and previously associated with infection severity were found to be associated with the treatment response, amongst others.

4.1. Genes Related to Response to Immunomodulators

Several genetic variants in the genes involved in immune response or coding for transporter proteins were associated with the response to immunomodulators.

GWAS analyses revealed several ankyrin3 (ANK3) polymorphisms associated with the type of ventilation in the IMM cohort. ANK3 is a gene mostly related to neuronal development [39,40,41]. A previous study identified the potential role of ANK3 in the PPARα/PPARγ signalling pathway and immune infiltration [39] (see Figure 4). Minor alleles of ANK3 genetic variants were associated with a lower probability of needing invasive ventilation.

Candidate gene association studies revealed several statistically significant associations after Bonferroni corrections. An ATP-binding cassette subfamily B member 11 (ABCB11) polymorphism was associated with admission to the ICU. ABCB11 encodes a protein belonging to the ATP-binding cassette (ABC) transporter superfamily, which facilitates the movement of various molecules across cellular membranes [42]. The ABCB11 rs3770585-A allele was associated with a higher probability of being admitted to the ICU of patients treated with immunomodulators. Several genetic polymorphisms in the family of Toll-like receptors (TLRs) were associated with survival at 90 days. It is well known the function of TLRs in inflammation [43,44,45,46], and TLRs have been suggested as possible targets for treatments against COVID-19 [47,48]. TLRs regulate cytokine expression and indirectly trigger the adaptive immune system through the secretion of pro-inflammatory cytokines such as IL-1, IL-6, and tumour necrosis factor-alpha (TNF-α) [45] (see Figure 4). Minor alleles of TLR5 variants were associated with a lower probability of survival in patients treated with immunomodulators and may help to identify patients requiring alternative treatments.

It is important to note that just a limited number of significant associations were observed in response to the immunomodulators, probably due to the moderate sample size of the cohort of patients treated with Tocilizumab or similar molecules. However, the involvement of the genes associated with IMM treatment with the modulation of the immune response and transport suggests that these could be plausible findings. Nevertheless, confirmation of these associations in a larger sample of IMM-treated patients is required.

4.2. Genes Related to the Corticoid Response

GWAS and candidate gene studies showed that most of the genes associated with the response to corticoids are involved in inflammation and the immune response.

An RNA-binding fox-1 homolog 1 (RBFOX1) polymorphism was associated with radiological affectation in patients treated with corticoids in the GWAS study. This gene regulates tissue-specific alternative RNA splicing. A previous review proposed that overexpression of RBFOX1 inhibits inflammation and oxidative stress-related factors repressing (NF-κB) [49] (see Figure 4). Another study analysing the influence of RBFOX1 in SARS-CoV-2 infection suggested that RBFOX1 may act as an upstream regulator for ACE2 [50]. In our study, the “C” allele of the rs551128984 variant was associated with a lower probability of developing radiological affectation, contributing evidence to the involvement of this gene in SARS-CoV-2 infection and treatment. Other genetic variants not related to inflammation or the immune response were also observed in the GWAS study to be associated with the corticoid response. Two zinc finger matrin-type 3 (ZMAT3) polymorphisms were associated with radiological affectation. This gene is implicated in the regulation of alternative splicing processes, influencing the stability and translation function of RNA [51]. Minor alleles of the rs74370746 and rs78451671 variants were associated with a lower probability of developing radiological affectation. A previous GWAS study comparing symptomatic and asymptomatic patients of COVID-19 suggested a potential correlation between genetic variability in ZMAT3 and COVID-19 severity in the Chinese population [52], a finding that would endorse a possible role in relation to treatment response, although different risk polymorphisms were identified in both studies. A MIR924 host gene (MIR924HG) polymorphism was also associated with radiological affectation in corticoid-treated patients. MIR924HG is a IncRNA that regulates the expression of CELF4, a gene involved in alternative mRNA splicing [53]. The “T” allele of the rs36036468 variant was associated with a lower probability of developing radiological affectation. However, the relation between MIR924HG and the response to corticoids and/or SARS-CoV-2 infections is still to be discerned.

Several polymorphisms in the TLRs family of genes were associated with survival at 90 days in the candidate gene analyses. As explained before, TLRs regulate the secretion of pro-inflammatory cytokines [45] (Figure 4). Minor alleles in the TLR1, TLR6, and TLR10 genetic variants were associated with a lower probability of survival at 90 days in patients treated with corticoids. The “C” allele in the rs12377632 variant of the TLR4 gene was associated with a decreased likelihood of requiring more invasive ventilation, whereas the “G” allele of the rs7868859 variant of the same gene was associated with an increased likelihood of needing more invasive ventilation. Those two variants are in linkage disequilibrium, where the major allele on one variant is correlated with the minor allele of the other. Although no previous study has related TLR and response to corticoids, Dexamethasone inhibits important pathways in the host defence against SARS-CoV-2, such as TLR7 and IFIH1/MDA5 [54], exposing TLRs as a possible target for this treatment. Our results suggest that the TLR1, TLR4, TLR6, and TLR10 variants may influence the outcome during corticoid treatment in COVID-19 patients.

A polymorphism in the cytochrome P450 family 2 subfamily C member 19 (CYP2C19) gene was associated with admission to the ICU. CYP2C19 was selected for its role in the metabolization of many xenobiotics, including anticonvulsive drugs, Clopidogrel, Omeprazole, Diazepam, some barbiturates, and certain corticoids [55]. The CYP2C19 rs12258243-A allele was associated with a higher probability of corticoid-treated patients being admitted to the ICU. Several studies have shown the influence of corticoids on the expression of CYP2C19 and other CYPs [56,57]. The rs12258243 variant is in linkage disequilibrium with the rs12248560 variant (r^2^ = 0.84) that predicts CYP2C19 ultrarapid activity. Our results suggest that the CYP2C19 rs12258243 variant is associated with corticoid metabolism alterations that contribute to variability in the treatment response. Another gene related to drug availability may influence the corticoid response: An UDP glucuronosyltransferase 1A1 (UGT1A1) polymorphism was associated with the type of ventilation. The main function of UGT1A1 is degrading bilirubin, a hormone that activates the PPARα receptor and reduces the inflammation by reducing the production of PCR, TNF-α, and IL-6 [58] (see Figure 4). The “T” allele of the rs6742078 variant was associated with a higher probability of needing invasive ventilation. Previous studies have related UGT1A1 with the response to antivirals [59]. One study has related *UGT1A1**6 (rs4148323) with the response to (CDE-11), a treatment for lymphoma that includes Dexamethasone, Irinotecan, and other compounds [60]. Thus, these results contribute further evidence of the relation of UGT1A1 variants with response to corticoid treatment.

A genetic variant in angiotensin-converting enzyme 2 (ACE2) was associated with admission to the ICU. The ACE2 rs62578917-A allele was linked to a higher probability of ICU admission of patients treated with corticoids. However, the effect of corticoids on the expression of ACE2 remains unclear. The ACE2 gene was selected for study due to the role of its encoded protein in facilitating the entry of the SARS-CoV-2 virus in the host cells [61]. Regarding the relationship between ACE2 and the response to corticosteroids, some researchers have observed reduced ACE2 expression in chronic obstructive pulmonary disease [62]. However, another study reported increased ACE2 expression in asthmatic patients using inhaled corticoid therapies [63]. Nevertheless, our findings suggest that the ACE2 rs62578917 variant may influence the response to corticoid treatments in COVID-19 patients.

Several candidate genes involved in the cytokine storm induced by the infection were associated with the response to corticoid treatments. Three IL1A polymorphisms were associated with radiological affectation in patients treated with corticoids. The IL-1 family comprises various pro- and anti-inflammatory proteins, including IL-1α and IL-1β, which exert pro-inflammatory effects by binding to active and inactive receptors. IL-1α-mediated inflammation likely contributes to the pathogenesis of COVID-19, leading to various pathological alterations through the activation of inflammatory cascades, myeloid cell sensing, and inflammasome activation [29] (see Figure 4). Minor alleles of IL1A variants were associated with a lower probability of developing radiological affectation. The effect of the genetic variants on IL-1α expression or functioning is unknown. However, our results suggest that IL1A variants merit investigation as possible predictors of response to corticoid treatments and their effect on inflammation. Two IFNG Antisense RNA 1 (IFNG-AS1) polymorphisms were associated with survival at 90 days. IFNG-AS1 acts as a positive regulator of interferon-gamma (IFNγ) secretion [64]. IFNγ plays a crucial role in the body’s defence against viruses [35] and is one of the cytokines involved in the cytokine storm [29] (see Figure 4). Our results showed that minor alleles of the IFNG-AS1 genetic variants were associated with a higher probability of survival at 90 days. Although no previous study investigated INFG-AS1 polymorphisms in relation to treatment response, the IFNγ levels have been found increased in some patients resistant to corticoids [36]. Our results suggest that IFNG-AS1 variants may also contribute to response to corticoid treatments in COVID-19 patients.

4.3. Genes Related to Response to Corticoids and/or Immunomodulators

Several genes involved in the immune response and corticoid metabolism were found associated with the response in the combined sample of patients treated with corticoids and/or immunomodulators. However, these findings mainly reflect the results obtained in the CORT cohort, suggesting that its larger sample size influenced the outcome of the analyses.

GWAS analyses revealed two RBFOX1 polymorphisms associated with radiological affectation in the COMB cohort. In our study, minor alleles of RBFOX1 variants were associated with a lower probability of developing radiological affectation, thus providing evidence for the role of RBFOX1 in corticoid and the immunomodulator response in COVID-19 patients. Additionally, two ATP-binding cassette subfamily G member 1 (ABCG1) polymorphisms were associated with radiological affectation in the GWAS study. ABCG1 is responsible for transporting a lipidic across cellular membranes of macrophages and is implicated in regulating cellular lipid homeostasis in other cell types [65]. In our study, minor alleles of ABCG1 variants were associated with a lower probability of developing radiological affectation. A study on murine models revealed a connection between the deficiency of ABCG1 in alveolar macrophage and pulmonary granulomatous inflammation [66]. Another study suggested that ABCG1 expression is downregulated by TLR4, contributing to inflammation and lipid accumulation in vascular smooth muscle cells, mitigating the PPARγ/LXRα signalling pathway [44] (see Figure 4).

Within the selected candidate genes, several IFNG-AS1 polymorphisms were associated with survival at 90 days. As explained before, IFNG-AS1 acts as a positive regulator of IFNγ secretion [64], a key component of the immune system [29] (see Figure 4). Minor alleles of the IFNG-AS1 genetic variants were associated with a higher probability of survival at 90 days. Interestingly, the association was of a higher magnitude than that observed in the CORT cohort, suggesting that there is a relationship between IFNG-AS1 variants and the response to COVID-19 treatments and to corticoids in particular. However, the association of these IFNG-AS1 variants with the response to immunomodulators cannot be confirmed due to the limited number of participants treated with this medication. Three IL-1α polymorphisms were also associated with radiological affectation in the COMB cohort. Minor alleles of the IL1A genetic variants were associated with a lower probability of developing radiological affectation. As explained before, IL-1α participates in the immune response, and it is also related to the cytokine storm [29,34]. It is plausible that, in addition to the corticoid response, IL-1α influences the response to the immunomodulators used to reduce the cytokine storm in COVID-19 patients. The results of this study would support this hypothesis, although no significant association was detected in the subgroup of patients treated with immunomodulators, probably due to its limited sample size. Finally, the rs149895872 variant in the TLR10 gene was associated with survival at 90 days in the COMB sample. TLR10 has been reported to be the only TLR that exhibits anti-inflammatory properties [43] (see Figure 4). The TLR10 rs149895872-C allele was associated with an increased mortality risk in treated patients. However, the properties and mechanisms of action of TLR10 are still not clear. Our results suggest that the TLR10 variant rs149895872 may contribute to the response to corticoid and immunomodulator treatments in COVID-19 patients. However, as in the previous case, the possible association of TLR10 variants with the response to immunomodulators needs to be investigated in a larger sample.

4.4. Pathway Enrichment Analysis

The gene enrichment analyses of the 5000 most significantly associated genes in the GWAs studies revealed enrichment mostly in genes involved in neuronal activity in the biological process and cellular component databases. Previous studies have related COVID-19 infection with neurological inflammation and, consequently, dysregulation of neural cell types [67]. We propose that variants in the neurogenesis can create neurons more susceptible to inflammation. Another study suggested that neurological complications are common in COVID-19 patients [68]. These complications can be considered part of the lack of response of treatments.

However, when analysing molecular functions, overrepresentation of the genes involved in molecular signalling was observed. Interestingly, the most significant cause of severe clinical complications and lack of response to treatment in COVID-19 is the cytokine storm, which is essentially a dysregulation of molecular signalling.

4.5. Study Limitations

Our study has several limitations. In addition to the moderate sample size of the IMM cohort, no detailed clinical data were available on the specific symptoms experienced by the patients after hospital discharge, except survival after 90 days. Furthermore, no information on the virus variants present in the recruited patients was available, as the information was not routinely collected during the first and second waves of COVID-19 in Spanish hospitals. It would be interesting to investigate the correlation between the virus characteristics and treatment outcome, which was not possible in our cohorts.

5. Conclusions

Many of the genes we found related to COVID-19 treatment response interact with the NF-κB factor (RBFOX1, TLR10, TLR2, TLR6, TLR1, TLR5, TLR4, ABCG1, ANK3, UGT1A1, IFNG-AS1, and IL1A) (see Figure 4). This factor is key in the regulation of the expression of cytokines that produce the cytokine storm and one of the main targets of corticoids. Furthermore, the Immunomodulators used for the treatment of COVID-19 are antibodies that usually target IL-6, IL-6R, and other interleukins that are also components of the cytokine storm. Our results suggest that genetic variants in the pathway of the pro-inflammatory NF-κB factor and related to cytokine storms may constitute predictors of the response to treatments for severe coronavirus infections.

In summary, a number of genetic variants in proteins involved in immune response and cell transport were found associated with the response to corticoids and IMM treatments. Although several different genetic variants were found associated with the response to corticoids and with the response to immunomodulators that may help to select the most adequate treatment, further studies are required to confirm their specificity. If replicated, these findings may help to personalise treatments in future severe coronavirus infections.