key: cord-0292571-q7fqq33u
authors: Huoman, J.; Sayyab, S.; Apostolou, E.; Karlsson, L.; Porcile, L.; Rizwan, M.; Sharma, S.; Das, J.; Rosen, A.; Lerm, M.
title: Mild SARS-CoV-2 infection modifies DNA methylation of peripheral blood mononuclear cells from COVID-19 convalescents
date: 2021-07-06
journal: nan
DOI: 10.1101/2021.07.05.21260014
sha: fc4a00cbff5f32c37ba6103a1a03ac9aaa84dce7
doc_id: 292571
cord_uid: q7fqq33u

Background: Coronaviruses such as SARS-CoV-2 may circumvent host defence mechanisms by hijacking host proteins, possibly by altering DNA methylation patterns in host cells. While most epigenetic studies have been performed in severely ill COVID-19 patients, studies on individuals who have recovered from mild-to-moderate disease remain scarce. The aim of this study was to assess epigenome-wide DNA methylation patterns in COVID-19 convalescents compared to uninfected controls from before and after the pandemic outbreak began. Methods: DNA was extracted from peripheral blood mononuclear cells originating from uninfected controls before (Pre20, n=5) and after (Con, n=18) 2020, COVID-19 convalescents (CC19, n=14) and symptom-free individuals with a SARS-CoV-2-specific T cell response (SFT, n=6), as well as from Pre20 (n=4) samples stimulated in vitro with SARS-CoV-2. Subsequently, epigenome-wide DNA methylation analyses were performed using the Illumina MethylationEPIC 850K array, and statistical and bioinformatic analyses comprised differential DNA methylation, pathway over-representation and module identification network analyses. Results: DNA methylation patterns of COVID-19 convalescents were altered as compared to uninfected controls, with similar results observed in in vitro stimulations of PBMC with SARS-CoV-2. Differentially methylated genes from the in vivo comparison constituted the foundation for the identification of a possibly SARS-CoV-2-induced module, containing 66 genes of which six could also be identified in corresponding analyses of the in vitro data (TP53, INS, HSPA4, SP1, ESR1 and FAS). Pathway over-representation analyses revealed involvement of Wnt, cadherin and apoptosis signalling pathways amongst others. Furthermore, numerous interactions were found between the obtained differentially methylated genes from both settings and the network analyses when overlaying the data unto the SARS-CoV-2 interactome. Conclusions: Epigenome-wide DNA methylation patterns of individuals that have recovered from mild-to-moderate COVID-19 are different from those of non-infected controls. The observed alterations during both in vivo and in vitro exposure to SARS-CoV-2 showed involvement in interactions and pathways that are highly relevant to COVID-19. The present study provides indications that DNA methylation is one of several epigenetic mechanisms that is altered upon SARS-CoV-2 infection. Further studies on the mechanistic underpinnings should determine whether the observed effects are reflecting host-protective antiviral defence or targeted viral hijacking to evade host defence.

1 8 1 closely together ( Figure 2C , Figure S1 ). The observed methylome differences prompted us to 1 8 2 identify differentially methylated CpGs (DMCs), which we defined as CpG sites with a Table S2a ). This identified DMC signature could furthermore accurately distinguish the . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint 

. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi. org/10.1101 org/10. /2021 As similar pathways were revealed in the findings from the clinical study and the SARS-CoV- comparison. These analyses found a module consisting of six genes (TP53, INS, HSPA4, 2 9 0 SP1, ESR1 and FAS), which were among the previously identified module genes from the in 2 9 1 vivo setting and also were identical to those that had been identified by more than two 2 9 2 module identification methods ( Figure 4A ). Furthermore, explorations of the overlap between 2 9 3 identified genes in the differential DNAm analyses and network module analyses to the 2 9 4 genes from the SARS-CoV-2 interactome identified numerous interactions in the in vivo 2 9 5 (n=11/54), in vitro (n=100/542) and network module setting (n=33/66) ( Figure S3 ). . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint The epigenetic events triggered during a mild COVID-19 infection are largely unknown, interferons were shown to be transcriptionally inhibited by hypermethylation, in severely ill 3 1 0 COVID-19 patients compared to controls, while genes originating from inflammatory 3 1 1 responses were granted transcriptional accessibility through general hypomethylation (17).

