Godoy‑Tena et al. Genome Medicine (2022) 14:134 https://doi.org/10.1186/s13073-022-01137-4 RESEARCH Open Access Epigenetic and transcriptomic reprogramming in monocytes of severe COVID‑19 patients reflects alterations in myeloid differentiation and the influence of inflammatory cytokines Gerard Godoy‑Tena1, Anis Barmada2,3, Octavio Morante‑Palacios1, Carlos de la Calle‑Fabregat1, Ricardo Martins‑Ferreira1, Anna G. Ferreté‑Bonastre1, Laura Ciudad1, Adolfo Ruiz‑Sanmartín4, Mónica Martínez‑Gallo5, Ricard Ferrer4, Juan Carlos Ruiz‑Rodriguez4, Javier Rodríguez‑Ubreva1, Roser Vento‑Tormo2 and Esteban Ballestar1,6*    Abstract Background: COVID-19 manifests with a wide spectrum of clinical phenotypes, ranging from asymptomatic and mild to severe and critical. Severe and critical COVID-19 patients are characterized by marked changes in the myeloid compartment, especially monocytes. However, little is known about the epigenetic alterations that occur in these cells during hyperinflammatory responses in severe COVID-19 patients. Methods: In this study, we obtained the DNA methylome and transcriptome of peripheral blood monocytes from severe COVID-19 patients. DNA samples extracted from CD14 + CD15- monocytes of 48 severe COVID-19 patients and 11 healthy controls were hybridized on MethylationEPIC BeadChip arrays. In parallel, single-cell transcriptomics of 10 severe COVID-19 patients were generated. CellPhoneDB was used to infer changes in the crosstalk between mono‑ cytes and other immune cell types. Results: We observed DNA methylation changes in CpG sites associated with interferon-related genes and genes associated with antigen presentation, concordant with gene expression changes. These changes significantly over‑ lapped with those occurring in bacterial sepsis, although specific DNA methylation alterations in genes specific to viral infection were also identified. We also found these alterations to comprise some of the DNA methylation changes occurring during myeloid differentiation and under the influence of inflammatory cytokines. A progression of DNA methylation alterations in relation to the Sequential Organ Failure Assessment (SOFA) score was found to be related to interferon-related genes and T-helper 1 cell cytokine production. CellPhoneDB analysis of the single-cell transcrip‑ tomes of other immune cell types suggested the existence of altered crosstalk between monocytes and other cell types like NK cells and regulatory T cells. *Correspondence: eballestar@carrerasresearch.org 1 Epigenetics and Immune Disease Group, Josep Carreras Research Institute (IJC), 08916 Badalona, Barcelona, Spain Full list of author information is available at the end of the article © The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​ mmons.​org/​publi​cdoma​in/​zero/1.​0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 2 of 22 Conclusion: Our findings show the occurrence of an epigenetic and transcriptional reprogramming of peripheral blood monocytes, which could be associated with the release of aberrant immature monocytes, increased systemic levels of pro-inflammatory cytokines, and changes in immune cell crosstalk in these patients. Keywords: COVID-19, Monocytes, Epigenomics, DNA methylation, Single-cell transcriptomics, Immune cell crosstalk Background complex (MHC-II) antigens [23] and increased activation Severe acute respiratory syndrome coronavirus 2 (SARS- of apoptotic pathways [24]. CoV-2) causes the well-known Coronavirus disease Differentiation and activation of monocytes and other 2019 (COVID-19), which has become a major global myeloid cells are directly associated with epigenetic health burden. SARS-CoV-2 infection occurs through mechanisms [25]. The functional plasticity of these cells the nasopharyngeal mucosa [1]. Subsequent immune is also reflected at the epigenetic level, and several stud- responses occur at the local mucosa and at a systemic ies have shown that DNA methylation profiles, among level. An effective response to SARS-CoV-2 infection other epigenetic marks, vary in response to inflamma- requires coordination between the innate and adap- tory cytokines, hormones, and other factors [26, 27], tive immune systems, including granulocytes, mono- depending on their functionality. Cytosine methylation cytes, macrophages, and T and B cells [2, 3]. The range (5mC) occurs at CpG dinucleotides and is generally asso- of immune responses to SARS-CoV-2 infection is diverse, ciated with transcriptional repression [28], although its from asymptomatic or mild upper-respiratory illness to relationship with transcription depends on the genomic severe viral pneumonia, acute respiratory distress syn- location of the affected CpG sites [29]. In some cases, drome, and death [4]. The most severe forms of COVID- DNA methylation changes occur as a result of upstream 19 are caused by dysregulation of immune homeostasis, environmental effects that link cell membrane recep- which leads to hyperinflammation in the lungs [5]. This tors, signaling pathways, and transcription factors (TFs) has been shown to be more pronounced in the elderly that can either directly recruit DNA methyltransferases and in individuals with pre-existing comorbidities [6, 7]. (DNMT) and ten–eleven translocation (TET) enzymes, Nevertheless, despite the numerous studies performed or indirectly influence their binding to specific genomic in the field, the impact of exacerbated immune responses sites. associated with severe COVID-19 at the systemic level The characterization of the epigenetic and transcrip- remains unclear. tomic reprogramming in monocytes, given their central Various studies have demonstrated that peripheral role in inflammatory responses, is essential if we are to pathogenic T cells and inflammatory monocytes can understand the specific dysregulated pathways involved induce a cytokine storm in severe COVID-19 patients in severe forms of COVID-19. In this study, we obtained [8]. This takes the form of excessive production of inflam- the DNA methylation profiles of peripheral blood mono- matory mediators, specifically, interleukin (IL)-6, IL-1β, cytes of severe COVID-19 patients and studied their granulocyte–macrophage colony-stimulating factor relationship with transcriptomic changes, obtained by (GM-CSF), tumor necrosis factor-alpha (TNFα), and generating droplet-based single-cell RNA sequencing interferon gamma (IFNγ) [8–11]. IFN is essential for (scRNA-seq) data from peripheral blood. inducing the innate immune response during viral infec- tion through different interferon regulatory factors (IRFs) Methods [12]. Further, in COVID-19 patients, type I IFN defi- Human samples ciency appears to be a hallmark of severe cases [13–19] in Our study included a selection of 58 severe COVID- association with persistent blood viral load and an exac- 19 patients from the intensive care unit (ICU) of Vall erbated inflammatory response [14]. d’Hebron University Hospital (Barcelona) recruited dur- Single-cell omics studies have identified specific tran- ing the second wave of infection in Spain (October to scriptional features in monocytes, natural killer (NK) November 2020). Peripheral blood samples were taken cells, dendritic cells (DCs), and T cells associated with at different times following admission of the patient to the severity of COVID-19 [13, 20–22]. These studies the ICU, as specified in Additional file 1. Table S1 (Days have revealed that severe COVID-19 is marked by a dys- in ICU). Ninety-four percent of the patients required regulated myeloid cell compartment [13]. It has also been intubation and all enrolled cases were confirmed to be shown that monocytes from severe COVID-19 patients infected with SARS-CoV-2 using real-time RT-PCR are characterized by a tolerogenic phenotype with at the time of collection. For all enrolled patients, the reduced expression of class II major histocompatibility date of enrollment, clinical classification, or treatment Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 3 of 22 Table 1 Summary of patient cohort for DNA methylation analysis Healthy controls COVID19 severe patients p value Number 11 48 - Age (mean ± SD) 50 ± 11.16 60 ± 11.96 0.0042 Sex (% female) 36.4 25 0.710 SOFA 0 5 ± 2.97 2.4e − 07 IL-6 (pg/ml) (mean ± SD) NA 316.94 ± 1238.82 - Days in ICU (mean ± SD) NA 6 ± 5.93 - Treated with dexamethasone (%) NA 52.08 - Obesity (%) NA 27.03 - Hypertension (%) NA 56.25 - Death (%) NA 33.33 - Mechanical ventilation (%) NA 93.75 - was obtained from the clinical records. From all these bovine serum (FBS) until monocyte purification. The patients, 48 of the 58 patients were selected for DNA monocyte population was isolated by flow cytometry methylation analysis (Additional file 1. Table S1) and (FacsAria Fusion, BD, Beckton Dickinson, San Jose, CA, peripheral blood mononuclear cells (PBMCs) from 10 of USA). PBMCs were stained with CD14-APC-Vio770 the 58 patients were used for droplet-based scRNA-seq (Miltenyi Biotec) and CD15-FITC (Miltenyi Biotec) in analysis (Additional file 2. Table S2). The control popu- staining buffer (MACS) for 20 min. A gating strategy lation for the DNA methylation analysis comprised 11 was employed to eliminate cell debris, doublets, and healthy donors (HDs) recruited at the Blood Bank of Vall DAPI + cells. CD14 and CD15 antibodies were used to d’Hebron University Hospital. Table 1 summarizes the isolate CD14 + CD15 − . Purified cells were pelleted and characteristics and clinical data from patients included stored at − 80 °C. in the DNA methylation analysis. We included an addi- After monocyte isolation, DNA was isolated using the tional group of 14 patients from the same hospital for AllPrep DNA/RNA/miRNA Universal Kit (Qiagen) fol- DNA methylation and expression validation, includ- lowing the manufacturer’s instructions. ing 9 severe COVID-19 patients and 5 mild COVID-19 patients, together with an additional group of 6 HDs. The DNA methylation profiling DNA Methylation™ Kit (Zymo Research, CA, USA) validation cohort was collected during February 2022, Bisulfite (BS) conversion was performed using EZ-96 applying the same selection criteria as for the discov- ery cohort. For the validation cohort, we only included according to the manufacturer’s instructions. Five hun- non-vaccinated patients, to match the vaccination status dred nanograms of BS-converted DNA was hybridized with that of the patients collected in the initial phase of on Infinium Methylation EPIC BeadChip arrays (Illu- the study. Clinical information corresponding to the new mina, Inc., San Diego, CA, USA). These were used to cohort is also included in Additional file 1. Table S1 (vali- analyze DNA methylation. They enable > 850,000 methyl- dation cohort). This study was approved by the Clinical ation sites per sample to be assessed at single-nucleotide Research Ethics Committees of Hospital Universitari resolution, which corresponds to 99% of the reference Germans Trias i Pujol (PI-20–129) and Vall d’Hebron sequence (RefSeq) genes. University Hospital (PR(AG)282/2020), both of which Each methylation data point was obtained from a com- adhered to the principles set out in the WMA Declara- bination of the Cy3 and Cy5 fluorescent intensities from tion of Helsinki. Informed consent was obtained from all the methylated and unmethylated alleles. Background patients before their inclusion. intensity computed from a set of negative controls was subtracted from each data point. For representation and Monocyte purification and DNA isolation further analysis, we used beta (b) and M values. Beta is PBMCs were obtained from peripheral blood by Ficoll the ratio of methylated probe intensity to overall inten- gradient using Lymphocyte Isolation Solution (Rafer, sity (the sum of the methylated and unmethylated probe Zaragoza, Spain) from 48 of the severe COVID-19 intensities). M is calculated as the ­log2 ratio of the inten- patients and 11 HDs. Once PBMCs were isolated, all sities of the methylated and unmethylated probes. For samples were stored at − 150 °C in 10% DMSO in fetal statistical purposes, the use of M is more appropriate Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 4 of 22 since b-values are severely heteroscedastic for highly coefficient (rho) to measure the association of two varia- methylated and unmethylated CpG sites. Raw DNA bles and thereby identify CpG sites in which DNA meth- methylation data are available at GEO, with accession ylation was correlated with SOFA in patients with severe number GSE188573 [30]. COVID-19. We selected the CpG sites for which Spear- man’s rho was greater than 0.4 and had an associated Quality control, data normalization, and statistical analysis value of p < 0.01. Principal component analysis (PCA) of of DMPs b-values from ShinyÉPICo was used to determine the Quality control and analysis of EPIC arrays were per- correlations of PCs with clinical variables such as dexa- formed using ShinyÉPICo [31], a graphical pipeline that methasone treatment, obesity, and hypertension. Pearson uses minfi (v1.36) [32] for normalization, and limma correlation coefficients between numerical variables and (v3.46) [33] for analyzing differentially methylated posi- PCs were calculated. Categorical variables were entered tions. ShinyÉPICo is available as an R package at the in a linear model together with the PCs, which were con- Bioconductor (http://​bioco​nduct​or.​org/​packa​ges/​shiny​ sidered as a function of the variable. epico/) and GitHub (https://​github.​com/​omora​nte/​shiny​ epico) sites. We used the BS conversion control probe Gene ontology, transcription factor enrichment, test included in ShinyÉPICo to determine whether the and chromatin state discovery and characterization conversion rate was above the quality threshold of 2 The GREAT (v3.0.0) online tool (http://​great.​stanf​ord.​ established by Illumina. The threshold was calculated edu/​public/​html) was used for gene ontology (GO) anal- from the information of the BS conversion control probes ysis, in which genomic regions were annotated by apply- of the EPIC arrays. When the BS conversion reaction is ing adapted basal and extension settings (5 kb upstream, successful, control probes display strong signal in the red 5 kb downstream, 1000 kb plus distal). GRCh37 (UCSC channel, whereas if the sample has unconverted DNA, hg19, Feb. 2009) was used as the alignment genome ref- control probes have a strong signal in the green channel. erence. Annotated CpGs in the EPIC array were used The red/green ratio for each control position was calcu- as background. GO terms were considered significant lated for each sample. for a > twofold change and an FDR < 0.05. Enrichment is CpH and SNP loci were removed by the Noob method, represented as − ­log2(FDR). GO categories with p < 0.05 followed by quantile normalization. Sex chromosomes were considered significantly enriched. GO analysis of (X and Y) were also excluded from the analysis to avoid differentially expressed genes (DEGs) was carried out discordant information among samples. Even when data using the online Enricher gene ontology analysis tool were generated in a single batch and randomized, we (https://​maaya​nlab.​cloud/​Enric​hr/). GO categories with applied the batch effect correction. Sex and age of the a > twofold change and an FDR < 0.05 were considered donors were included as covariates, to minimize con- significantly enriched. founding effects due to differences between the median We used the findMotifsGenome.pt tool from the motif age of the patient and control cohorts, and the Trend and discovery HOMER software (v4.10.3) to analyze motif Robust options were implemented in the eBayes moder- enrichment [35]. A flanking window of ± 250 bp from ated t-test analysis. To compare healthy donors with the each CpG was applied, and CpGs annotated in the EPIC entire severe COVID-19 patient cohort, we identified array were used as background. To determine the loca- differentially methylated CpG sites by using t-tests and tion relative to a CpG island (CGI), we used “hg19_cpgs” a method with defined empirical array weights, included annotation in the annotatr (v1.8) R package. The statis- in the limma package [33], and selecting CpGs with a tical test used for the enrichment in these analyses was false discovery rate (FDR) of < 0.05 and a Δβ of > 0.15. To Fisher’s exact test. Chromatin functional state enrich- test the effects of potential changes in monocyte sub- ment of DMPs was measured using public PBMC data set proportions, we also included this information as a from the Roadmap Epigenomics Project (http://​www.​ covariate, and performed the same analysis as above, but roadm​apepi​genom​ics.​org/) generated with ChromHMM including only those samples for which such information (v1.23) [36]. Enrichments were calculated with Fisher’s was available. exact test using array annotation as background regions. We used the iEVORA package (v1.9.1) [34] to iden- Only significantly enriched states are shown. tify differentially variable positions (DVPs). This algo- rithm identifies differences in variance using Bartlett’s Heatmaps and PCA plots test (FDR < 0.001), followed by the comparison of means Heatmaps of DMPs were generated with functions avail- using t-test (p < 0.05) to regularize the variability test, able in the ComplexHeatmap (v2.11.1) and gplots (v3.1.3) which is overly sensitive to single outliers. For the anal- R packages. We used PCA for the low-dimensional analy- ysis in Fig. 2, we calculated Spearman’s correlation ses. PCA projection matrices were calculated with R’s Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 5 of 22 prcomp function, and visual representations of PCs were of each initial cluster produced by the Leiden algo- plotted with the ggfortify package (v4.1.4). rithm. Data were not denoised because no significant contamination or ambient RNA was present. Previously Whole‑genome bisulfite sequencing (WGBS) analysis described scRNA-seq datasets of HDs [41] were then DNA methylation values of Ensembl Regulatory Build integrated for comparison using single-cell variational regions of progenitor cells such as hematopoietic stem inference (scVI) [42] with a generative model of 64 latent cell (HSC), multipotent progenitor (MPP), common variables and 500 iterations. More specifically, scVI myeloid progenitor (CMP), granulocyte macrophage pro- employs a negative binomial model using raw counts, genitor (GMP), and control monocytes were extracted selecting 5000 highly variable genes to produce the latent from public whole-genome bisulfite sequencing (WGBS) variables. Defined cell-cycle phase-specific genes in the (GSE87197) [37]. Using GenomicRanges (v1.42.0) and Seurat package (v4.1.0) [43] were excluded from these to based on genomic location, the overlap of the hyper- reduce the dependence of clustering on cell-cycle effects. methylated DMPs observed in COVID-19 compared Data were subsequently analyzed using Scanpy (v1.9.1) with HD was determined with the Ensembl Regulatory [44] following the recommended standard practices. Regions from the hematopoietic precursors and mono- For quality control, genes expressed in fewer than three cytes. For this analysis, all DNA methylation data were cells, and cells with fewer than 200 genes or more than annotated with respect to the GRCh38 human genome 20% mitochondrial gene content, were removed prior reference. to downstream analysis. Data were normalized (scanpy. pp.normalize_per_cell, scaling factor = 10,000) and Single‑cell capture ­log2-transformed (scanpy.pp.log1p). For gene expression PBMCs from 10 ICU patients were used to generate sin- visualization (e.g., heat maps), data were further scaled gle-cell gel beads-in-emulsion (GEMs) (Additional file 2. (scanpy.pp.scale, maximum value = 10). Table S2). Cells were then washed three times and counted. For samples with low viability (< 90%), we performed Cell type clustering and annotation Ficoll separation in an Eppendorf tube to eliminate dead The resulting latent representation from the integrated cells and increase cell viability. For samples with greater datasets was used to compute the neighborhood graph than 90% viability, we filtered using a Flowmi strainer and (scanpy.pp.neighbors), then the Louvain clustering algo- counted the cells before loading into 10X chromium to rithm (scanpy.tl.louvain, resolution = 3) and Uniform generate single-cell GEMs, following the manufacturer’s Manifold Approximation and Projection (UMAP) visual- instructions. We loaded 50,000 cells per pool, including ization (scanpy.tl.umap) were employed. Cell type anno- a total of 4 patients per pool. Datasets from patients and tations were manually refined using literature-driven, HDs are available as h5ad files (https://​www.​COVID-​19cel​ cell-specific marker genes. Identified residual RBCs from latlas.​org/​index.​patie​nt.​html (Additional file 2. Table S2). incomplete PBMC isolation were excluded before further In parallel, genomic DNA was isolated from the same 10 analysis, as recommended [45]. lution (as described in [38]) using a Maxwell® 16 Blood PBMCs for genotyping and subsequent donor deconvo- Differential gene expression and transcription DNA Purification Kit from Promega following the manu- factor‑enrichment analysis facturer’s instructions. Differential gene expression between COVID-19 patients and healthy individuals (FDR < 0.05) was analyzed using scRNA‑seq cell type identification and annotation the limma package [46]. To predict transcription factor Single-cell transcriptome data from COVID-19 patients (TF) involvement in transcriptomic changes, we used were quantified and aligned using Cell Ranger (v3.1) DoRothEA (Discriminant Regulon Expression Analysis) with the GRCh38 genome concatenated to SARS-Cov-2 v2 tool [47]. Regulons with a confidence score of A–C genome as a reference. Thereafter, cells from pooled sam- were analyzed, and cases with p < 0.05 and a normalized ples were deconvolved and demultiplexed using Soupor- enrichment score (NES) of ± 2 were considered signifi- cell (v3.0) [39], yielding a genotype variant that allows cantly enriched. donor identity to be matched across different samples. This additionally enabled the removal of doublet cells that Cell–cell communication could not be explained by any single genotype. Scrublet Based on the differential expression analysis, Cell- (v0.2.3) [40] was subsequently employed to further filter PhoneDB [48] v3 (www.​CellP​honeDB.​org) was used to out other doublets based on computed doublet scores. infer changes in ligand/receptor interactions between the Specifically, Student’s t-test (p < 0.01) after Bonferroni identified cell types in COVID-19 versus HD. Specifically, correction was used within fine-grained sub-clustering instead of random shuffling, as used in the previously Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 6 of 22 described statistical method, differentially expressed (Miltenyi Biotec), CD366-PEvio615 (Miltenyi Biotec), genes (FDR < 0.05) were used to select interactions that CD62L-APC (Miltenyi Biotec), CD58-PE (Miltenyi Bio- were significantly enriched in either severe COVID-19 tec), CD191-PEvio770 (Miltenyi Biotec), CD52-PEvio615 patients or healthy individuals relative to the other group. (Miltenyi Biotec), CD48-APC (Miltenyi Biotec). Cells An interaction was considered enriched if at least one of were analyzed in a BD FACSCanto-II flow cytometer. the two partners (ligand or receptor) was differentially expressed, and if both partners were expressed by at least Statistical analysis 10% of the interacting cells. All statistical analyses were done with R v4.0.2. Box, bar, violin, bubble, and line plots were generated using func- Bisulfite pyrosequencing tions from the ggplot2 (v3.3.6) and ggpubr (v4.0) pack- EZ DNA Methylation-Gold kit (Zymo Research) was ages. Mean normalized DNA methylation values were used to BS-converted 500 ng of genomic DNA follow- compared using two-tailed test. Multivariate frequency ing the manufacturer’s instructions. BS-treated DNA was distributions were calculated using Fisher’s exact test. PCR-amplified using IMMOLASE DNA polymerase kit The levels of significance are indicated as: * p < 0.05, ** (Bioline). Primers used for the PCR were designed with p < 0.01, *** p < 0.001, and **** p < 0.0001. PyroMark Assay Design 2.0 software (Qiagen) (Addi- tional file 3. Table S3). PCR amplicons were pyrose- Results quenced with the PyroMark Q24 system and analyzed DNA methylome remodeling in peripheral blood with PyroMark Q48 Autoprep (Qiagen). monocytes of severe COVID‑19 patients To directly inspect epigenetic alterations in periph- Real‑time quantitative polymerase chain reaction eral blood monocytes in severe COVID-19, we isolated (RT‑qPCR) CD14 + CD15 − cells from 59 blood samples, com- The Transcriptor First Strand cDNA Synthesis Kit prising 48 severe COVID-19 patients and 11 healthy (Roche) was used to convert 250 ng of total RNA to donors (HDs), and performed DNA methylation profil- cDNA following the manufacturer’s instructions. RT- ing (Fig. 1A, Table 1, and Additional file 1. Table S1). For qPCR primers were designed with Primer3 software [49] cell sorting, we first separated live cells from debris, then (Additional file 3. Table S3). RT-qPCR reactions were extracted singlets and isolated CD14 + CD15 − cells to prepared with LightCycler 480 SYBR Green I Master avoid neutrophil contamination (Fig. 1B) [50]. Since we (Roche) according to the manufacturer’s instructions and selected CD14 + cells, the purification procedure only analyzed with a LightCycler 480 instrument (Roche). included classical (CM) (CD14 + CD16 −) and interme- diate monocytes (IM) (CD14 + CD16 +), excluding the Flow cytometry non-classical monocyte (NCM) (CD14lowCD16 +) sub- To study the surface cell markers on monocytes population, which in healthy individuals corresponds (CD14 +), PBMCs from the 10 patients used for single- to around 5% of the total monocyte compartment [51]. cell analysis and 10 HDs were defrosted and washed once Negative selection using CD15 was necessary, as there is with PBS. After blocking for non-specific binding with Fc a significant increase in the frequency of neutrophils in block (BD Pharmingen) for 5 min on ice, cells were incu- severe COVID-19 patients, as activated neutrophils are bated for 20 min on ice using staining buffer (PBS with not separated in the Ficoll step [52] (Additional file 4. 4% fetal bovine serum and 0.4% EDTA). Antibodies used Figure S1A-S1C). To confirm the purity of our mono- included the following: CD14-FitC (Miltenyi Biotec), cytes, we performed FACS analysis and obtained an CD85-PEvio770 (Miltenyi Biotec), CD172a-APC (Milte- average purity of 98% (example in Additional file 4 Fig- nyi Biotec), CD97-PEvio770 (Miltenyi Biotec), CD31-PE ure S1D). Studies in various other inflammatory diseases (See figure on next page.) Fig. 1 Analysis of DNA methylation in blood monocytes of severe COVID-19 patients. A Scheme depicting the cohort and workflow for monocyte purification of severe COVID-19 patients and controls and DNA methylation analysis. B Representative flow cytometry profile, indicating sorting gates used to purify monocytes from HD and COVID-19 patients’ peripheral blood. C Scaled DNA methylation (z-score) heatmap of differentially methylated positions (DMPs) between HDs (blue bar above) and COVID-19 patients (red bar above). Significant DMPs were obtained by applying a filter of FDR > 0.05 and a differential of beta value (Δß) > 0.15. A scale is shown on the right, in which blue and red indicate lower and higher levels of methylation, respectively. Clinical and treatment data of COVID-19 patients are represented above the heatmap. SOFA, IL-6 level, and days in the ICU scales are shown on the right of the panel D Principal component analysis (PCA) of the DMPs. HDs and severe COVID-19 patients are illustrated as blue and red dots, respectively. E Gene ontology of hypermethylated and hypomethylated DMPs. Selected significant functional categories (FDR < 0.05) are shown. F Bubble plot of TF motifs enriched on hypermethylated and hypomethylated DMPs. Bubbles are colored according to their TF family; their size corresponds to the FDR rank. G Box plot of individual DNA methylation values of CpG from hypermethylated and hypomethylated clusters (b-values), with the name of the closest gene and the position relative to the transcription start site Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 7 of 22 Fig. 1 (See legend on previous page.) Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 8 of 22 have shown that the proportions of monocytes can shift related to type I interferons (IFN) signaling and MHC between the three major subsets, i.e., CM, IM, and NCM. class II (Fig. 1E, bottom panel). For instance, it has been shown that severe COVID-19 Transcription factor (TF) binding motif enrichment patients feature reduced NCM and IM populations [53]. analysis, in 250-bp windows surrounding DMPs, revealed The analysis of monocyte subpopulations in our cohort overrepresentation of TFs of significance to the immune showed a significant increase in the CM population and response. The hypermethylated cluster CpGs displayed a decrease in the NCM population (Additional file 4. Fig- enrichment of binding motifs of IRFs and ETS TF fami- ure S1E-S1F). Since we purified CD14 + monocytes, our lies, which are linked to IFN changes (Fig. 1F, left panel). study only included CM and IM. Motifs of the bZIP TF family like AP-1, Jun, Fosl2, Fra1, We performed DNA methylation profiling of isolated and Fra2 were enriched in the hypomethylated cluster. monocytes and identified 2211 differentially methylated DMPs of the hypomethylated cluster were also enriched positions (DMPs) of CpGs in severe COVID-19 patients in motifs of the signal transducer factor and activator of compared with HDs (FDR < 0.05 and absolute Δß > 0.15). transcription factor (STAT) members STAT1 and STAT3. Of these, 1773 were hypermethylated (hypermethylated We also detected enrichment of the glucocorticoid cluster) and 438 were hypomethylated (hypomethylated response element (GRE) in the hypomethylated cluster cluster) (Fig. 1C and Additional file 5. Table S4). PCA of (Fig. 1F, right panel). Given these results, we hypoth- these DMPs showed that the two groups of monocytes esized that pharmacological treatment with glucocorti- (COVID-19 and HD) separated along the first principal coids (GCs) in severe COVID-19 patients in the intensive component axis (Fig. 1D). We obtained similar results care unit (ICU) might influence DNA methylation in when we included monocyte subpopulation propor- monocytes. To test this possibility, we performed limma tions as a covariate in the analysis (overlap, p < 0.0001) analysis and subsequent binding motif enrichment after (Additional file 6. Figure S2A). No significant differences separating COVID-19 patients into two groups, with and (FDR < 0.05) were observed within COVID-19 patients without GC treatment. Both groups of patients exhibited separated by their condition (obesity, hypertension, days significant enrichment of GRE motifs in the hypometh- admitted to the ICU, and exitus/death) or treatment with ylated cluster (Additional file 6. Figure S2E), suggest- dexamethasone (Additional file 1. Table S1). None of the ing that the endogenous production of GCs in severe abovementioned conditions was significantly correlated COVID-19 patients could participate in the hypometh- with the DNA methylation changes (Additional file 6. ylation through GRE. However, given the size of the Figure S2B). This was also apparent from the PCA show- cohort, we cannot rule out the possibility that pharmaco- ing the overlap of patients with different clinical param- logical treatment could also influence DNA methylation eters (Additional file 6. Figure S2C). changes and therefore remains as a potential confounder The analysis of the genomic functional features of the factor. DMPs in the hypermethylated and hypomethylated clus- Inspection of the individual genes within or in the ters (Additional file 6. Figure S2D) using public data from vicinity of the DMPs revealed several genes with func- monocytes [36] revealed an enrichment in promoters tions essential to the viral immune response. The list of and enhancers. This is consistent with their proposed relevant genes included IRF8, RUNX3, CD226, and CD83 roles for DNA methylation in regulatory elements [54]. in the hypermethylated cluster, and STAT1, FOXO3, Gene ontology analysis (GO) of the two DMP clusters IL1R1, and OAS1 in the hypomethylated cluster (Fig. 1G). revealed several functional categories associated with We validated these results using bisulfite pyrosequenc- the immune response to viral infection (Fig. 1E). In the ing in a new cohort of severe COVID-19 patients (Addi- hypermethylated cluster, we observed enrichment of cat- tional file 6. Figure S2F). Interestingly, these changes egories such as natural killer-mediated immunity, leuko- were also observed in mild COVID-19 patients (Addi- cyte migration, adaptive immune response, and positive tional file 6. Figure S2F). IRF8, IL1R1, and CD83 are regulation of interferon gamma production. We also associated with the IFN response. CD226 encodes a gly- observed hypermethylation in the MHC-II protein com- coprotein related to monocyte, NK, and T cell adhesion. plex that was related to antigen presentation. In addi- This glycoprotein has been shown to be involved in the tion, we found an enrichment of the positive regulation cytotoxicity of these cells and is known to be altered in of MAP kinase activity category (Fig. 1E, top panel). In COVID-19 patients [13]. STAT1 is associated with the the hypomethylated cluster, we observed enrichment of cytokine response, which, in turn, is related to IL1R1. The functional categories relevant to viral infection, including latter is the receptor of interleukin 1, which participates defense response to virus and negative regulation of viral in the inflammatory response and is strongly expressed genome replication. Importantly, the hypomethylated in severe COVID-19 patients [14]. OAS1 is induced by cluster also featured enrichment of functional categories interferons and activates latent RNase, causing viral RNA Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 9 of 22 degradation, which could be related to the identification this end, we first estimated the DNA methylation values of the category negative regulation of viral genome repli- of DMPs corresponding to the sepsis relative to the HD cation in the GO analysis. comparison from our previous sepsis study (accession number GSE138074) [27] using the data from the severe Monocytes from severe COVID‑19 patients display COVID-19 methylation dataset. Overall, we found sig- increased DNA methylation variability nificant enrichment in the hypermethylation and hypo- Overall, our DNA methylation analysis showed greater methylation clusters (Fig. 3A). We also calculated the heterogeneity (different variable positions, DVPs) in odds ratio of the overlap between these two datasets the profiles from COVID-19 patient monocytes than and found a strong enrichment of the hyper-DMPs in in those from HDs (Additional file 6. Figure S2G). We COVID-19 relative to those in sepsis (FDR ≤ 2.22·10−16) then examined the relationship between the DNA meth- and in the hypo-DMPs (FDR ≤ 2.22·10−16) (Fig. 3B). We ylation profiles and the Sequential Organ Failure Assess- also confirmed an enrichment in introns and depletion ment (SOFA) score, which is used in ICUs to calculate in promoters relative to the background when testing the organ damage. The score ranges from 0 to 24, with values genomic location of the DMPs common to both COVID- greater than 6 being associated with a significant increase 19 and sepsis (Fig. 3C and Figure S3A). DMPs located in in the risk of mortality [55]. Using Spearman’s correlation introns are often localized in enhancer regions involved coefficient to assess specific hypermethylated or hypo- in long-distance regulation [54]. methylated CpGs with SOFA, we identified 1375 CpG We then determined that the two datasets had 362 sites whose methylation levels positively correlated with hypermethylated and 92 hypomethylated CpGs in com- SOFA (increased methylation) (rho < 0.4 and p < 0.01) mon (Fig. 3D), corresponding to 51% of the total DMPs and 1497 CpG sites with an inverse correlation with of the sepsis patients (Additional file 8. Figure S3B). GO SOFA (decreased methylation) (rho <  − 0.4 and p < 0.01) analysis of the shared DMPs revealed significant enrich- (Fig. 2A and Additional file 7. Table S5). The mean nor- ment in functional terms related to host response, malization DNA methylation profiles of increased and including regulation of NK cells, inflammatory response, decreased methylation CpG sites were similar in patients and leukocyte chemotaxis (Additional file 8. Figure S3C). with low SOFA (< 6) and in healthy controls in an unsu- Shared hypermethylated CpGs were enriched in func- pervised representation but differed between the low tional categories related to cell signaling, such as the JAK- and high SOFA score groups (Fig. 2B). These results sug- STAT and MAPK pathways, that could be involved in the gest that changes in DNA methylation are concomitantly reduction of the inflammatory response and the IL15- exacerbated for higher SOFA scores, which is associated and IL12-mediated signaling pathways, which are related with bad prognosis. Several CpGs correlating with SOFA to cytokine production and Th1 proliferation (Fig. 3E, were associated with genes, such as IL17R, SOCS5, and left panel). Shared hypomethylated CpGs were enriched PCDHA5, that are involved in T cell-mediated inflam- in functional categories responsible for regulating the matory responses (Fig. 2C). Others, like FOXG1 and inflammatory response, such as negative regulation of CDC20B, are associated with DNA damage. GO analy- IL-1 production and positive regulation of macrophage sis revealed that changes in DNA methylation that are activation. In concordance with the hypermethylated concomitant with SOFA show an overrepresentation of cluster, we also observed negative regulation of IFNα terms associated with IFNγ, production of the molecular production (Fig. 3E, right panel). It is of note that severe mediator involved in inflammatory response, viral gene COVID-19-specific DMPs were enriched in functional expression, the B cell proliferation involved in immune categories related to virus infection, such as the defense response, and Th1 cell cytokine production (Fig. 2D). response to virus, and impairment of the antigen-pre- senting process, which seems to be specific to COVID-19 DNA methylation alterations in monocytes of severe infection [13, 23] (Additional file 8. Figure S3D). COVID‑19 patients significantly associate with those Inspection of TF binding motifs corresponding to the derived from patients with bacterial sepsis, myeloid DMPs shared between the two groups, separating the differentiation, and the influence of inflammatory shared hypermethylated and hypomethylated CpG sets cytokines revealed IRF family transcription factors like IRF1, IRF2, To better characterize the impact of DNA methylation IFR3, and IRF8 in the shared hypermethylated CpG set, changes in COVID-19, we compared the DMPs from which are well established regulators of the type I IFN severe COVID-19 patients with those obtained from system, being common in viral and bacterial infections monocytes derived from patients with bacterial sepsis [57]. We also detected enrichment of the ETS transcrip- in a previous study by our team [27], given that severe tion factors that are regulated by MAPK proteins, which COVID-19 can be considered a form of sepsis [56]. To were enriched in the GO analysis (Fig. 3F). In the shared Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 10 of 22 Fig. 2 DNA methylation changes in COVID-19 monocytes parallel organ damage. A Heatmap of severe COVID-19 patients with DNA methylation ordered by SOFA score, including all CpG-containing probes significantly correlated with the SOFA score (Spearman correlation coefficient rho > 0.4, p < 0.01). Clinical and treatment data of COVID-19 patients are shown above the heatmap. SOFA, IL-6 level, and days in the ICU scales are shown on the right of the panel B. Normalized methylation values from heatmap showing overall group methylation of HD. Patients with SOFA ≤ 6 are indicated as SOFA LOW; those with SOFA > 6 are indicated as SOFA HIGH. C DNA methylation levels (b-values) of selected individual CpGs (and closest genes) in hypermethylated and hypomethylated sets and their position relative to the transcription start site. D Gene ontology (GO) analysis of hypermethylated and hypomethylated DMPs, analyzed with the GREAT online tool, in which CpG annotation in the EPIC array was used as background. Statistical significance: * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001 hypomethylated set, we noted enrichment of STAT3 and several genes among the shared hypomethylated set, like TFs from bZIP AP-1, like Jun, and other bZIPs, like CEBP. CD163, SOCS1, and IL10, that have been associated with Interestingly, GRE was also present in the shared hypo- the acquisition of tolerogenic properties in monocytes methylated cluster (Fig. 3F). This suggests the influence [58] (Fig. 3G). of GC in the acquisition of aberrant methylation profiles In both infections, systemic inflammation could be in COVID-19 and sepsis. Individual genes associated responsible for part of the DNA methylation changes that with the COVID-19/sepsis shared hypermethylated and arise in monocytes. To address this possibility, we exam- hypomethylated CpG genes include type I IFN-related ined the DNA methylation levels of the hypomethylated genes, like IRF2, and others, such as IL1A and CCR2, that and hypermethylated CpGs of severe COVID-19 and are involved in inflammatory processes and monocyte sepsis patients in monocytes isolated from healthy donor chemotaxis, respectively (Fig. 3G). We also identified PBMCs that had been treated in vitro with inflammatory Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 11 of 22 cytokines like IFNα, IFNγ, and TNFα [26] (accession and compared them with those of 10 HDs from a pub- number GSE134425). This analysis revealed several sig- lic dataset [41] (Additional file 2. Table S2 and Addi- nificant changes following the trends for both COVID- tional file 9. Figure S4A-S4B). This analysis enabled us 19 and sepsis (Additional file 8. Figure S3E), suggesting to identify 24 cell populations based on specific mark- that these inflammatory cytokines, which are elevated ers (Fig. 4A and Additional file 9. Figure S4C-S4D), and in these patients, could influence the monocyte DNA thereby not only to determine the alterations in gene methylomes. expression in monocytes, but also to inspect altera- An alternative explanation for the observed changes in tions in additional immune cell subsets. Strikingly, severe COVID-19 monocyte methylomes could be that the monocyte fraction comprised solely CD14 + cells DNA methylation changes reflect alterations during mye- (CD14 mono: CD14) (Fig. 4B). loid/monocyte differentiation or the release of immature In the CD14 + monocyte cluster, we identified 10,440 or aberrant monocytes. This has been described in severe differentially expressed genes (DEGs) between COVID- COVID-19 cases [13, 59–62]. It is worth noting that 19 patients and HDs (Additional file 10. Table S6). The immature cells are also released from the bone marrow in top DEGs (based on the fold change (FC)) included pro- sepsis [63]. To test this hypothesis, we used public whole- inflammatory molecules (IL1B, CCL3), surface mark- genome bisulfite sequencing (WGBS) data (GSE87197) ers (CD163, CD63, AREG, CD74, S100A12, S100A12, of progenitor cells including HSC, MPP, CMP, and GMP S100A8, S100A9), and transcription factors (JUN, MAFB, cells and monocytes as references. We compared the NF-KB) (Fig. 4C). We observed upregulation of mono- 1773 hypermethylated CpGs based on their genomic cyte-derived cell markers like S100A12, S100A8, and location and obtained 1511 unique Ensembl Regions, S100A9. S100A8 is already known to contribute to the which grouped in two clusters. Cluster 1 showed low- cytokine storm in severe COVID-19 [41, 64]. Pro-inflam- level demethylation in monocytes compared with all matory genes like IL1B of IRF1 were downregulated, as hematopoietic precursor cell types, whereas cluster 2 well as HLA genes, in agreement with previous stud- showed clear demethylation in monocytes (Fig. 3H). ies, suggesting decreased antigen presentation in severe These results are compatible with the possibility that a COVID-19 patients. Finally, we observed downregulation proportion of the DMPs in severe COVID-19 result from of the NF-κB inhibitor zeta-encoding gene NFKBIZ, con- aberrant myeloid differentiation or the release of imma- sistent with activation of this pro-inflammatory pathway ture monocytes, which display higher methylation levels, [65]. Since type I IFNs are essential for antiviral immu- and are not demethylated to the extent they are during nity, and the DNA methylation analysis had indicated normal differentiation. the potential occurrence of epigenetic alterations in IFN- stimulated genes (ISGs), we checked the expression levels Aberrant DNA methylation is associated with changes of genes regulated by type I IFNs and found downregula- in gene expression of COVID‑19 patient monocytes tion of several ISGs, such as STAT1, BST2, PTPN6, and To study the relationship between the DNA methyla- IRF1 (Additional file 11. Figure S5A). In addition, given tion changes and aberrant gene expression of mono- that some of the observed DNA methylation changes cytes derived from severe COVID-19 patients, we were associated with genes involved in antigen presen- obtained single-cell (sc) RNA-seq data of peripheral tation, we inspected HLA genes in our expression data blood mononuclear cells (PBMCs) from 10 additional and found this gene set to be significantly downregulated, severe COVID-19 patients from the same hospital (See figure on next page.) Fig. 3 Comparative analysis of DNA methylation in blood monocytes of severe COVID-19 and bacterial sepsis patients. A Violin plot representing the mean methylation state of the DMPs found in the comparison between HDs and sepsis patients with b-values obtained from severe COVID-19 patients. B Fisher’s exact test showing the odds ratio ± 95% confidence interval of the overlap between DMPs found in monocytes from bacterial sepsis patients and DMPs in monocytes from COVID-19 patients. C Proportions of the genomic locations (in relation to genes) of DMPs in COVID-19 and sepsis; Bg., background, EPIC probes. D Venn diagram of the overlap of COVID-19 DMPs identified by the comparison of HDs and severe COVID-19 patients with DMPs identified by the comparison between HDs and sepsis patients, separating hypermethylated and hypomethylated DMPs. E Gene ontology analysis of hypermethylated and hypomethylated overlapping DMPs identified in the previous comparison. Selected significant categories (p < 0.05) are shown. F TF binding motif analysis of shared hypermethylated and hypomethylated DMPs comparing HDs and COVID-19 patients, and by HDs and sepsis patients. The panel shows the fold change (FC), TF family. Boxes with black outlines indicate TF binding motifs with FDR < 0.05. G Box-plot showing the DNA methylation values of individual CpGs (together with the name of the closest gene and its position relative to the transcription start site) from the hypermethylated and hypomethylated clusters from both COVID-19 and sepsis. H Scaled DNA methylation heatmap of regions from the whole-genome bisulfite sequencing (WGBS) data of hematopoietic stem cells (HSCs), multipotent progenitors (MPPs), common myeloid progenitors (CMPs), and granulocyte macrophage progenitors (GMPs) that overlap with the genomic location of the 1772 hypermethylated DMPs identified in the COVID-19 vs. HDs comparison. Statistical significance: * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001 Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 12 of 22 Fig. 3 (See legend on previous page.) Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 13 of 22 consistent with dysfunction in antigen processing and expression levels showed enrichment of functional cat- presentation (Additional file 11. Figure S5B). egories like negative regulation of T cells, IFNα, and anti- GO analysis of both DEG sets revealed enrichment gen presentation (Additional file 11. Figure S5C-S5D). in functional terms coincident with those from DNA This analysis reinforced the relationship between DNA methylation analysis. We observed functional catego- methylation changes and expression changes related ries such as cytokine-mediated signaling, IL-12-medi- to the acquisition of a more tolerogenic phenotype ated signaling, negative regulation of T cell activation, in monocytes in COVID-19 patients. Some examples negative regulation of IFNγ production, and defense include IL10, a tolerogenic cytokine whose expression response to the virus in the upregulated cluster genes is increased in COVID-19, and NFKBIz, whose level of (Fig. 4D). Conversely, functional categories such as anti- expression is decreased (Fig. 3H). We validated these gen processing and presentation by MHC-I and MHC- results using bisulfite pyrosequencing and qRT-PCR with II and IFNγ-mediated signaling were enriched among a new cohort of severe COVID-19 patients (Additional the downregulated gene set (Fig. 4D). We then studied file 11. Figure S5E-F). The analysis also included mild TFs potentially involved in the transcriptomic changes COVID-19 that showed partial or total DNA methylation observed in COVID-19 monocytes, using Discrimi- changes to the extent seen in severe COVID-19 cases nant Regulon Expression Analysis (DoRothEA), and (Additional file 11. Figure S5E-F). found that MAF family members, GATA3, STAT4, and IRF4, were associated with upregulated genes in severe Potential relationship between transcriptional COVID-19 (Fig. 4E). Conversely, STAT6, STAT2, IRF2, and epigenetic reprogramming and altered immune cell– IRF3, and LYL1 were associated with downregulated cell communication genes (Fig. 4E). TF enrichment of upregulated and down- Given the overrepresentation of genes associated with regulated genes was also consistent with the results from cytokine activity, MHC class II-mediated antigen pres- DNA methylation analysis, in which binding motifs for entation among the observed DNA methylation, and several of these TFs were overrepresented among the gene expression alterations in severe COVID-19, we regions neighboring the DMPs. explored the potential correlation of these changes in We determined the significance of a negative correla- monocytes with their pattern of communication with tion between DMPs and the expression levels of their other immune cell types. To systematically analyze the closest genes (rho =  − 0.31; p = 9.8e − 16) (Fig. 4F). To effect of cell–cell communication on monocytes, we study the relationship between DNA methylation and used CellPhoneDB (www.​cellp​honedb.​org), a repository expression changes further, we performed Gene Set of ligands, receptors, and their interactions integrated Enrichment Analysis (GSEA) of the genes associated within a statistical framework that predicts enriched cel- with hypermethylated and hypomethylated CpG clus- lular interactions between two cell types using scRNA- ters. Genes associated with hypermethylated CpGs were seq datasets. This allowed us to infer potentially altered generally downregulated (NES = 1669; FDR = 0.0005), interactions between monocytes and other immune whereas those associated with hypomethylated CpGs cell subsets in severe COVID-19. In particular, we were upregulated (NES =  − 1187; FDR = 0.0596) in inspected cell–cell communication alterations between COVID-19 patients (Fig. 4G). GO analysis of genes CD14 + and CD4 + memory, CD4 + naïve, CD8 + mem- with an inverse relationship between methylation and ory, and CD8 + naïve T cells; B cell subsets including (See figure on next page.) Fig. 4 Correlation between DNA methylation and gene expression. A UMAP visualization showing the immune cell populations identified from Louvain clustering and cell-specific marker gene expression. B Dot plot representing the expression of selected marker genes identified in the cell population. The scale represents the mean gene expression level in the cell subset and the circle size represents the percentage of cells in the subset of expressing cells. C Heatmap representing differentially expressed genes (DEGs) with a l­og2(FC) > 0.6, above, and ­log2(FC) <  − 0.6, below. Genes overexpressed and downregulated in COVID-19 patients in relation to HDs are depicted in red and blue, respectively. D Gene ontology (GO) overrepresentation of GO Biological Process categories comprising the upregulated and downregulated DEGs. The odds ratios for each group and the − ­log2(FC) are shown. Selected significant categories (FDR < 0.05) are shown. E Discriminant Regulon Expression Analysis (DoRothEA) of COVID-19 severe patients compared with HDs. Normalized enrichment score (NES) and l­og2(FC) of transcription factor expression are depicted. F Correlation of average DNA methylation levels of DMPs with average gene expression of DEGs in the HDs vs. COVID-19 severe patients. ­Log2(FC) of expression is plotted on the y-axis, higher numbers representing a higher level of expression in COVID-19 and lower numbers a higher level of expression in HDs. DNA methylation is depicted on the x-axis as Δβ, lower numbers representing a lower level of methylation in COVID-19 monocytes, and higher numbers a lower level of methylation in HDs. Points are colored according to their genomic context. G Gene set enrichment analysis (GSEA) of HD vs. COVID-19, using hypomethylated-associated genes and hypermethylated-associated genes as genesets. The running enrichment score is represented, and the normalized enrichment score (NES) is shown above (FDR < 0.01). H Representation of individual DNA methylation values of DMPs from the hypermethylated and hypomethylated clusters (beta values), the position in respect to the transcription start site, and the relative expression of the closely related DEGs. Statistical significance: * p < 0.05, ** p < 0.01, *** p < 0.0001, **** p < 0.00001 Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 14 of 22 Fig. 4 (See legend on previous page.) Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 15 of 22 Fig. 5 Cell–cell communication analysis. Dot plot of selected receptor/ligand pair (A) and ligand/receptor (B) interactions between CD14 + monocytes and other cell components in the COVID-19 patient group. Gene expression is indicated as ­log2(FC) for differentially expressed genes (FDR < 0.05), which, in both cases (A and B), are the molecules presented on the left. The percentage expression of the differentially expressed genes in each cell type is indicated by the circle size. Molecules shown in blue are those expressed in CD14 + monocytes. Molecules expressed in the immune cell partner are shown in red memory, naïve, and plasma B cells; natural killer cells LILRB1, LILRB2, and PECAM1 (CD31), the products of (NK ­CD56dim: NK ­CD56bright) (Fig. 5A,B). Our analysis which are involved in the inhibition of immune response, revealed 4483 ligand/receptor pairs, in which the expres- were downregulated in monocytes. Their corresponding sion levels of ligands and receptors of CD14 + and/or ligand-encoding genes, CD99, HLA-F, and CD38, were interacting partners in the aforementioned cell types expressed in all the analyzed cell types, except for CD38, were significantly different between severe COVID-19 which is only expressed in NK and plasma B cells. Addi- patients and HDs, suggesting changes in the interaction tionally, the gene encoding for receptor LAIR1, which of the corresponding immune cells (Additional file 12. inhibits IL-2 expression, was upregulated in mono- Table S7). The aberrant levels of the proteins encoded by cytes [66], which might influence the interaction with these genes in monocytes were validated by flow cytom- cells expressing its corresponding ligand, i.e., plasma B etry (Figure S5G), supporting a potential impact on cell– cells and monocytes. Our analysis also revealed changes cell communication. in the expression of TNF receptor genes (TNFRSF14, Figure 5A illustrates the significant ligand-receptor TNFRSF1B, TNFRSF1A) in monocytes, which could interactions that may be affected when the expression of affect the interaction with T cells through the ligands receptor in monocytes is altered, revealing their poten- encoded by TNF and LTA. This is compatible with the tial impact on other cell types. In general, there was a possibility that TNF-associated DNA methylation altera- high frequency of interactions involving different types tions in monocytes could arise from altered interactions of NK cells, consistent with the terms observed in the with T cells through these ligand-receptor pairs. We also GO analysis performed with DMPs (Fig. 1E). PILRA, noted downregulation of the receptor TNFRSF14, which Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 16 of 22 interacts with CD160 in NK cells. Some studies have analysis suggests that pro-inflammatory cytokines, the argued that CD160 is essential for NK-mediated IFNγ release of immature or aberrant monocytes, and spe- production [67], a conclusion that is consistent with the cific dysregulated immune cell–cell communication results obtained in our gene ontology analysis of the events may be responsible for some epigenetic changes. DNA methylation data. ADGRE5 (CD97) was downregu- To date, there have been very few DNA methylation lated in monocytes. This receptor interacts with CD55, studies addressing the involvement of COVID-19 DNA which is expressed in all the analyzed cell types. This methylation in regulating the angiotensin-converting interaction is involved in leukocyte migration [68]. The enzyme 2 (ACE2) type I membrane receptor gene [72], potential alteration of this interaction could be linked to which is present in arterial, lung type II alveolar cells, the observed hypomethylation of CpGs close to genes where it acts as a SARS-CoV-2 receptor. There is a sug- related to leukocyte migration (Fig. 1E, top). gestion that the host epigenome may represent a risk We also examined DEGs corresponding to ligands factor for COVID-19 infection. Very few studies have expressed in all immune cell types, whose corresponding reported alterations in DNA methylation in relation to receptors are expressed in monocytes, to identify poten- immune responses [73–75]. Our study aimed to explore tial cell–cell communication alterations that might affect the involvement of DNA methylation in relation to a monocytes (Fig. 5B). In general, we detected upregula- severe COVID-19 outcome in the myeloid compart- tion of ligands in regulatory T cells (Treg) and down- ment, which is directly related to systemic inflamma- regulation of ligands in plasma B cells. We also observed tion. We specifically studied monocytes because it is increased levels of CCL5 and CCL3, expressed in NK the cell type that undergoes the most dramatic tran- cells, that interact with the CD191 receptor (CCR1), and scriptomic reprogramming during COVID-19 infection whose inhibition potentially suppresses immune hyper- [13, 21, 23, 76]. In this regard, our study provides the activation in critical COVID-19 patients [69]. In the con- first instance of DNA methylome profiling in a specific text of antigen presentation, there was upregulation of immune cell type in COVID-19 patients. HLA-F from Treg and NKT, which interacts with LILRB1 Our data revealed that most DNA methylation in monocytes. Recent studies have associated LILRB1 changes in monocytes derived from severe COVID-19 with the development of tolerance [70]. Our analysis also patients occurred in genomic sites enriched in PU.1 revealed low levels of CD99, expressed in CD4 + mem- binding motifs, consistent with earlier studies show- ory and naïve T cells, Treg and memory B cells, and the ing its role as a pioneer TF directly recruiting TET2 receptor PILRA, which is expressed in monocytes. The and DNMT3b [77]. In our case, most DNA methylation opposite occurs with CD8 memory and naïve T cells changes occurred in genes related to cytokines, MHC and NK CD56(bright), which enhances T cell migration class II proteins, and IFN signaling. Similar results [71]. There was a similar trend between CD74 and the about the defective function of MHC-II molecules and receptor APP expressed in monocytes, which is involved activation of apoptosis pathways were obtained in sin- in antigen processing and presentation. This could be gle-cell atlas studies of PBMCs from severe COVID-19 related to the impaired antigen presentation previously patients [6, 21, 78, 79] and in sepsis [80, 81]. highlighted in our data. We found that DNA methylation changes in severe In brief, the potential alteration of cell–cell communi- COVID-19 patients share some features with sepsis, cation events, through increased or decreased levels of especially those associated with the expression of toler- ligands and receptors involving inflammatory cytokines, ogenic cytokines like IL-10 [82]. The acute phase of antigen presentation-related factors, and cell activation these infections suggests a dysregulated inflammatory regulators, in severe COVID-19 patients could affect host response, resulting in an imbalance between pro- downstream cell-signaling pathways and TFs and per- inflammatory and anti-inflammatory mediators [14]. haps influence DNA methylation profiles in monocytes, Some studies have suggested that viral components thereby perpetuating aberrant immune responses. induce STAT1 dysfunction and compensatory hyper- activation of STAT3 in SARS-CoV-2-infected cells [83]. Discussion We noted the involvement of kinases like JNK, and Our results reveal that peripheral blood monocytes earlier studies had shown that COVID-19 infection from severe COVID-19 patients display aberrant DNA activates the JNK and ERK pathways that end in the methylomes and transcriptomes associated with func- AP-1-dependent gene expression of pro-inflammatory tions related to IFN type I signaling and antigen pres- cytokines [84]. One of the most strongly affected TFs entation, among others. The changes are significantly is STAT2, together with STAT6, which could be linked associated with organ damage and with DNA methyla- to the aberrant IFN signaling in monocytes in COVID- tion changes occurring in bacterial sepsis. Finally, our 19 [83]. The presence of STAT2 downregulation also Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 17 of 22 suggests a deficiency in the ability to cross-present to The relationship between DNA methylation and gene CD8 + T cells [85]. expression is complex. DNA methylation patterns are We also identified GRE binding sites in association with cell-type-specific and are established during dynamic dif- DNA methylation changes. Generally, the glucocorticoid ferentiation events by site-specific remodeling at regula- receptor (GR) is activated when patients are treated with tory regions [97]. In general, methylation of CpGs located GC. However, we also noted significant GRE enrich- in gene promoters, first exons, and introns is negatively ment in patients who were not treated with GC, suggest- correlated with gene expression [98]. The analysis of our ing that endogenous production of GC in COVID-19 data shows that there is an inverse correlation between patients could regulate GR and affect DNA methyla- the CpG methylation changes and the expression lev- tion at its genomic binding sites. GC is also produced els of the closest genes. The comparison of the inferred endogenously in sepsis patients in whom cytokines like TFs associated with DNA methylation changes and gene IL-1β, TNFα, and IL-6 induce its production from the expression changes shows common factors like IRF2 and adrenal cortex using cholesterol as a substrate to reduce IRF3, which regulate downregulated genes and hyper- inflammatory responses [86, 87]. These cytokines were methylated CpGs. In this context, it is possible that hypomethylated and overexpressed in our dataset, con- reduced levels of IFN regulatory factor IRF3 or defective sistent with the results of other studies that have reported IRF7 function reduces the level of IFNα/β gene expres- increased levels in the serum of COVID-19 patients [88, sion, increasing the sensitivity to viral infection [12, 99]. 89]. GRE binding sites are enriched in the DMPs com- Finally, analysis of cell–cell communication has mon to COVID-19 and sepsis. GR is a nuclear receptor revealed potential relationships between DNA methyla- expressed in most cell types that can trigger the expres- tion changes and altered communication of monocytes sion of anti-inflammatory genes through direct DNA and other immune cells (e.g., T, plasma B and NK cells). binding. Furthermore, GRE represses the action of other Our data suggest the potential reduction of interac- inflammation-related TFs, including members of the tions between monocytes and NK cells through CD160, NF-KB and AP-1 families [90, 91], which are also known which mediates the antibody-dependent cell-mediated to be downregulated in our cohort. Taken together, our cytotoxicity that it is essential for IFNγ production [67]. results suggest the existence of a relationship between The potentially greater interaction between monocytes extracellular factors associated with the cytokine storm and Treg through multiple ligand and receptor pairs is occurring in severe COVID-19 and DNA methylation an interesting finding, since Tregs are immunosuppres- changes. Several studies have shown an increase in the sive cells responsible for maintaining immune homeosta- levels of inflammatory cytokines in severe COVID-19, sis [100]. In any case, the use of CellPhone DB is useful which may contribute to the severity of the disease [92]. for inferring cell–cell communications events; however, However, it is also possible that the DNA methyla- additional validation experiments would be necessary to tion changes are partly due to the release of immature validate interactions and activation of downstream sign- or altered monocytes from myelopoiesis, as reported aling pathways. for severe COVID-19 [13, 20, 93, 94] and sepsis [63]. In our study, we could not determine whether the Release of immature myeloid cells from the bone marrow observed DNA methylation alterations in COVID-19 in severe COVID-19 is reminiscent of emergency mye- were the cause or the consequence of the changes in lopoiesis [95]. This is a well-known phenomenon, char- gene expression. The analysis of mild COVID-19 cases, acterized by the mobilization of immature myeloid cells in which the DNA methylation and expression level of to restore functional immune cells, and by its contribu- a few genes showed differences in their similarities with tion to the dysfunction of innate immunity [96]. In fact, severe COVID-19 cases, suggests that there are cases a proportion of the hypermethylated CpGs in monocytes where expression changes might anticipate DNA meth- from severe COVID-10 patients overlap with regions that ylation changes. In any case, it is reasonable to propose become demethylated during myeloid differentiation. that some DNA methylation changes help perpetuate This suggests that part of the hypermethylated CpG sites dysregulated immune responses. in isolated peripheral blood CD14 + might be associated Some limitations of our study include the size of the with aberrantly differentiated monocytes released into cohort, and the unequal numbers of individuals admin- the bloodstream in severe COVID-19 patients. However, istered particular drugs in the different patient groups, the small numbers of CD34 + cells in the PBMC frac- which could have affected the COVID-19 data. How- tion of COVID-19 patients and the lack of CD14 + cells ever, despite these limitations, we found no significant in this subset suggest no interference with our results for differences among severe COVID-19 patients with CD14 + CD15- cells, isolated with our method. respect to the time they were admitted to the ICU or began to receive treatment. This suggests that DNA Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 18 of 22 methylation is quite a general occurrence in the con- and Projection; DoRothEA: Discriminant Regulon Expression Analysis; NES: Normalized enrichment score; RT-qPCR: Real-time quantitative Polymerase text of COVID-19. Another limitation concerns the cell Chain Reaction; CM: Classical monocytes; IM: Intermediate monocytes; NCM: population analyzed, since the method for monocyte Non-classical monocytes. isolation comprises two populations, CM and IM, one of which (CM) is expanded in the patient group. How- Supplementary Information ever, the analysis including the monocyte subsets as a The online version contains supplementary material available at https://​doi.​ covariate indicates that there are no major differences. org/​10.​1186/​s13073-​022-​01137-4. Finally, in the comparison with DNA methylation of Additional file 1: Table S1. Clinical features of severe COVID-19 patients progenitor cells, it is important to note that the DMPs and healthy donors (HD) included in the DNA methylation analysis. were overlapped with genomic regions, and not single- Additional file 2: Table S2. Clinical features of severe COVID-19 patients base data, and further analyses would be required. and healthy donors (HD) included in scRNAseq analysis. Future studies would benefit from having access to a Additional file 3: Table S3. Primer sequences. wider cohort in which it is possible to identify signifi- Additional file 4: Figure S1. Flow cytometry profiles of a representative cant links between alterations and drug treatments. sample for each group; HD (A) and COVID-19 (B) indicating the sorting Incorporating mild and asymptomatic cases would strategy and gates used in the study. C Boxplot representing the mean percentage of CD14+ and CD15+ cells in HD and patients. D Flow improve our ability to dissect drug- and severity-related cytometry profile of the CD14+CD15- purified monocytes using the same specificity in relation to DNA methylation changes. As gates used in the study. E HD and COVID-19 indicating the gates used for is the case for other medical conditions, the analysis 1of monocyte subtype analysis (classical monocytes, CM; intermediate mono‑ cytes, IM; and non-classical monocytes, NCM) (F) Boxplot representing the DNA methylation changes would be very likely to help mean percentage of CM, IM and NCM in HD and COVID-19 patients. predict disease severity, progression, and recovery. Additional file 5: Table S4. List of hypermethylated and hypomethylated genes in COVID-19 monocytes (FDR < 0.05; ∆β ≥ 0.2). Additional file 6: Figure S2. DNA methylation analysis in blood mono‑ Conclusions cytes of severe COVID-19 patients. A Venn diagram of the overlap of Our study provides unique insights into the epigenetic hyper- and hypomethylated DMPs identified with a subcluster of samples including (CM/IM_YES) or not (CM/IM_NO) the CM and IM percentage as a alterations of monocytes in severe COVID-19. We have covariable in the comparison between HD and severe COVID-19 patients shown that peripheral blood monocytes from severe (B) Heatmap representation of beta values of first 10 Principal compo‑ COVID-19 patients undergo changes in their DNA nent (PC) analysis that correlations PCs with different clinical variables. Numerical variables were correlated to PCs using Pearson correlation, methylomes, in parallel with changes in expression, whereas categorical variables were entered in a linear model together and that these significantly overlap with those found in with the PCs. C Principal component analysis (PCA) of the DMPs. The HDs patients with sepsis. We have also shown DNA meth- are illustrated in grey, and the severe COVID-19 patients are illustrated as blue and red in function of the clinical parameter or treatment with dexa‑ ylation changes are associated with organ dysfunc- methasone. D Enrichment analysis of different chromatin states for CpGs tion. Finally, our results suggest a relationship between sites corresponding the Hyper- and Hypomethylated clusters. The FDR is DNA methylation changes in COVID-19 patients and represented with the size of the bubble, as shown. The relative enrich‑ ment is represented as Odds Ratio. TssA, Active TSS; TssBiv, Bivalent/Poised changes that occur during myeloid differentiation TSS; BivFlnk, Flanking Bivalent TSS/Enh; EnhBiv, Bivalent Enhancer; ReprPC, and others that can be induced by pro-inflammatory Repressed PolyComb; ReprPCWk, Weak Repressed PolyComb; Quies, Qui‑ cytokines. CellPhoneDB analysis also suggests that escent/Low; TssAFlnk, Flanking Active TSS; TxFlnk, Transcr. at gene 5’ and 3’; Tx, Strong transcription; TxWk, Weak transcription; EnhG, Genic enhancers; alterations in immune cell crosstalk can contribute to Enh, Enhancers; ZNF/Rpts, ZNF genes & repeats; Het, Heterochromatin. E transcriptional reprogramming in monocytes, which TF binding motif analysis of hypomethylated DMPs comparing patients involves dysregulation of interferon-related genes no treated with dexamethasone vs. HDs (DEX_no) and patients treated with dexamethasone vs. HDs (DEX_yes). The panel shows the fold change and genes associated with antigen presentation and (FC), TF family. Black outlined boxes indicate TF binding motifs with FDR chemotaxis. values < 0.05. F Box plot of individual DNA methylation values of CpG from the hypermethylated and hypomethylated clusters with the name of the closest gene and the position in respect to the transcription start site Abbreviations calculated using pyrosequencing in the validation cohort that include HD SARS-CoV-2: Severe acute respiratory syndrome coronavirus 2; COVID-19: and patients with mild and severe infection of COVID-19. G Volcano plot Coronavirus disease 2019; IL: Interleukin; GM-CSF: Granulocyte-macrophage showing the p value vs the variance ratio for HD and COVID-19 associated colony-stimulating factor; TNF-a: Tumor necrosis factor-alpha; IFN: Inter‑ differentially variable CpG positions (DVPs). DVPs were identified using the feron; IRF: Interferon regulatory factor; NK: Natural killer; DC: Dendritic cell; algorithm iEVORA. Statistical significance: * p < 0.05, ** p < 0.01, *** p < MHC-II: Class II major histocompatibility complex; 5mC: 5-Methylcytosine; 0.0001, **** p < 0.00001. TF: Transcription factor; DNMT: DNA methyltransferase; TET: Ten–eleven Additional file 7: Table S5. List of hypermethylated and hypomethylated translocation; scRNA-seq: Single-cell RNA sequencing; ICU: Intensive care CpGs related to SOFA (rho < 0.4/ rho > 0.4 and p < 0.01). unit; PBMC: Peripheral blood mononuclear cell; HD: Healthy donor; FBS: Fetal Additional file 8: Figure S3. DNA methylation comparative analysis bovine serum; BS: Bisulfite; FDR: False discovery rate; DVPs: Differentially vari‑ between blood monocytes of severe COVID-19 patients and bacterial able positions; PCA: Principal component analysis; PC: Principal component; sepsis patients (A) Proportions of the genomic locations (in relation to HSC: Hematopoietic stem cell; MPP: Multipotent progenitor; CMP: Common genes) of hyper- and hypomethylated DMPs in COVID-19 and sepsis; Bg., myeloid progenitor; GMP: Granulocyte macrophage progenitor; WGBS: background, EPIC probes. B Venn diagram of the overlap of COVID-19 Whole-genome bisulfite sequencing; UMAP: Uniform Manifold Approximation Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 19 of 22 Availability of data and materials DMPs identified by the comparison between HD and severe COVID-19 DNA methylation data for this publication have been deposited in the NCBI patients with DMPs identified by the comparison between HD and septic Gene Expression Omnibus and are accessible through GEO SuperSeries patients. Gene ontology (GO) analysis of the shared (C) and not shared (D) accession number GSE188573 [30]. Single-cell data generated in this study are DMPs from the previous representation. E Violin plot of the mean methyla‑ available as h5ad files at the European Genome-Phenome Archive (EGA) with tion status of the identified DMPs with b-values obtained from mono‑ Dataset ID EGAD00001007982 [101]. cytes derived from healthy donor PBMCs exposed in vitro for 4 days to interferon-alpha (IFNα) (100ng/mL), interferon gamma (IFNγ) (100ng/mL), tumor necrosis factor-alpha (TNFα) (10ng/mL) and untreated (Control) Declarations (n=3), [26]. Statistical significance: * p < 0.05, ** p < 0.01, *** p < 0.0001, **** p < 0.00001. Ethics approval and consent to participate This study was approved by the Clinical Research Ethics Committees of Hos‑ Additional file 9: Figure S4. A UMAP visualization showing the two pital Universitari Germans Trias i Pujol (PI-20–129) and Vall d’Hebron University groups identified from Louvain clustering. B Barplot representation of Hospital (PR(AG)282/2020), which adhered to the principles set out in the the proportion of the two study groups in the different cell type found. C WMA Declaration of Helsinki. All samples were managed in compliance with Barplot representation of proportion of each cell type in each study group. participants’ written informed consent to participate in the study. D Barplot representation of proportion of each cell type in each sample used for the scRNA-seq analysis. Consent for publication Additional file 10: Table S6. List of Differentially Expressed Genes (DEGs) Written informed consent for publication was provided by the participants. in CD14+ cells between severe COVID-19 patients and HD (FDR < 0.05). Additional file 11: Figure S5. Analysis of DEG in monocytes derived Competing interests from severe COVID-19 patients. A Heatmap representation of expression The authors declare that they have no competing interests. levels of IFN genes (ligands and receptor) in HD and COVID-19 CD14+ monocytes. Gene ontology (GO) analysis of the upregulated. B Dot plot of Author details 1 major histocompatibility complex (MHC) genes. LogFc and percentage of Epigenetics and Immune Disease Group, Josep Carreras Research Institute expression is representeate. C and downregulated (D) DEG that present (IJC), 08916 Badalona, Barcelona, Spain. 2 Wellcome Sanger Institute, Wellcome a negative correlation with their close CpG. E Box plot of individual DNA Genome Campus, Hinxton CB10 1RQ, UK. 3 Department of Medical Genetics, methylation values of CpG from the hypermethylated and hypomethyl‑ University of Cambridge, Cambridge CB2 0QQ, UK. 4 Intensive Care Depart‑ ated clusters with the name of the closest gene and the position in ment, Vall d’Hebron University Hospital, Shock, Organ Dysfunction and Resus‑ respect to the transcription start site. Calculated using pyrosequencing in citation (SODIR) Research Group, Vall d’Hebron Research Institute (VHIR), the validation cohort that include HD and patients with mild and severe Universitat Autònoma de Barcelona, 08035 Barcelona, Spain. 5 Immunology infection of COVID-19. F Box plot of relative expression of individual genes Division, Vall d’Hebron University Hospital and Diagnostic Immunology performed by real-time quantitative polymerase chain reaction (RT-qPCR) Research Group, Vall d’Hebron Research Institute (VHIR), 08035 Barcelona, in the validation cohort that include HD and patients with mild and severe Spain. 6 Epigenetics in Inflammatory and Metabolic Diseases Laboratory, infection of COVID-19. G Box plot of mean florescence intensity (MFI) Health Science Center (HSC), East China Normal University (ECNU), Shang‑ of cell surface markers in CD14+ cells in the cohort used for single cell hai 200241, China. analysis. Statistical significance: * p < 0.05, ** p < 0.01, *** p < 0.0001, **** p < 0.00001. Received: 29 April 2022 Accepted: 7 November 2022 Additional file 12: Table S7. List of ligand-receptor pairs with different expression in COVID-19. References Acknowledgements 1. Lai C-C, Shih T-P, Ko W-C, Tang H-J, Hsueh P-R. Severe acute respiratory We thank the CERCA Programme/Generalitat de Catalunya and the Josep syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 Carreras Foundation for institutional support. We are very grateful to the (COVID-19): The epidemic and the challenges. Int J Antimicrob Agents. Cytometry Laboratory and the Scientific and Technological Centers, Univer‑ 2020;55:105924. sity of Barcelona (CCiT-UB), for their help with the gating strategy and cell 2. Schultze JL, Aschenbrenner AC. COVID-19 and the human innate sorting. immune system. Cell. 2021;184:1671–92. 3. Sette A, Crotty S. Adaptive immunity to SARS-CoV-2 and COVID-19. Cell. Authors’ contributions 2021;184:861–80. G.G.-T. and E.B. conceived and designed the study; G.G.-T., A.G.F.-B., and 4. Huang C, Wang Y, Li X, Ren L, Zhao J, Hu Y, et al. Clinical features of L.C. prepared and purified the samples; G.G.-T., O.M.-P., C.C.-F., and R.M.-F. patients infected with 2019 novel coronavirus in Wuhan, China. Lancet. performed bioinformatic analyses; A.B. analyzed the single-cell datasets; 2020;395:497–506. A.R.-S., M.M.G., R.F., and J.C.R.-R. provided the patient samples and analyzed 5. Tay MZ, Poh CM, Rénia L, MacAry PA, Ng LFP. The trinity of COVID-19: immu‑ the clinical data; G.G.-T., J.R.-U., R.V.-T., and E.B. analyzed and interpreted the nity, inflammation and intervention. Nat Rev Immunol. 2020;20:363–74. data; G.G.-T. and E.B. wrote the manuscript; all authors read and approved the 6. Zhang J-Y, Wang X-M, Xing X, Xu Z, Zhang C, Song J-W, et al. Single-cell final manuscript. landscape of immunological responses in patients with COVID-19. Nat Immunol. 2020;21:1107–18. Funding 7. Zhang C, Shi L, Wang F-S. Liver injury in COVID-19: management and This study has been funded by R+D+i project of the Spanish Ministry of challenges. Lancet Gastroenterol Hepatol. 2020;5:428–30. Science and Innovation (grant number PID2020-117212RB-I00/ MICIN/ 8. Zhou Y, Fu B, Zheng X, Wang D, Zhao C, Qi Y, et al. Pathogenic T-cells AEI/10.13039/501100011033). We also thank the Chan Zuckerberg Initiative and inflammatory monocytes incite inflammatory storms in severe (grant 2020–216799) and Wellcome Sanger core funding (WT206194). COVID-19 patients. Natl Sci Rev. 2020;7:998–1002. This publication has also been supported by the Unstoppable campaign 9. Henry BM, de Oliveira MHS, Benoit S, Plebani M, Lippi G. Hematologic, of the Josep Carreras Leukaemia Foundation. Additional funding comes biochemical and immune biomarker abnormalities associated with from Project PI18/00346 (Instituto de Salud Carlos III and co-funded by severe illness and mortality in coronavirus disease 2019 (COVID-19): a European Union (ERDF/ESF ). A.B. received additional support from a Gates meta-analysis. Clin Chem Lab Med. 2020;58:1021–8. Cambridge Scholarship. 10. Lee JS, Park S, Jeong HW, Ahn JY, Choi SJ, Lee H, et al. Immunophe‑ notyping of COVID-19 and influenza highlights the role of type Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 20 of 22 I interferons in development of severe COVID-19. Sci Immunol. 32. Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, 2020;5:eabd1554. Hansen KD, et al. Minfi: a flexible and comprehensive Bioconductor 11. Gustine JN, Jones D. Immunopathology of hyperinflammation in package for the analysis of Infinium DNA methylation microarrays. COVID-19. Am J Pathol. 2021;191:4–17. Bioinformatics. 2014;30:1363–9. 12. Sato M, Suemori H, Hata N, Asagiri M, Ogasawara K, Nakao K, et al. 33. Carvalho BS, Irizarry RA. A framework for oligonucleotide microarray Distinct and essential roles of transcription factors IRF-3 and IRF-7 in preprocessing. Bioinformatics. 2010;26:2363–7. response to viruses for IFN-alpha/beta gene induction. Immunity. 34. Teschendorff AE, Gao Y, Jones A, Ruebner M, Beckmann MW, Wachter 2000;13:539–48. DL, et al. DNA methylation outliers in normal breast tissue identify field 13. Schulte-Schrepping J, Reusch N, Paclik D, Baßler K, Schlickeiser S, Zhang defects that are enriched in cancer. Nat Commun. 2016;7:10478. B, et al. Severe COVID-19 is marked by a dysregulated myeloid cell 35. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple compartment. Cell. 2020;182:1419-1440.e23. combinations of lineage-determining transcription factors prime cis- 14. Hadjadj J, Yatim N, Barnabei L, Corneau A, Boussier J, Smith N, et al. regulatory elements required for macrophage and B cell identities. Mol Impaired type I interferon activity and inflammatory responses in Cell. 2010;38:576–89. severe COVID-19 patients. Science. 2020;369:718–24. 36. Bernstein BE, Stamatoyannopoulos JA, Costello JF, Ren B, Milosavljevic 15. Blanco-Melo D, Nilsson-Payant BE, Liu W-C, Uhl S, Hoagland D, Møller R, A, Meissner A, et al. The NIH Roadmap Epigenomics Mapping Consor‑ et al. Imbalanced host response to SARS-CoV-2 drives development of tium. Nat Biotechnol. 2010;28:1045–8. COVID-19. Cell. 2020;181:1036-1045.e9. 37. Farlik M, Halbritter F, Müller F, Choudry FA, Ebert P, Klughammer J, et al. 16. Trouillet-Assant S, Viel S, Gaymard A, Pons S, Richard J-C, Perret M, et al. DNA methylation dynamics of human hematopoietic stem cell dif‑ Type I IFN immunoprofiling in COVID-19 patients. J Allergy Clin Immu‑ ferentiation. Cell Stem Cell. 2016;19:808–22. nol. 2020;146:206-208.e2. 38. Rodríguez-Ubreva J, Arutyunyan A, Bonder MJ, Del Pino-Molina L, Clark 17. Bastard P, Rosen LB, Zhang Q, Michailidis E, Hoffmann H-H, Zhang Y, SJ, de la Calle-Fabregat C, et al. Single-cell Atlas of common variable et al. Autoantibodies against type I IFNs in patients with life-threatening immunodeficiency shows germinal center-associated epigenetic COVID-19. Science. 2020;370:eabd4585. dysregulation in B-cell responses. Nat Commun. 2022;13:1779. 18. Zhang Q, Bastard P, Liu Z, Le Pen J, Moncada-Velez M, Chen J, et al. 39. Heaton H, Talman AM, Knights A, Imaz M, Gaffney DJ, Durbin R, et al. Inborn errors of type I IFN immunity in patients with life-threatening Souporcell: robust clustering of single-cell RNA-seq data by genotype COVID-19. Science. 2020;370:eabd4570. without reference genotypes. Nat Methods. 2020;17:615–20. 19. Arunachalam PS, Wimmers F, Mok CKP, Perera RAPM, Scott M, Hagan T, 40. Wolock SL, Lopez R, Klein AM. Scrublet: computational identification et al. Systems biological assessment of immunity to mild versus severe of cell doublets in single-cell transcriptomic data. Cell Syst. 2019;8:281- COVID-19 infection in humans. Science. 2020;369:1210–20. 291.e9. 20. Silvin A, Chapuis N, Dunsmore G, Goubet A-G, Dubuisson A, Derosa 41. Stephenson E, Reynolds G, Botting RA, Calero-Nieto FJ, Morgan MD, L, et al. Elevated Calprotectin and Abnormal Myeloid Cell Subsets Tuong ZK, et al. Single-cell multi-omics analysis of the immune Discriminate Severe from Mild COVID-19. Cell. 2020;182:1401-1418.e18. response in COVID-19. Nat Med. 2021;27:904–16. 21. Guo C, Li B, Ma H, Wang X, Cai P, Yu Q, et al. Single-cell analysis of 42. Lopez R, Regier J, Cole MB, Jordan MI, Yosef N. Deep generative mod‑ two severe COVID-19 patients reveals a monocyte-associated and eling for single-cell transcriptomics. Nat Methods. 2018;15:1053–8. tocilizumab-responding cytokine storm. Nat Commun. 2020;11:3924. 43. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell 22. Wilk AJ, Rustagi A, Zhao NQ, Roque J, Martínez-Colón GJ, McKechnie JL, transcriptomic data across different conditions, technologies, and spe‑ et al. A single-cell atlas of the peripheral immune response in patients cies. Nat Biotechnol. 2018;36:411–20. with severe COVID-19. Nat Med. 2020;26:1070–6. 44. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expres‑ 23. Saichi M, Ladjemi MZ, Korniotis S, Rousseau C, Ait Hamou Z, Massenet- sion data analysis. Genome Biol. 2018;19:15. Regad L, et al. Single-cell RNA sequencing of blood antigen-presenting 45. Yao C, Bora SA, Chen P, Goodridge HS, Gharib SA. Sample processing cells in severe COVID-19 reveals multi-process defects in antiviral and single cell RNA-sequencing of peripheral blood immune cells from immunity. Nat Cell Biol. 2021;23:538–51. COVID-19 patients. STAR Protoc. 2021;2:100582. 24. Bader SM, Cooney JP, Pellegrini M, Doerflinger M. Programmed cell 46. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers death: the pathways to severe COVID-19? Biochem J. 2022;479:609–28. differential expression analyses for RNA-sequencing and microarray 25. Álvarez-Errico D, Vento-Tormo R, Sieweke M, Ballestar E. Epigenetic studies. Nucleic Acids Res. 2015;43:e47. control of myeloid cell differentiation, identity and function. Nat Rev 47. Garcia-Alonso L, Holland CH, Ibrahim MM, Turei D, Saez-Rodriguez J. Immunol. 2015;15:7–17. Benchmark and integration of resources for the estimation of human 26. Rodríguez-Ubreva J, de la Calle-Fabregat C, Li T, Ciudad L, Ballestar ML, transcription factor activities. Genome Res. 2019;29:1363–75. Català-Moll F, et al. Inflammatory cytokines shape a changing DNA 48. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. Cell methylome in monocytes mirroring disease activity in rheumatoid PhoneDB: inferring cell–cell communication from combined expres‑ arthritis. Ann Rheum Dis. 2019;78:1505–16. sion of multi-subunit ligand–receptor complexes. Nat Protoc Springer, 27. Lorente-Sorolla C, Garcia-Gomez A, Català-Moll F, Toledano V, Ciudad L, US. 2020;15:1484–506. Avendaño-Ortiz J, et al. Inflammatory cytokines and organ dysfunction 49. Koressaar T, Remm M. Enhancements and modifications of primer associate with the aberrant DNA methylome of monocytes in sepsis. design program Primer3. Bioinformatics. 2007;23:1289–91. Genome Med. 2019;11:66. 50. Agashe C, Chiang D, Grishin A, Masilamani M, Jones SM, Wood RA, et al. 28. Luo C, Hajkova P, Ecker JR. Dynamic DNA methylation: in the right place Impact of granulocyte contamination on PBMC integrity of shipped at the right time. Science. 2018;361:1336–40. blood samples: Implications for multi-center studies monitoring regula‑ 29. Deaton AM, Bird A. CpG islands and the regulation of transcription. tory T cells. J Immunol Methods. 2017;449:23–7. Genes Dev. 2011;25:1010–22. 51. Kapellos TS, Bonaguro L, Gemünd I, Reusch N, Saglam A, Hinkley ER, 30. Godoy-Tena G, Barmada A, Morante-Palacios O, de la Calle-Fabregat et al. Human monocyte subsets and phenotypes in major chronic C, Ciudad L, Ruiz-Sanmartín A, Martínez-Gallo M, Frank-Bertoncelj M, inflammatory diseases. Front Immunol. 2019;10:2035. Ferrer R, Ruiz-Rodriguez JC, Rodríguez-Ubreva J, Vento-Tormo R BE. 52. Kuhns DB, Priel DAL, Chu J, Zarember KA. Isolation and func‑ Epigenetic and transcriptomic remodeling in monocytes of severe tional analysis of human neutrophils. Curr Protoc Immunol. COVID-19 patients mirrors systemic changes in cytokines and immune 2015;111:7.23.1–7.23.16. cell crosstalk. Datasets. Gene Expression Omnibus. 2022. Available from: 53. Gatti A, Radrizzani D, Viganò P, Mazzone A, Brando B. Decrease of non- https://​www.​ncbi.​nlm.​nih.​gov/​geo/​query/​acc.​cgi?​acc=​GSE18​8573. classical and intermediate monocyte subsets in severe acute SARS- 31. Morante-Palacios O, Ballestar E. shinyÉPICo: a graphical pipeline to CoV-2 infection. Cytometry A. 2020;97:887–90. analyze Illumina DNA methylation arrays. Bioinformatics. 2021;37:257–9 54. Jones PA. Functions of DNA methylation: islands, start sites, gene bod‑ ies and beyond. Nat Rev Genet. 2012;13:484–92. Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 21 of 22 55. Ribas VJ, Lopez JC, Ruiz-Rodriguez JC, Ruiz-Sanmartin A, Rello J, Vellido 76. Yao C, Bora SA, Parimon T, Zaman T, Friedman OA, Palatinus JA, et al. A, On the use of decision trees for ICU outcome prediction in sepsis Cell-type-specific immune dysregulation in severely ill COVID-19 patients treated with statins. IEEE Symp Comput Intell Data Min. IEEE. patients. Cell Rep. 2021;34:108590. 2011;2011:37–43. 77. de la Rica L, Rodríguez-Ubreva J, García M, Islam ABMMK, Urquiza JM, 56. López-Collazo E, Avendaño-Ortiz J, Martín-Quirós A, Aguirre LA. Hernando H, et al. PU.1 target genes undergo Tet2-coupled demethyla‑ Immune response and COVID-19: a mirror image of sepsis. Int J Biol Sci. tion and DNMT3b-mediated methylation in monocyte-to-osteoclast 2020;16:2479–89. differentiation. Genome Biol. 2013;14:1–21. 57. Pietras EM, Saha SK, Cheng G. The interferon response to bacterial and 78. Xu G, Qi F, Li H, Yang Q, Wang H, Wang X, et al. The differential immune viral infections. J Endotoxin Res. 2006;12:246–50. responses to COVID-19 in peripheral and lung revealed by single-cell 58. Liu Q, Ou Q, Chen H, Gao Y, Liu Y, Xu Y, et al. Differential expression and RNA sequencing. Cell Discov. 2020;6:73. predictive value of monocyte scavenger receptor CD163 in popula‑ 79. Villani A-C, Satija R, Reynolds G, Sarkizova S, Shekhar K, Fletcher J, et al. tions with different tuberculosis infection statuses. BMC Infect Dis. Single-cell RNA-seq reveals new types of human blood dendritic cells, 2019;19:1006. monocytes, and progenitors. Science. 2017;356:eaah4573. 59. Townsend L, Dyer AH, Naughton A, Imangaliyev S, Dunne J, Kiersey R, 80. Delano MJ, Ward PA. Sepsis-induced immune dysfunction: can immune et al. Severe COVID-19 is characterised by inflammation and immature therapies reduce mortality? J Clin Invest. 2016;126:23–31. myeloid cells early in disease progression. Heliyon. 2022;8:e09230. 81. Venet F, Demaret J, Gossez M, Monneret G. Myeloid cells in sepsis- 60. Gibellini L, De Biasi S, Paolini A, Borella R, Boraldi F, Mattioli M, et al. acquired immunodeficiency. Ann N Y Acad Sci. 2021;1499:3–17. Altered bioenergetics and mitochondrial dysfunction of monocytes in 82. Boks MA, Kager-Groenland JR, Haasjes MSP, Zwaginga JJ, van Ham SM, patients with COVID-19 pneumonia. EMBO Mol Med. 2020;12:e13001. ten Brinke A. IL-10-generated tolerogenic dendritic cells are optimal for 61. Brauns E, Azouz A, Grimaldi D, Xiao H, Thomas S, Nguyen M, et al. functional regulatory T cell induction–a comparative study of human Functional reprogramming of monocytes in patients with acute and clinical-applicable DC. Clin Immunol. 2012;142:332–42. convalescent severe COVID-19. JCI insight. 2022;7. 83. Matsuyama T, Kubli SP, Yoshinaga SK, Pfeffer K, Mak TW. An aber‑ 62. Kvedaraite E, Hertwig L, Sinha I, Ponzetta A, Hed Myrberg I, Lourda rant STAT pathway is central to COVID-19. Cell Death Differ. M, et al. Major alterations in the mononuclear phagocyte land‑ 2020;27:3209–25. scape associated with COVID-19 severity. Proc Natl Acad Sci U S A. 84. Varshney B, Agnihothram S, Agnihotram S, Tan Y-J, Baric R, Lal SK. SARS 2021;118:e2018587118. coronavirus 3b accessory protein modulates transcriptional activity of 63. Reyes M, Filbin MR, Bhattacharyya RP, Sonny A, Mehta A, Billman K, et al. RUNX1b. PLoS ONE. 2012;7:e29542. Plasma from patients with bacterial sepsis or severe COVID-19 induces 85. Xu J, Lee MH, Chakhtoura M, Green BL, Kotredes KP, Chain RW, et al. suppressive myeloid cell production from hematopoietic progenitors STAT2 is required for TLR-induced murine dendritic cell activation and in vitro. Sci Transl Med. 2021;13:eabe9599. cross-presentation. J Immunol. 2016;197:326–36. 64. Ren X, Wen W, Fan X, Hou W, Su B, Cai P, et al. COVID-19 immune 86. Vandewalle J, Libert C. Glucocorticoids in sepsis: to be or not to be. features revealed by a large-scale single-cell transcriptome atlas. Cell. Front Immunol. 2020;11:1318. 2021;184:1895-1913.e19. 87. Dantzer R. Neuroimmune interactions: from the brain to the immune 65. Hariharan A, Hakeem AR, Radhakrishnan S, Reddy MS, Rela M. The role system and vice versa. Physiol Rev. 2018;98:477–504. and therapeutic potential of NF-kappa-B pathway in severe COVID-19 88. Conti P, Ronconi G, Caraffa A, Gallenga C, Ross R, Frydas I, et al. Induc‑ patients. Inflammopharmacology. 2021;29:91–100. tion of pro-inflammatory cytokines (IL-1 and IL-6) and lung inflamma‑ 66. Spolski R, Li P, Leonard WJ. Biology and regulation of IL-2: from molecu‑ tion by Coronavirus-19 (COVI-19 or SARS-CoV-2): anti-inflammatory lar mechanisms to human therapy. Nat Rev Immunol. 2018;18:648–59. strategies. J Biol Regul Homeost Agents. 2020;34:327–31. 67. Tu TC, Brown NK, Kim T-J, Wroblewska J, Yang X, Guo X, et al. 89. Belfer I, Buzas B, Hipp H, Dean M, Evans C, Lorincz I, et al. Haplotype CD160 is essential for NK-mediated IFN-γ production. J Exp Med. structure of inflammatory cytokines genes (IL1B, IL6 and TNF/LTA) in US 2015;212:415–29. Caucasians and African Americans. Genes Immun. 2004;5:505–12. 68. Safaee M, Clark AJ, Ivan ME, Oh MC, Bloch O, Sun MZ, et al. CD97 is a 90. Adcock IM, Nasuhara Y, Stevens DA, Barnes PJ. Ligand-induced differen‑ multifunctional leukocyte receptor with distinct roles in human cancers tiation of glucocorticoid receptor (GR) trans-repression and transac‑ (Review). Int J Oncol. 2013;43:1343–50. tivation: preferential targetting of NF-kappaB and lack of I-kappaB 69. Chua RL, Lukassen S, Trump S, Hennig BP, Wendisch D, Pott F, et al. involvement. Br J Pharmacol. 1999;127:1003–11. COVID-19 severity correlates with airway epithelium-immune cell inter‑ 91. Cain DW, Cidlowski JA. After 62 years of regulating immunity, dexa‑ actions identified by single-cell analysis. Nat Biotechnol. 2020;38:970–9. methasone meets COVID-19. Nat Rev Immunol. 2020;20:587–8. 70. Nowak I, Malinowski A, Barcz E, Wilczyński JR, Wagner M, Majorczyk 92. Vabret N, Britton GJ, Gruber C, Hegde S, Kim J, Kuksin M, et al. E, et al. Possible role of HLA-G, LILRB1 and KIR2DL4 gene polymor‑ Immunology of COVID-19: current state of the science. Immunity. phisms in spontaneous miscarriage. Arch Immunol Ther Exp (Warsz). 2020;52:910–41. 2016;64:505–14. 93. Dress RJ, Ginhoux F. Monocytes and macrophages in severe COVID-19 - 71. Dufour EM, Deroche A, Bae Y, Muller WA. CD99 is essential for leukocyte friend, foe or both? Immunol Cell Biol. 2021;99:561–4. diapedesis in vivo. Cell Commun Adhes. 2008;15:351–63. 94. Szabo PA, Dogra P, Gray JI, Wells SB, Connors TJ, Weisberg SP, et al. Longi‑ 72. Albini A, Di Guardo G, Noonan DM, Lombardo M. The SARS-CoV-2 tudinal profiling of respiratory and systemic immune responses reveals receptor, ACE-2, is expressed on many different cell types: implications myeloid cell-driven lung inflammation in severe COVID-19. Immunity. for ACE-inhibitor- and angiotensin II receptor blocker-based cardiovas‑ 2021;54:797-814.e6. cular therapies. Intern Emerg Med. 2020;15:759–66. 95. Schultze JL, Mass E, Schlitzer A. Emerging principles in myelopoiesis 73. Corley MJ, Pang APS, Dody K, Mudd PA, Patterson BK, Seethamraju H, at homeostasis and during infection and inflammation. Immunity. et al. Genome-wide DNA methylation profiling of peripheral blood 2019;50:288–301. reveals an epigenetic signature associated with severe COVID-19. J 96. Mann ER, Menon M, Knight SB, Konkel JE, Jagger C, Shaw TN, et al. Lon‑ Leukoc Biol. 2021;110:21–6. gitudinal immune profiling reveals key myeloid signatures associated 74. Balnis J, Madrid A, Hogan KJ, Drake LA, Chieng HC, Tiwari A, et al. with COVID-19. Sci Immunol. 2020;5:eabd6197. Blood DNA methylation and COVID-19 outcomes. Clin Epigenetics. 97. Ehrlich M, Lacey M. DNA methylation and differentiation: silencing, 2021;13:118. upregulation and modulation of gene expression. Epigenomics. 75. Castro de Moura M, Davalos V, Planas-Serra L, Alvarez-Errico D, Arribas C, 2013;5:553–68. Ruiz M, et al. Epigenome-wide association study of COVID-19 severity 98. Li L, Gao Y, Wu Q, Cheng ASL, Yip KY. New guidelines for DNA methy‑ with respiratory failure. EBioMedicine. 2021;66:103339. lome studies regarding 5-hydroxymethylcytosine for understanding transcriptional regulation. Genome Res. 2019;29:543–53. Godoy‑Tena et al. Genome Medicine (2022) 14:134 Page 22 of 22 99. Salman AA, Waheed MH, Ali-Abdulsahib AA, Atwan ZW. Low type I interferon response in COVID-19 patients: interferon response may be a potential treatment for COVID-19. Biomed reports. 2021;14:43. 100. Kondĕlková K, Vokurková D, Krejsek J, Borská L, Fiala Z, Ctirad A. Regulatory T cells (TREG) and their roles in immune system with respect to immunopathological disorders. Acta medica (Hradec Kral). 2010;53:73–7. 101. Single-cell multi-omics analysis of COVID-19 patients with pre-existing autoimmune diseases.Dataset. European Genome-Phenome Archive. 2022. Available from: https://​ega-​archi​ve.​org/​datas​ets/​EGAD0​00010​ 07982. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub‑ lished maps and institutional affiliations. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions