Functional reprogramming of monocytes in acute and convalescent severe COVID-19 patients

Severe COVID-19 disease is associated with dysregulation of the myeloid compartment during acute infection. Survivors frequently experience long-lasting sequelae but little is known about the eventual persistence of this immune alteration. Herein, we evaluated Toll-like receptor-induced cytokine responses in a cohort of mild to critical patients during acute or convalescent phases (n=97). In the acute phase, we observed impaired cytokine production by monocytes in the most severe patients. This capacity was globally restored in convalescent patients. Yet, we observed increased responsiveness to TLR1/2 ligation in patients that recovered from severe disease, indicating that these cells display distinct functional properties at the different stages of the disease. We identified a specific transcriptomic and epigenomic state in monocytes from acute severe patients that can account for their functional refractoriness. The molecular profile of monocytes from recovering patients was distinct and characterized by increased chromatin accessibility at AP-1 and MAF loci. These results demonstrate that severe COVID-19 infection has a profound impact on the differentiation status and function of circulating monocytes both during the acute and the convalescent phases in a completely distinct manner. This could have important implications for our understanding of short and long-term COVID19-related morbidity. In this study, we explore the functional features of circulating monocytes in acute and convalescent COVID-19 patients. Stimulation with Toll-like receptor ligands associated with bacterial or viral patterns established the functional impairment of circulating mononuclear phagocytes in acutely ill and severely affected patients. This was not observed in convalescent patients but we observed a distinct pattern of cytokine production. We further evaluated the transcriptional and epigenomic profiles of CD14 + monocytes from patients that experienced severe COVID-19. Taken together, these data support an underlying epigenetic basis for functional reprograming of monocytes during and after severe COVID-19.


Introduction
COVID-19, caused by the SARS-CoV-2 virus primarily infects the respiratory tract. There is a very broad spectrum of disease presentation. Some people who get COVID-19 have only mild symptoms. But for others, infection leads to pneumonia, respiratory failure, and, in some cases, death. The major complication of severe COVID-19 infection is acute respiratory distress syndrome (ARDS) presenting with dyspnoea and acute respiratory failure that requires mechanical ventilation. Severe COVID-19 appears to be associated with coagulopathy presenting as thrombosis in various organs and it is proposed that SARS-CoV-2 causes lesions to endothelial cells which then triggers an inflammation throughout the body and fuels ARDS (1,2). Multiple studies indicate that the mononuclear phagocyte system is strongly perturbed during acute infection and contributes to this hyperinflammatory complications (3)(4)(5). In patients with mild COVID-19, there is a specific increase of inflammatory monocytes that display a strong interferon (IFN)-stimulated gene signature. This was shown to be absent in severe disease. Instead, in these patients, there are signs of emergency myelopoiesis, marked by the occurrence of immunosuppressive neutrophils and HLA-DR Lo monocytes (6). Impaired type I IFN production and decreased expression of HLA-DR by plasmacytoid and conventional dendritic cells (DCs), respectively, was also reported in these patients (7,8).
In view of the novel nature of this disease, little is known about the long-term residual deficits of these patients. Initial assessment of COVID-19 survivors indicate that this disease could have multiple long-term effects leading to respiratory but also cardiovascular, renal or neurological sequelae (9,10). Here, in the light on the complex immune dysregulation that occurs during acute infection (11,12), we hypothesize that COVID-19 could have long-term impact on immune functions, such as the ones observed in patients that experienced bacterial sepsis (13). Indeed, immune profiling of SARS-CoV-2 recovered patients indicates persistent changes in the phenotype of innate (monocytes and granulocytes) but most importantly in adaptive (T, B, and NKT) immune cells (14)(15)(16).
In this study, we explore the functional features of circulating monocytes in acute and convalescent COVID-19 patients. Stimulation with Toll-like receptor ligands associated with bacterial or viral patterns established the functional impairment of circulating mononuclear phagocytes in acutely ill and severely affected patients. This was not observed in convalescent patients but we observed a distinct pattern of cytokine production. We further evaluated the transcriptional and epigenomic profiles of CD14 + monocytes from patients that experienced severe COVID-19. Taken together, these data support an underlying epigenetic basis for functional reprograming of monocytes during and after severe COVID-19.

