Increased expression of hypoxia-induced factor 1α mRNA and its related genes in myeloid blood cells from critically ill COVID-19 patients

Abstract Background COVID-19 counts 46 million people infected and killed more than 1.2 million. Hypoxaemia is one of the main clinical manifestations, especially in severe cases. HIF1α is a master transcription factor involved in the cellular response to oxygen levels. The immunopathogenesis of this severe form of COVID-19 is poorly understood. Methods We performed scRNAseq from leukocytes from five critically ill COVID-19 patients and characterized the expression of hypoxia-inducible factor1α and its transcriptionally regulated genes. Also performed metanalysis from the publicly available RNAseq data from COVID-19 bronchoalveolar cells. Results Critically-ill COVID-19 patients show a shift towards an immature myeloid profile in peripheral blood cells, including band neutrophils, immature monocytes, metamyelocytes, monocyte-macrophages, monocytoid precursors, and promyelocytes-myelocytes, together with mature monocytes and segmented neutrophils. May be the result of a physiological response known as emergency myelopoiesis. These cellular subsets and bronchoalveolar cells express HIF1α and their transcriptional targets related to inflammation (CXCL8, CXCR1, CXCR2, and CXCR4); virus sensing, (TLR2 and TLR4); and metabolism (SLC2A3, PFKFB3, PGK1, GAPDH and SOD2). Conclusions The up-regulation and participation of HIF1α in events such as inflammation, immunometabolism, and TLR make it a potential molecular marker for COVID-19 severity and, interestingly, could represent a potential target for molecular therapy. Key messages Critically ill COVID-19 patients show emergency myelopoiesis. HIF1α and its transcriptionally regulated genes are expressed in immature myeloid cells which could serve as molecular targets. HIF1α and its transcriptionally regulated genes is also expressed in lung cells from critically ill COVID-19 patients which may partially explain the hypoxia related events.


Introduction
The coronavirus disease 2019 (COVID-19) epidemic caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has rapidly developed into a devastating pandemic [1]. As of today, the World Health Organization (WHO) has reported more than 46 million persons diagnosed with COVID-19 and over 1.2 million deaths worldwide [2]. It is associated with significant mortality in high risk patients, with poor prognostic features upon admission, such as advanced age, as well as a history of obesity, diabetes and hypertension. The spectrum of the disease is broad, and includes pneumonia, sepsis, and acute respiratory distress syndrome (ARDS) [3]. This complicated course occurs in the context of a cytokine storm characterized by overproduction of TNF, IL6, IL1ß, and G-CSF and generalized vascular hyperpermeability [3]. Viral particles spread through the respiratory mucosa infecting other cells, unleashing a series of immune responses, characterized by a reduction of T and B lymphocytes and an increment of monocytes and neutrophils [4][5][6][7]. The pathogenesis of this severe form of COVID-19 is poorly understood.
Hypoxaemia, defined as a decrease in the partial pressure of oxygen is an ominous sign of COVID-19, and it is usually an indicator of disease severity [8,9]. An oxygen saturation above 90% is associated with better outcomes [10]. Over 80% of COVID-19 patients in the intensive care unit have severe hypoxaemia [11]. A kind of "silent hypoxia" in which COVID-19 patients deteriorate rapidly without warning and develop respiratory failure has been described [12]. Hypoxia indicates an imbalance of oxygen delivery to tissues and leads to compromised function, which is quantitatively related to organ, tissue and even cell type. The hypoxia-inducible factors (HIF) are considered master regulators of oxygen homeostasis and are oxygen level sensitive [13,14]. HIF1a is a heterodimeric transcription factor that bind to hypoxia response elements, which participates through the regulation of the expression of several genes in numerous cellular events such as O 2 sensing, glucose metabolism, lipid metabolism, angiogenesis and other aspects of endothelial biology [14]. Currently, there is scarce information regarding the expression of HIF in patients with severe COVID-19 and its potential involvement in the immunopathogenesis of this condition. Therefore, in the present work we carried out scRNAseq to identify the cell populations present in critically ill COVID-19 patients and to determine the expression of hypoxiainduced factor 1a (HIF1a) and its related genes.