Other studies reported on DNAm patterns in whole blood of COVID-19 sufferers, comparing several antiviral immunity-related pathways. Whether the changes we found are reflecting an 3 1 6 antiviral defence mechanism or reflect viral manipulation of the host epigenome warrants 3 1 7 further studies. The main finding of our study was that a number of genes in the networks deriving from both 3 1 9 the DNA methylomes of mildly ill COVID-19 subjects and the in vitro stimulated PBMCs 3 2 0 methylomes were shared and consistently found by several module identification methods. This may indicate the importance of these genes as hubs for protein-protein interactions in 3 2 2 the course of SARS-CoV-2 infection and recovery. One of these genes was tumor protein 53 3 2 3 (TP53), an evolutionarily conserved protein that is one of the most well-studied hub genes in 3 2 4 . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint 16 cell signalling due to its central role in cancer (21) and that interacts with a variety of viral 3 2 5 proteins from different classes (https://thebiogrid.org)(22). The ability of mutual inhibition and 3 2 6 downregulation has been shown for TP53 and one of the previously identified SARS 3 2 7 coronaviruses -SARS-CoV (23). Furthermore, TP53 has in several other studies been 3 2 8 identified as a hub gene, in whole blood from COVID-19 patients (24), and interacting with 3 2 9 ACE2 in SARS-CoV-2-infected human induced pluripotent stem cell-derived cardiomyocytes 3 3 0 (25). In line with findings from our study, transcriptomic analyses of PBMCs from a small Furthermore, apoptosis of T cells in PBMCs induced by FAS was reported to be increased in Interestingly, reports on differentially expressed genes overlapping between acute respiratory (32), suggesting that this gene may be of particular importance in the mounting of a CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint shock protein, a member of the HSP70 family, also binds proteins of the Human Herpes Virus 4 and HIV. Although our study does not provide any evidence for a protective role, 3 5 3 HSP70 family members have been discussed as antiviral defence components (33, 34). In correlations of the inactivating delta F508 polymorphism, which is protective against chloride 3 5 7 ion secreting diarrhoeas, with prevalence and mortality in . This is particularly

interesting as SARS-CoV-2-induced diarrhoea has been suggested to involve Ca2+ activated 3 5 9 chloride channels (37). Similarly, it has been hypothesised that transport of chloride ions over 3 6 0 CFTR may be implicated in Covid-19-induced lung oedema (38). Furthermore, the 3 6 1 muscarinic acetylcholine receptor 1 and 3 signalling pathway was present in over- was not investigated in our study, this could suggest that these pathways found may be 3 6 7 implicated in the development of for instance post-acute COVID-19 syndrome, as the effects 3 6 8 we observe may have persisted for months after the initial exposure to the virus. Altogether, 3 6 9 the network centrality of the hub genes that we derived from the in vivo and in vitro data 3 7 0 suggests that they may be of particular importance in the interaction with epigenetically Although an obvious limitation of the study is the lack of validation of the DNAm findings on a 3 7 4 transcriptional level, it serves as a pilot study that generates hypotheses for further studies 3 7 5 within the field. Hence, whether the observed DNAm patterns are indeed associated or even . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint 18 possibly also in closer proximity to the time of infection with SARS-CoV-2, it will be possible 3 7 9

to study the role of DNAm alterations in anti-viral defence and in viral manipulation of the 3 8 0 same.

An advantage of the investigation of epigenetic modifications in in mild to moderately ill 3 8 2 patients, is that we may be able to discern DNAm differences that otherwise would have 3 8 3 been masked due to an overriding inflammatory response. These subtle changes may not 3 8 4 only be relevant to how a less severe immune response is mounted towards SARS-CoV-2, In conclusion, we found epigenome-wide differences in DNAm patterns of individuals that 3 9 2 had recovered from a mild-to-moderate disease course of COVID-19 compared to non- analyses. The study suggests that DNAm is one of several epigenetic mechanisms that are 3 9 6 altered upon SARS-CoV-2 infection. However, whether the effects are reflecting targeted 3 9 7 viral hijacking to evade host defence or host-protective antiviral defence mechanisms 3 9 8 remains to be determined. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. In this study, participants were enrolled between May 29 th and July 10 th 2020 during the first CoV-2 infection and/or other infectious disease symptoms, as well as being younger than 18 4 1 7 years. The study participants voluntarily entered the study in a consecutive manner. The Additionally, blood samples from anonymous healthy blood donors from the blood bank at 4 2 2

Linköping University Hospital before 2020 were included as a separate group in the analyses participants presented with either mild or asymptomatic initial infection, and none was . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. well. From the included individuals, the following information was retrieved using health 4 3 0 questionnaires: self-reported COVID-19 symptoms (if applicable, one or several of the 4 3 1 following: fever, headache, shortness of breath, loss of smell/taste, cough, fatigue, muscle 4 3 2 pain, nausea, sinusitis/congestion), date of self-reported symptoms, weeks between 4 3 3 symptoms and sampling, age, sex, smoking, weight, height, comorbidities as well as 4 3 4 medications. The blood and saliva from the study participants was processed in a Biosafety informed consent, and the present study was approved by the Regional Ethics Committee for Sweden) further on termed as complete culture medium, prior to further processing. Up to 10 4 5 4 . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint 21 ml of whole blood was used for plasma separation by centrifugation (2000g for 15min, 4°C) 4 5 5 and aliquots were stored at -80°C till further analysis. contains 166 peptides consisting of 15-mers, overlapping with 11 amino acids, covering the 4 6 2 S1 domain of the spike S1 protein (amino acid 13-685). The peptides were combined into 4 6 3 one pool. IFN-γ ELISpot Plus kit was purchased from Mabtech (3420-4HST-10, Sweden).

Briefly, the pre-coated wells were plated with unfractionated PBMCs at counts of 300 000 4 6 5 cells/well, and the cells were cultured with peptides for the S protein of SARS-CoV-2 at a 4 6 6 final concentration of 2 µg/ml (diluted in complete culture medium) for 20 to 22 hrs in a 37°C, 4 6 7 5% CO 2 incubator. Cells cultured with medium alone were used as negative controls. when the stimulation index was >2, and the increment value was >10. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. Prior to saliva collection, participants were required to rinse their mouth with water and 4 8 3 confirmed they did not show documented oral disease or injury, that they had fasted, 4 8 4 refrained from smoking, chewing a gum, taking oral medication, tooth brushing for a 4 8 5 minimum of 1 hour before sampling and that no dental work had been performed within 24 4 8 6 hours prior to sample collection. Donors were asked to provide a 5 ml sample of saliva in a 4 8 7 50 ml sterile conical tube by passive drooling. processing to preserve sample integrity. Samples were centrifuged (2500g for 20 minutes at MagPlex-C microspheres (Luminex Corp., Austin, TX, USA) were used for the coupling of . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. coupled beads were stored at 4°C in the dark until further use. bead mixture suspended in PBS-T + 1% BSA (~50 beads/µl) was then added to each well.

The plate was incubated in the dark at 600 rpm for 1h at RT. The wells were then washed 100 events for each bead number was set to read and the median value was obtained for the 5 2 2 analysis of the data. All sample analyses were repeated three times. A naked, non-antigen-5 2 3 coupled bead was included as a blank along with PBS-T + 1% BSA as a negative control. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. FBS and 100 μ g/ml gentamicin) was added, and the plate was incubated for 8 hours at 37°C in presence of 5% CO 2 . After incubation, the supernatant was discarded, and the cells were fixed for 2 hours with 4% formaldehyde. Next, Triton-X (1:500 in phosphate buffered saline, . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint 25 (Scions, Code: J2 at 1:100 dilution) for 1.5 h followed by detection using horseradish lysate was titrated to be 5×10 6 per ml.

Cells were monitored in the IncuCyte S3 live cell analysis system (Sartorious, Göttingen,

Germany) to allow quantification of cell death in SARS-CoV-2 infected wells versus controls.

After 48h incubation the cell culture media was collected from each well and centrifugated at For the performance of epigenome-wide DNA methylation analyses, DNA was extracted from . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. 850K array. 200 ng of DNA from each sample was analysed. Initial descriptive analyses of demographic variables were performed on the available 5 9 4 information about age, gender, smoking and BMI (kg/m 2 ). Continuous variables were 5 9 5 compared using an unpaired two-tailed t-test and categorical variables were examined using 5 9 6

the Pearson χ 2 test or Fisher's exact test (if the number of observations was smaller than 5 9 7 five), see Table S1 . The resulting raw IDAT-files from the MethylationEPIC array analyses were processed in R The resulting raw IDAT-files containing the raw DNA methylation profiles for each cell type were pre-processed in several steps. The following filters were applied: i) removal of probes probes, iv) removal of all probes in X and Y chromosomes. We removed the sex 6 0 8 chromosomes from our data set, as female X-inactivation skews the distribution of beta values ( Figure S4 ). Of the initial 865 918 probes, 841 524 probes remained upon filtering. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. assessed before and after the normalisation ( Figure S5 ). Thereafter, we performed singular 6 1 4

value decomposition (SVD) analyses using the ChAMP package (45) identify underlying components of variation within the filtered and normalised data set ( Figure   6 1 6 S6). Significant components consisted of slide, batch and sample groups that contributed to 6 1 7 variation within the data set. Corrections were performed for the identified components using 6 1 8

ComBat from the SVA package (46) (version 3.38.0). As PBMCs consist of multiple 6 1 9 nucleated cell types in peripheral blood, we utilised the Houseman method to infer cell type proportions between any of the individuals or between sample groups (Table S7) , motivating 6 2 2 our choice of not correcting for these cell type proportions. Differential DNA methylation analysis in vivo 6 2 5

As we were interested in studying CpGs that were differentially methylated between CC19s 6 2 6

and non-infected controls from both before and after the start of the COVID-19 pandemic, we 6 2 7

performed differential DNA methylation analyses, using the limma package (version 3.46.0).

A linear model was fitted to the filtered, normalised and SVD-corrected DNA methylation 6 2 9

data. Identified sources of variation that were still present upon SVD correction provided the 6 3 0 basis for the inclusion of these variables as co-variates in the models, in this case gender 6 3 1 and BMI ( Figure S6 ). For each investigated probe, moderated t-statistics, log2 Fold Change between the CC19s vs. non-infected controls (Cons + Pre20). Differentially methylated CpGs 6 3 5

(DMCs) were defined as CpG sites having a nominal p-value of less than 0.01 along with an 6 3 6

MMD of > 0.2. As a means to ascertain the quality of the identified DMCs, genomic inflation 6 3 7

. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint 28 and pertaining bias were estimated using the BACON package (version 1.18.0). As the 6 3 8 estimated genomic inflation for the comparison was close to 1 (genomic inflation: 1.20, bias: 6 3 9 0.01, Figure S7 ), this suggested that no major genomic inflation was present in the 6 4 0 comparisons, and no correction for this was deemed necessary. The distribution of the The resulting raw IDAT-files containing the raw DNA methylation profiles for each cell type were pre-processed in several steps. The following filters were applied: i) removal of probes 6 5 1 with detection p-values above 0.01, ii) removal of non-CpG probes, iii) removal of multi-hit 6 5 2 probes, iv) removal of all probes in X and Y chromosomes. In this dataset, we did not have 6 5 3 any information on demographic variables, as the samples derived from anonymous donors. However, we still removed the sex chromosomes from our data set, as female X-inactivation 6 5 5 skews the distribution of beta values. Of the initial 861 728 probes, 837 694 probes remained 6 5 6 upon filtering. After filtering, quality control was performed, and normalisation of the data was (Table S7) , motivating our choice of not correcting for these cell type proportions. The β - quality of the data was assessed before and after the normalisation ( Figure S8 ). SVA 6 6 3 . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint

Epigenetic Landscape during Coronavirus Infection

Epigenetic mechanisms influencing COVID-191

Epigenetic mechanisms regulating COVID

package (version 3.40) was applied to correct the batch effect. Cell deconvolution was 6 6 4 performed using FlowSorted.Blood.EPIC package (version 1.11). 6 6 5 6 6 6 Differential DNA methylation in vitro 6 6 7 To evaluate the difference between the MOCK and INFECTION, the fold change was 6 6 8 calculated using the cut-off obtained from the density plot (M-value >|2|; Figure S9 ) for each 6 6 9CpG site. Only those CpGs with higher values than the cut-off, were selected for further 6 7 0 analysis. Venn analysis was performed among the samples using the ggVennDiagram To make biological sense of the putatively SARS-CoV-2-induced DNA methylation 6 7 4 differences, we performed PANTHER pathway over-representation test analyses using the 6 7 5 PANTHER database (version 16.0). The Fisher's exact test was used for generation of 6 7 6 nominal p-values (significance level set to p-value of < 0.05), in case false discovery rate 6 7 7 correction was too stringent. The significantly enriched pathways were displayed in dot plots 6 7 8 generated in R using ggplot2 package (version 3.3.3). DMGs generated in the in vivo setting. An input object was constructed using the pre-2020 6 8 2 (Pre20, n=5) and post-2020 (Con, n=18) non-exposed controls and COVID-19 convalescents 6 8 3 (CC19, n=14), as a two-column data frame containing gene annotation and P-value of the 6 8 4 significant DMGs (n=54). The graph clustering algorithm MCODE (48) was used to identify 6 8 5 molecular complexes and create a large disease module, which was then fitted to a protein- . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint 30 in the network. Centrality measurements of degree, betweenness and closeness were used 6 8 9to expose the most central nodes in the network. Finally, a functional enrichment of the 6 9 0 genes present within the module was carried out using StringDB (49). In addition, the 6 9 1 inference of modules was performed with two other methods from the MODifieR package 6 9 2 (DIAMOnD and WGCNA)(50) to study whether it was possible to condense the module 6 9 3 genes to fewer genes of particular interest within the network, for both the in vivo and the in 6 9 4 vitro setting. Overlap to SARS-CoV-2 interactome 6 9 7A publicly available protein-protein interaction (PPI) network of SARS-COV-2 and human . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint 32 Pre20 -Control (uninfected, pre-pandemic 2020) T-cell responses associated with natural HIV control. PLoS pathogens. 2020;16(8). macrophages from TB-exposed individuals predicts exposure to mycobacteria. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 6, 2021. Availability of data and materials 9 0 5The datasets used and/or analysed in the presented work will be available upon publication 9 0 6 due to a pending patent, for reference: GeneExpression Omnibus (GEO-ID: GSE178962). Competing interests 9 0 8None of the authors declare any competing interests. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)The copyright holder for this preprint this version posted July 6, 2021. ; https://doi.org/10.1101/2021.07.05.21260014 doi: medRxiv preprint