Characteristics of the enrolled individuals
To assess the eventual persistence of immune dysregulation, we enrolled patients during acute infection and convalescence phase (2 to 9 months after the onset of the symptoms, secondarily divided into "early recovery" for days 56 to 120 and "late recovery" for days 160 to 268 post symptoms onset). Patients were categorized according to the severity of the disease, "mild disease" represents patients who weren't hospitalized and "severe disease" corresponds to patients who required hospitalization due to COVID-19 (in regular ward or intensive care).
Severe patients were more frequently male and older (Table 1). Most patients (89,7%) had a nasopharyngeal RT-PCR positive for SARS-CoV-2 at the hospital admission; otherwise, the diagnosis of COVID-19 was based on clinical and/or radiological features. During the acute phase, none of the mild patients died while death occurred in 38% of severe cases.
Comorbidities and other demographic and disease features are shown in extended data (Table S.1).

Strong alteration of TLR-induced cytokine production by classical monocytes in severe patients during acute infection.
To evaluate the functionality of the innate immune cells in COVID-19, we performed ex vivo stimulation of whole blood with different pathogen-associated molecular patterns (PAMPs): bacterial LPS, a Toll-like receptor (TLR)4 ligand, R848, a synthetic ligand for TLR7/8 and Pam3CSK4, a synthetic ligand for TLR1/2. Medium alone served as unstimulated condition.
We performed multiplex determination of cytokine levels in the supernatants and assessed cytokine production in monocytes at the single-cell level by multiparametric flow cytometry.
In acutely ill patients that did not require hospitalization, we observed few significant modulations of cytokine production in comparison to healthy controls: increased IFNβ levels in response to R848 and increased LPS-elicited IL10 and IL6 production. In sharp contrast, in severely ill patients, we observed strongly reduced levels of almost all cytokines in response to LPS or R848 stimulation with the exception of IFNβ which was significantly increased ( Fig.1). Of note, in comparison to other stimuli, Pam3CSK4 induced low levels of some cytokines (CCL2, IL6, IL10) but to a similar extent in the different groups. Because whole blood is a complex mixture of cells that can directly or indirectly produce these cytokines in response to TLR stimulation, we also evaluated cytokine production at the single-cell level using flow cytometry. Consistent with our Luminex data, the proportion of CD14 high monocytes expressing IL1β, IL6, TNFα and IL12/23 in response to LPS or R848 was found to be decreased in severe patients in comparison to healthy subjects or mild disease patients. In contrast to these conditions, stimulation with Pam3CSK4 was weaker but was comparable in the 3 groups ( Fig.2

.A).
To fully apprehend the extent of this functional alteration, we were interested in evaluating the capacity of each cell to produce more than 1 cytokine simultaneously (i.e.polyfunctionality). The capacity to produce 3 or 4 cytokines concurrently was impaired in severe patients during acute COVID-19 in response to stimulation by LPS or R848 but not by Pam3CSK4 (Fig.2B). This latter parameter was increased in monocytes from mild COVID-19 patients as compared to controls. Taken together, these experiments demonstrate that the ability of monocytes to mount an appropriate pro-inflammatory response is impaired in severe COVID-19 patients upon stimulation through TLR4 or TLR7/8. When considering the most informative immune parameters in all acute COVID-19 patients, we observed a very strong correlation between the different responses to LPS and R848 as compared to the same analysis in naive subjects. This indicates that during SARS-CoV2 infection, modulation of the capacity of monocytes to produce different cytokines is highly coordinated and that the effect is global (Fig. S.1A).
Based on these data for our 24 severe patients, we defined whether these immune parameters were correlated with clinical outcome. As shown in Fig. 2C, a low degree of cytokine production was significantly associated with a higher risk of death among hospitalized patients indicating that innate immune paralysis is a key feature of the most severe forms of COVID-19.
Functional modulation of monocyte function in severe COVID-19 convalescent patients.
In order to define whether this immune dysfunction persists after the resolution of the acute infectious episode, we prospectively recruited patients 56 to 268 days after the onset of the symptoms. We quantified cytokine production in whole blood assays. Except for increased IFNβ concentrations in previously infected and hospitalized subjects, cytokine levels were found to be comparable in non-infected controls and severe convalescent patients, indicating that the global dysfunctional state observed during severe COVID-19 infection is not longlasting (Fig.3). Of note, although it did not reach statistical significance in the whole group, we also observed high basal CCL2 production in a subset of these patients. In convalescent patients from mild disease, there was a trend for higher IL1β, IL6 and TNF production.
Interestingly, focusing on classical monocytes, patients that recovered from a severe disease showed an increased capability to produce IL1β, IL6 and TNFα upon engagement of TLR1/2. On the other hand, in response to LPS stimulation, we noted a reduced proportion of IL12/23 + cells in the same group of patients ( Fig.4A and B). This distinctive pattern appeared to concern mainly patients in the early phase of recovery (Fig.4C).
To visualize how all the different immune parameters globally vary across the clinical groups we used t-SNE, an analytical method that compares biological samples without taking into account sample classification (Fig.4D). We selected features that were significantly different between naïve controls, acute and convalescent patients (Fig.S.1B). This analysis shows an approximate 3-way separation with clusters enriched either for COVID-19 patients in the acute (cluster I) or convalescent stages (cluster III). Cluster I comprised a very clear subgroup of patients that were severely ill. In conclusion, alteration of TLR-responsiveness is a key feature of the most severe forms of COVID-19. Months after the recovery, the monocytes of these patients display a distinct pattern of cytokine production which suggests that the infectious event has induced a long-lasting functional reprograming.