Patients and tissue samples
Blood samples from five critically ill patients with COVID-19 were collected in EDTA-coated tubes. Tissues were collected from patients diagnosed, treated and followed at the Internal Medicine Department of the Hospital de Especialidades, Centro M edico Nacional Siglo XXI of the Instituto Mexicano del Seguro Social in April 2020. A family member of each participating patient signed an informed consent and the study protocol was approved by the Comisi on Nacional de Etica e Investigaci on Cient ıfica del Instituto Mexicano del Seguro Social in accordance to the Helsinki declaration. SARS-CoV-2 infection was corroborated by RT-qPCR at an official federal government reference laboratory.

Sample preparation, scRNAseq library generation and sequencing
Peripheral blood from the five critical COVID-19 patients was collected in EDTA-coated tubes, and immune cells were isolated according to standard centrifugation methods followed by red blood cell lysis.
Chromium Next GEM Single Cell 3 0 Reagent Kits v3.1 and the protocol from 10X Genomics was followed as recommended by the manufacturers. Briefly, immune cells were pooled in a single tube and were suspended in 1x phosphate buffered saline (PBS) to 700-1200 cells per ml. Cell suspension was loaded in Chromium Next GEM Chip G and sorted in the Chromium Controller from 10X Genomics. The Cell-Gel Beads in Emulsion (GEMs) were then incubated to generate the barcoded cDNA. cDNA was cleaned using Dynabeads and washed, followed by cDNA amplification and SPRI selection. The retrieved cDNA was enzymatically fragmented, end-repaired, poly-A tailed and ligated. Size selection, adaptor ligation and amplification were done. Sequencing was done using NextSeq 550 System High-Output Kit (300 cycles) in NextSeq 500 system (Illumina) according to 10X Genomics specifications: Read 1 ¼ 28 cycles, Read 2 ¼ 91 cycles, Index 1 ¼ 8 cycles. All quality control steps were carried out using 4200 TapeStation System (Agilent) with High Sensitivity D1000 Screen Tape, whereas the concentration was calculated using Qubit 2.0 Fluorometer with Kit High Sensitivity assays.

scRNAseq bioinformatic analysis
Partek Flow software was used with scRNAseq toolbox. First the tags were trimmed and then the reads were aligned using STAR 2.7.3a algorithm to human genome hg38. UMI's were deduplicated and barcode filtered. Following criteria were then applied to each cell, i.e. gene number between 200 and 6000, UMI count above 300 and mitochondrial gene percentage below 20%. To quantify the transcriptome human hg38 Ensembl transcripts release 99 was used. Counts per million, add 1.0 Log 2.0 were the normalization parameters. Healthy donors' datasets were downloaded from 10X Genomics website and analyzed using Loupe Browser from 10X Genomics. Data has been deposited in Sequence Read Archive hosted by National Centre for Biotechnology Information under accession number PRJNA635580.

Dimensionality reduction and clustering
The filtered and normalized gene-barcode matrix was analyzed by principal components, then graph based and uniform manifold approximation and projection (UMAP) for dimension reduction using default parameters was carried out.

Cellular alterations analysis
ENRICHR (https://maayanlab.cloud/Enrichr/) was used to understand the biological meaning behind the resulting list of genes, to analyze the gene ontology and pathway information of significantly deregulated genes in the different cellular subsets.

In silico analysis of RNAseq files of cells derived from bronchoalveolar lavage fluid
Paired end fastq files were downloaded from BIGD-GSA and SRI from NCBI. Accession numbers are CRA002390 and SRR10571724, SRR10571730, and SRR10571732. Data files consisted in 3 COVID-19 patients bronchoalveolar lavage fluid (BALF) and 3 bronchoalveolar fluid lavage controls. Fastq files quality control of was performed using FastQC v. 0.11.9. Samples adapter were trimmed with the Cutadapt v tool. 2.10.
BALF samples duplicates were spotted using Picard v. 2.22.4, the RNAseq reads alignments were done using STAR tools v. 2.7.3a, with the parameters previously described [20]. The RNAseq reading count was carried out using the Rsubread package with the annotations of Homo sapiens GRCh.38.99 downloaded from Ensembl. Subsequent analyzes were performed using the stable Gene ID identifier with BioMart annotations: Ensembl Genes 100, Human genes (GRCh38.p13), the identifier was translated to HGNC symbol to facilitate biological interpretation. The BALF counts were subjected to quality control using the NOISeq package, biases associated with length, GC content, RNA composition and batch effect, were identified and removed. BALF samples were processed using full-quantile normalization for GC content and length, finally, the ARSyN function was used to eliminate batch effect. The normalization process was carried out using the NOISeq and EDASeq packages in R. Differentially expressed genes were calculated on the normalized counts using the noiseqbio function of the NOISeq package, using the predefined parameters and filtered by a probability greater than 0.99.

Clinical features of critically ill COVID-19 patients
Five male patients with a mean age of 47.8 ± 6.6 years (range 41-57) were studied. Interestingly, four of five patients presented hypoxaemia upon hospital admission to the hospital with O 2 saturation levels ranging from 40-84% of. The initial SOFA score of the five patients ranged between 2 and 4. All patients had radiological evidence of alveolar occupation and ground glass appearance on imaging studies and required intubation and invasive mechanical ventilation 3 to 7 days after admission as well as hemodynamic collapse that required vasopressor support. All five patients were treated with enoxaparin, hydroxychloroquine, azithromycin or clarithromycin as well as lopinavir and ritonavir; three received high dose glucocorticoids (hydrocortisone or dexamethasone). Three of the five patients died of ARDS and multiorgan failure (Table 1).

Immature myeloid cell populations in critically ill COVID-19 patients
We have previously reported that, as compared to healthy adults, lymphoid cell subsets, such as B and T lymphocytes as well as NK cells, are present in low quantities in critically ill COVID-19 patients, whereas cells of myeloid origin predominate. Interestingly, immature myeloid cell populations, such as band neutrophils, metamyelocytes, promyelocytes-myelocytes, monocytoid precursor, and immature monocytes prevailed. Mature lineages such as segmented neutrophils, mature monocytes and monocyte-macrophages were also observed ( Figure 1) [20].

Altered cellular pathways in immune cells from critically ill COVID-19 patients
Myeloid cell population subsets presented alterations in: (1) neutrophil-mediated immunity, (2) cellular response to type I interferon and (3) cell cycle control, DNA and RNA processing. The lymphoid cellular subset alterations were mostly related to ribosome biogenesis as well as RNA processing and modification. Along with the gene ontology analysis the algorithm compares the gene list input data with publicly available information and thus, our gene set lists were in accordance with COVID-19 data (Supplementary Figure 1).

HIF1a expression in leukocytes from critically ill COVID-19 patients
Once the blood cell populations were identified, we looked for HIF1a gene expression. As shown in Figure  2(A), HIF1a gene was expressed in all myeloid lineages to a greater extent than in lymphoid cells. This was particularly evident in the mature monocyte population. Since HIF1a is a transcription factor that acts as a trans regulator, we searched for HIF1a-regulated genes potentially involved in COVID-19 immunity. Among these immune related genes, we found an increased expression of CXCL8 or Interleukin-8, a chemokine involved in the migration of mature neutrophils to the site of infection in most myeloid cell subsets, and almost no expression in lymphoid populations and monocytoid precursors (Figure 2(B)). In keeping with the increased expression of CXCL8, the genes for chemokine receptors CXCR2 (Figure 2(C)), CXCR4 ( Figure  2(D)), as well as CXCR1, were also expressed at increased levels in most myeloid lineages. It is noteworthy that lymphoid cells did express the CXCR4 gene, but showed no expression of the CXCR2 gene, which can explain the exacerbated inflammatory response characteristic of these patients [21]. Interestingly, we found expression of Toll like receptor-2 and À4 (TLR2 and TLR4) in most myeloid populations, which could be related to SARS-CoV-2 sensing (Figure 2(E,F)).
We also found expression of superoxide dismutase 2 (SOD2), involved in the metabolism of reactive oxygen species; vimentin (VIM), a type III intermediate filament and plasminogen activator urokinase receptor (PLAUR), which is related to plasminogen activation.
Considering that HIF1a function is controlled by different factors, we also evaluated the expression of the Von Hippel Lindau (VHL) and prolyl-4-hydroxylase (P4HA1) genes, two of the main inhibitors of HIF1a function. Indeed, these two molecules are involved in the ubiquitination and degradation of HIF1a. In keeping with the increased expression of HIF1a and target genes observed, we found that neither VHL nor P4HA1 were expressed by the peripheral blood cells of critically ill COVID-19 patients (Figure 2(G,H)). We next evaluated the potential interaction between HIF1a and their transcriptional targets by assessing the simultaneous presence of their mRNAs in the same single cells ( Figure 4). As shown in Figure  4(A), a significant proportion of myeloid cells coexpressed HIF1a and CXCL8. We also observed a predominance of cells co-expressing HIF1a and TLR2 (Figure 4(B)) as well as HIF1a and SOD2 (Figure 4(C)). It is noteworthy that among the different myeloid subsets, mature monocytes were the ones that exhibited co-expression of most of the analyzed genes.
Finally, we evaluated the expression of HIF1a and its transcriptional targets in peripheral blood cells from healthy donors. To do so, we analyzed 10X Genomics publicly available datasets. We observed that HIF1a expression was lower as compared to the expression observed in COVID-19 patients ( Figure 5). Similar results were observed for genes such as CXCL8, CXCR2, PLAUR, TLR4 and SOD2 ( Figure 5).

HIF1a expression and its transcriptionally regulated genes in SARS-CoV-2 infected lungs
Lung is one of the main organs affected by SARS-CoV-2 infection. We have previously described the alteration of hyaluronan, glycosaminoglycan and mucopolysaccharides metabolism which could, in part, explain the presence of a viscous fluid found in the interstitial space that reduces lung compliance and thus, effective ventilation [23]. Lung cells from COVID-19 patients show expression of HIF1a, TLR2, TLR4, PFKFB3, CXCR1, À2 and À4, SOD2, PLAUR ( Figure 6) and other genes expressed in immature myeloid cells.

Discussion
Using a scRNAseq approach we have found a distinctive predominance of immature cell lineages of myeloid origin and a reduced number of lymphoid lineages in patients with life-threatening respiratory distress syndrome due to SARS-CoV-2 infection, which is in accordance with recently reported data [6][7][8].
Collectively, our results are in accordance with similar observations were severe COVID-19 patients presented marked changes in immune cell composition and phenotype, reflecting an emergency granulopoiesis [24]. It is noteworthy that a similar blood profile has  been observed in patients with acute parvovirus B19, cytomegalovirus and Epstein-Barr viral infection [20]. In line with this, it has been reported in severe COVID-19 patients, as well as, some bacterial infections that a high neutrophil-to-lymphocyte ratio is associated to multi-organ damage and mortality [25].
Hypoxia and inflammation are unequivocally linked [26] and are two of the main physiological consequences of SARS-CoV-2 infection, particularly in severe cases. In this study we present scRNAseq data regarding HIF1a-related gene expression in peripheral blood leukocytes from critically ill COVID-19 and characterized the different cell subpopulations. HIF1a is a heterodimeric transcription factor sensitive to oxygen and induced under hypoxic conditions [27]. The HIF1a trans element can regulate the expression of CXCL8 [28], CXCR1, CXCR2 [29] and CXCR4 [30]. CXCL8 expression can be stimulated by interleukin  (IL) 6, TNFa, hypoxia [31] and viral infection [32] in cells such as monocytes, neutrophils, epithelial cells and fibroblasts [33]. CXCL8 is a chemokine that exerts its pro-inflammatory functions throughout the CXCR1 and CXCR2 receptors. CXCL8 and its receptors contribute to pathogen elimination, through the transient activation of ERK, AKT, SRC and FAK leading to activation of neutrophils [32]. The expression of CXCL8, which is present in COVID-19 patients, is considered a potential prognostic factor in acute respiratory distress syndrome (ARDS) [34] and lung injury [35].
The SARS-CoV-2 viral entry depends upon binding of viral spike (S) protein to the host cell surface protein angiotensin-converting enzyme 2 (ACE2) [36].
Under physiological conditions, there is an equilibrium between ACE1 and ACE2. However, under hypoxic conditions the ACE2 gene is negatively regulated by HIF1a. It has been suggested that increased levels of ACE2 positively correlate with COVID-19, thus the stabilization of HIF1a which in turn, down regulates the ACE2 gene, could improve the outcomes in COVID-19 patients [37,38].
The immunopathological outcomes are most likely induced by the host-virus interaction. The interaction between viral antigen and host immune cells results in an exacerbated inflammatory response [39]. In the present study, we also found high expression of both TLR2 and TLR4 genes in peripheral blood leukocytes of severe COVID-19 cases. The viral Spike protein can be recognized by TLR2 [40] and TLR4 [39]. Upregulation of TLR2 and TLR4 has been found in the context of infections with the original SARS-CoV [41]. TLR4 constitutes one of the most efficient innate immune receptors, triggering pro-inflammatory responses after binding to the pathogenic ligand, and this interaction could be useful for developing drugs [39].
Metabolic reprograming of innate immune cells occurs during hyperinflammatory states. Immune cells contribute to systemic changes in metabolism by altering their metabolic profiles in response to the immunological state. Therefore, therapeutic modulation of immune cell metabolism could alter the inflammatory state and thus improve patient prognosis [42]. Inflammation and hypoxia are inherently linked, and hypoxia is a well-known glycolysis driver as oxygen deficit results in limited OXPHOS [43]. Previous studies have shown that the molecular mechanisms underlying the switch from OXPHOS to glycolysis during innate immune cell response require HIF1a [42]. Among the HIF1a responsive genes, we found expression of those related to carbohydrate metabolism, such as SLC2A3/ GLUT3, PFKFB3, PGK1 and GAPDH [44]. HIF signalling pathway activation in neutrophils results in an increased survival of these cells, b2 integrin expression, production of antimicrobial peptides and glycolysis. Neutrophils use high rates of Warburg-like glycolysis for ATP generation. The absence of HIF1a leads to reduced ATP pools resulting in a profound impairment of the inflammatory response [45]. HIF1a can also regulate nitric oxide production, pentose phosphate pathway, OXPHOS and arginase metabolism [42]. Overall, immunometabolism is now considered an indispensable regulator of immunity, with HIF1a playing a central role, modulating the function of various immune cell subsets [45]. The expression of HIF1a has been previously found to be a sepsis marker [42].
HIF1a participates in the regulation of a plethora of cellular events such as metabolism of ROS trough the regulation of SOD2 [46], the regulation of cytoskeleton trough VIM type III filament, which also participates in inflammation [47], and PLAUR which activates plasminogen and activates a cascade of extracellular proteases [48]. Interestingly, the expression of this gene could be used as a predictor of severe respiratory failure [49] which is consistent with our results.

Conclusion
In conclusion, in the present study, we have demonstrated the expression of HIF1a and its transcriptionally regulated genes, in myeloid cells, including both mature and immature subsets, present in peripheral blood of critically ill COVID-19 patients. The up-regulation and participation of HIF1a in relevant events such as inflammation, immunometabolism, and TLR supports its use as molecular marker for COVID-19 severity and as a potential candidate for targeted therapy.

Disclosure statement
At the time that this project was carried out RCG and CML worked for Analitek S.A. de C.V. which supplied research reagents, but also they were members of Dr. Mercado lab.
Analitek S.A. de C.V. did not influence the results nor design of the protocol. The rest of the authors declare not competing interests.

Ethical approval
All participating patients were recruited with signed informed consent and ethical approval from the Comisi on Nacional de Etica e Investigaci on Cient ıfica del Instituto Mexicano del Seguro Social in accordance with the Helsinki declaration.

Funding
This work was partially supported by grants 289499 from Fondos Sectoriales Consejo Nacional de Ciencia y Tecnologia Mexico and R-2015-785-015 from Instituto Mexicano del Seguro Social (MM). DMR is a recipient of the National Council for Science and Technology Fellowship "Catedra CONACyT" program.

Data availability statement
Data have been deposited in Sequence Read Archive hosted by National Centre for Biotechnology Information under accession number PRJNA635580.

Author contributions
MM, DMR and KTP conceived, designed and coordinated the project, performed scRNAseq experiments, analysed, discussed data and prepared the manuscript. EV, HM, CRGB, JT, AM, NWR, FBF, AFH, CRR, EPM, GSR, SVP, CML, RCG, LBA, RS, PPS and ACG performed scRNAseq experiments, analysed, discussed and interpreted biological data, wrote the paper. CLM, JCG, EFO, GFP, provided the biological sample, retrieved the immune cells and collect clinical data.