Monocytes from acute and convalescent COVID-19 patients display distinctive transcriptomic profiles
To gain further insight into the molecular features of CD14 + monocytes in these severely affected COVID-19 patients during acute and convalescent stages we selected only severe patients from our cohort and age-and sex-matched controls and performed global transcriptional profiling. Characteristics of these individuals are shown in Table S.2. Since we observed differences in terms of cytokine production between early (56 to 120 post symptom onset) and late recovery stages (160 to 268 post symptoms onset) (Fig.4B), we analysed these groups separately. We observed a clear separation between samples from acute and healthy controls upon principal component analysis (Fig.5A). For samples from convalescent patients, those from early recovery stage formed a distinct cluster while those from late recovery stage were found to be embedded within the control group. In acute samples, we identified 339 statistically differentially expressed genes (DEG: 184 up-and 155 downregulated genes compared to controls, with a fold change (FC) > 2 and a false discovery rate (FdR) <0.05, Fig.5B). Consistent with their altered functional response and recent reports (17,18), we observed decreased expression of genes related to key immune pathways such as antigenic presentation, innate immune responses and MAPK and NF-κB signalling (JUNB, ATF3, NFkB2) in monocytes from patients suffering from severe COVID-19 (Fig.5B). We also observed increased expression of genes involved in key metabolic processes, including lipid metabolism (Fig.5C). Using the same criteria for samples from convalescent patients in the early recovery period, we identified 521 DEGs (318 up-and 203 downregulated genes compared to controls). Expression of multiple genes encoding chemokines was upregulated along with important intracellular immunomodulatory proteins and transcription factors (PPARG, FOSL1, MAFB, MAFF, ATF4, FOXO3). Hardly any DEGs were identified for the patients at the latter stage. Importantly, very few genes that were up or downregulated in the acute stage were found to be also modulated in the recovery groups (Fig.5B), indicating that the transcriptomic program induced by the acute infection is not persistent.
We also analyzed publicly available gene sets from scRNASeq data of acute COVID-19 patients. In acute phase, we observed strong enrichment for genes identified in monocytes from the most severely affected patients (19). Schulte-Schrepping et al. described several blood monocyte subsets that arise in mild and severe patients (6). In monocytes from acute samples, we observed strong enrichment for marker genes of clusters C2 which corresponds to HLA-DR lo CD163 hi cells that are the most abundant in early stages of the disease in severe patients. In monocytes from recovered patients, the strongest enrichment was seen for marker genes of cluster C0 (HLA-DR lo S100A + cells) that dominates later in the course of the acute disease and of cluster C1 (HLA-DR hi CD83 + cells) that corresponds to activated cells.
Monocytes from acute patients were also enriched for the MS1 gene signature derived from immature monocyte state in sepsis patients (20). In contrast, genes that were modulated in the recovery phase (i.e. several weeks after the sepsis) (21) were significantly up-or downregulated in early recovered patients, indicating potential common underlying mechanisms between these two clinical situations. Altogether, our transcriptomic data on CD14 + monocytes from acute severe COVID-19 patients indicate that these cells display an altered profile that could account for their decreased responsiveness to PAMPs. The distinctive transcriptional state of monocytes identified in recovery patients is probably not related to the same processes that occur in acutely ill patients and is not persistent in later stages of convalescence. This unique profile displays similarities with monocytes from patients that recovered from sepsis.

Monocytes from acute and convalescent severe COVID-19 patients display distinct profiles of chromatin accessibility
To identify the epigenetic determinants of these distinct transcriptional and functional programs, we mapped chromatin accessibility by ATAC-Seq. As expected, we observed extensive modifications in monocytes from acute COVID-19 patients as 1617 and 1582 regions were found to be significantly more or less accessible in acute patients vs controls, respectively (Fig.6A). Strikingly, we also observed important changes in early but not late convalescing patients. A low number of differentially accessible regions (DARs) identified were common to both comparisons (Fig.6B). We used the Binding and Expression Target Analysis (BETA) package (22) to predict the activating or repressive function of these DARs. For acute patients, regulatory regions that were more or less accessible were clearly associated with genes that were up or downregulated at the transcriptional level, respectively (Fig.6C). For early convalescent patients, we also observed strong association between the regions that are more accessible and genes that are activated in this group. In contrast, the regions that were found to be less accessible were not clearly associated with down-regulated genes, suggesting that the most important regulatory features in this case are linked to gene activation rather than repression. These observations strongly suggest that epigenetic imprinting is responsible for the transcriptional signatures identified in acute and early convalescing patients. For example, in acute patients, less accessible regions were found in the loci of IL1B or IL1R1 genes (Fig. 6D). In early convalescent patients, we observed increased accessibility at the loci of the same genes.
Next, we performed gene-ontology analysis using Genomic Regions Enrichment of Annotations Tool (GREAT) (23). Consistent with our transcriptomic data, we observed that regulatory regions that are less accessible in monocytes from acute patients were associated with genes involved in inflammatory response, lipid metabolism or cytokine-mediated signalling pathways (Fig 6E). We observed a mirror image in recovery patients with increased accessibility associated to genes involved in innate immune response and TLR signalling pathway but also in wound healing. We then scanned for binding motifs at the center of ATAC peaks located in these sets of regions. Analysis of putative TF site enrichment in patients vs control-specific regions indicated a strong and significant enrichment for distinct motifs in both groups (Fig.7A). Consistent with the decreased expression of NFKB2 and JUNB (Fig 7B), we observed an overrepresentation for AP1 (Jun/Fos), REL binding motifs in regions that were less accessible in monocytes from acute patients. We also identified IRF1 and STAT motifs in these regions, suggesting that multiple inflammatory modules are deactivated at this stage of the disease. In sharp contrast, in regions that were more accessible in convalescent patients, we observed a strong enrichment for AP-1 and MARE-containing motifs (Fig.7A). We identified both motifs in a substantial proportion of these regions (845/3458 regions, linked to 121 significantly up-regulated genes in recovery patients in comparison to controls, termed AP1/MAF transcriptional module) For example, such motifs were identified in the locus of the PPARG gene (Fig.7C). Of note, we observed increased levels of MAFB, MAFF but also FOSL1, a partner of Jun family members that displays an important role in orchestrating the expression of genes related to wound response, Toll-like receptor activation, and interleukin signaling of macrophages (24). Given that JUN proteins also form heterodimers with members of the MAF family (25), it is tempting to speculate that upregulation of MAFB expression in these cells redirects AP-1 complexes to these genomic targets thereby promoting local chromatin remodeling. Consistent with this hypothesis, we observed a strong correlation between MAFB and FOSL1 levels and expression of the AP1/MAF transcriptional module in the whole cohort (Fig.7D). These correlations were less obvious with other members of the MAF, JUN, FOS or BATF families. Finally, we observed strong enrichment for MAF-responsive genes in monocytes from convalescent patients in comparison to controls (Fig.7E) (26). MafBdependent genes rather than cMAF-dependent genes were enriched, reinforcing the potential contribution of this specific transcription factor. Taken together, these data support the notion that the transcriptional program and functional properties of monocytes in acute or convalescing severe COVID-19 patients have strong epigenetic determinants.

Discussion
Multiple reports now indicate that COVID-19 infection is associated with drastic changes in the myeloid compartment, particularly in patients with a severe course of disease (18,(27)(28)(29).
It has been shown that severe and fatal COVID-19 leads to accumulation of HLA-DR lo monocytes with potentially suppressive and dysfunctional features. Here, to evaluate the function of these cells more thoroughly, we first assessed the capacity of circulating cells from mild and severe patients to produce cytokines in response to prototypical TLR ligands. We observed that the ability of monocytes from severe COVID-19 patients to produce inflammatory cytokines was severely impaired. Of note, stimulation with the TLR1/2 ligand Pam3CSK4 was generally less affected than with other ligands and we observed even higher levels of IFNβ in LPS-or R848-stimulated whole blood cultures, indicating complex functional changes. Moreover, among severely ill patients, decreased ability of monocytes to produce several cytokines simultaneously was associated with a higher risk of mortality. Our Next, to define whether this hematopoietic reprograming during acute severe COVID-19 infection could have long-term impact on monocyte function, we recruited patients that had recovered from the disease. The production of cytokines in patients convalescing from mild or severe disease was globally comparable to that of SARS-CoV2 naive controls. Yet, at the single-cell level, we observed subtle changes including decreased production of IL-12 in response to LPS and increased responsiveness to TLR1/2 ligation in patients that recovered from severe disease. This was more obvious in the first few months after the disease, suggesting that the impact of COVID-19 infection/hospitalization on monocyte function fades after a longer period. In these early convalescent patients (2-3 months post-infection), we observed striking molecular profiles characterized by epigenomic reprogramming reminiscent of trained immunity. This term has been proposed for the persistent enhanced state of the innate immune response following exposure to certain infectious agents or vaccines, which may result in increased resistance to related or unrelated pathogens (31). However, similar processes may also result in hypo-responsiveness and be deleterious, for example in the context of chronic metabolic and inflammatory diseases such as liver cirrhosisa condition associated with increased susceptibility to infections (32). More recently, Wimmers et al showed that administration of seasonal influenza vaccine induced persistent epigenomic changes in myeloid cells leading to innate refractoriness associated with decreased AP-1 activity (33). Multiple mechanisms may account for these long-lasting effects on innate immune cells. In the context of BCG immunization, exposure to LPS or sepsis, it involves modulation of myeloid progenitors in the bone marrow (34)(35)(36). In the context of COVID-19, they could also be the consequence of persistent viral antigen expression (detected up to 4 months after the onset of the disease) (37) or result from the modulation of the adaptive compartment as the majority of SARS-CoV2 specific CD8 T cells acquire a terminally differentiated phenotype (38).
Single-cell epigenomic profiling identified distinct states among classical monocytes and indicated that "activated" or "trained" subsets were enriched in convalescing COVID-19 patients (39). In line with this report, we observed that the epigenomic state of monocytes from patients that recovered from severe COVID-19 was associated with heightened AP-1 and MAF activities. We uncovered the potential contribution of FOSL1 and MafB in this functional reprograming. FOSL1 has been shown to regulate pro-and anti-inflammatory cytokine expression in macrophages, modulating profibrotic responses (40) and promoting lung or joint inflammation (41,42). MafB is also a key regulator involved in functional programing of macrophages in the context of tissue imprinting (43), lipid/cholesterol-rich environments or wound healing (44). Of note, the balance between MAF and MafB in alveolar macrophages could shape the response to SARS-CoV2 (45).
This functional reprograming of monocytes in convalescing COVID-19 patients could have both beneficial and detrimental consequences. In particular, an emerging complication of COVID-19 infection is a prolonged period of lingering symptoms post-infection (10). An hyperinflammatory state seems to be a cardinal feature of this syndrome that is also associated with increased risk for sudden cardiovascular events (46). Hence, the epigenomic changes in circulating monocytes we describe here could participate to cytokine-driven endothelial dysfunction (47).
We are mindful of the limitations of the present study. It would be important to follow these patients longitudinally as they represent a very heterogeneous population. Hence, it is possible that more subtle variations, for example in mild patients, have been overlooked. In addition, the differences we identified here using bulk approaches probably reflect the changes in the most abundant subsets of monocytes. Furthermore, we cannot exclude a selection bias in convalescent patients that survived from an acute infection. Finally, multiple factors that could influence the function of these cells were not taken into account, including treatment modalities and co-morbidities.
In conclusion, it is now clear that immune responses are strongly influenced by past and present interactions of the host immune system with its environment (48). We showed here that severe COVID-19 infection has a profound impact on the differentiation status and function of circulating monocytes both during the acute and the convalescent phases in a completely distinct manner. This could have important implications for our understanding of short and long-term COVID19-related morbidity and mortality.

Patient recruitment
For this study, between June 2020 and January 2021, we recruited 97 patients during acute infection or at early and late recovery phases. Patients were categorized according to the severity of the disease, "mild disease" represents patients who weren't hospitalized and "severe disease" corresponds to patients who were hospitalized due to COVID-19 (in regular ward or intensive care). Among acute infections, 11 were mild and 24 severe. 19 patients recovered from mild and 43 from severe disease. As controls, we recruited 32 SARS-CoV2naïve individuals (negative nasopharyngeal RT-PCR and serology for SARS-CoV-2) among health-care workers and nursing home residents (49) that were age-and sex-matched to severe COVID-19 patients.

Blood collection
All blood draws were performed in the hospital by a trained phlebotomist. Peripheral blood was drawn via sterile venipuncture into sodium-heparin vacutainers. Blood samples were kept at room temperature (RT) and processed within 4 hours of the blood draw.
For the intracellular cytokine staining (ICS), 10 µg/ml of brefeldin A (final concentration -Sigma-Aldrich) was added to each well, and samples were incubated for 6h at 37°C in 5% CO2, then treated with 2 mM EDTA (final concentration) for 10 min at 37°C. The cells were collected and resuspended in BD FACS Lysing Solution, placed into fresh tubes, and stored at -80°C.
For multiplex analysis, samples were incubated with whole blood for 24h at 37°C in 5% CO2.
The supernatant was collected and stored at -80°C.
Flow cytometric data acquisition was performed on a Cytoflex LX Beckmann & Coulter and analyzed using FlowJo software (Fig.S.3).

Isolation of peripheral blood mononuclear cells (PBMCs) and sorting of CD14 + monocytes
PBMCs were isolated from sodium-heparinised whole blood by density gradient centrifugation using Lymphoprep TM (ProteoGenix) and Leuco-Sep tubes (Greiner) and PBMCs were frozen in FCS with 10% DMSO in liquid nitrogen. PBMCs were thawed, stained with trypan blue and counted to ensure a majority of live cells. CD14 + monocytes were isolated by fluorescenceactivated cells sorting (FACS) on a SONY SH800S Cell Sorter (Fig.S.4). Post-sort purity was higher than 90%.

RNA sequencing
15.000 to 200.000 CD14 hi monocytes from 11 healthy controls, 13 acute and 20 convalescent patients were isolated by FACS directly in TRIzol reagent (ThermoFisher scientific). After chloroform extraction, RNA isolation was performed using a RNeasy kit (Qiagen) and sample quality was tested on a Fragment Analyzer (Agilent). Indexed cDNA libraries were obtained using the Ovation Solo RNA-Seq System (Tecan) following manufacturer's recommendations.
The multiplexed libraries were loaded on a NovaSeq 6000 (Illumina) using a S2 flow cell and then mapped to the reference genome GRCh38 by using STAR_2.5.3a software with default parameters. We then sorted the reads from the alignment according to chromosome positions and indexed the resulting BAM-files. Read counts in the alignment BAM-files that overlap with the gene features were obtained using HTSeq-0.9.1 with "--nonunique all" option (if the read pair aligns to more than one location in the reference genome, it is counted in all features to which it was assigned and scored multiple times). Genes with no raw read count greater or equal to 20 in at least 1 sample were filtered out with an R script, raw read counts were normalized and a differential expression analysis was performed with DESeq2 by applying an adjusted p-value <0.05 and an absolute log2-ratio larger than 1. Regions obtained by MACS2 were merged to create an atlas containing all obtained peaks for all the populations using bedtools (54) with a minimum overlapping of 1bp. Merged regions were subject to differential analysis using csaw workflow (55) .For downstream visualisation, a scaling factor was calculated using deepTools package (56) to normalise peak intensity to FRiP (fraction of reads in peaks) and generate bigwig files.
For Gene ontology analysis, we introduced BED files from differential ATAC-seq peaks to GREAT tool with default parameters (23). For motif analysis, Ciiider algorithm (http://ciiider.com/) was used to perform motif enrichment in the differentially accessible regions.
We used BETA package (22) with default parameters to integrate ATAC-seq (differentially accessible regions) and RNA-seq (transcriptome) data and evaluate the regulatory potential of chromatin accessibility to promote/repress genes expression.

Data availability
RNA-seq and ATAC-seq data that support the findings reported in this study have been deposited in the GEO Repository with the accession code (pending).

Unsupervised analysis
First, the centered and scaled cellular phenotypic data corresponding to patients in the different clinical categories was visualized in two dimensions using t-distributed stochastic neighbor embedding (t-SNE). The observed clustering was stable across t-SNE technical replicates. To assess which cellular phenotypes, i.e. features, best explained the differences across the disease stages, we computed P-values from a Kruskal-Wallis test.