DNA methylation signatures of severe RSV infection in infants: evidence from non-invasive saliva samples

Abstract

Background

Respiratory syncytial virus (RSV) poses significant morbidity and mortality risks in childhood, particularly for previously healthy infants admitted to hospitals lacking predisposing risk factors for severe disease. This study aimed to investigate the role of the host epigenome in RSV infection severity using non-invasive buccal swabs from sixteen hospitalized infants admitted to the hospital for RSV infection. Eight patients had severe symptoms, and eight had mild to moderate symptoms. For DNA methylation analyses, the Illumina EPIC BeadChip was used with DNA isolated from saliva samples. To evaluate the basal DNA methylation level of the identified biomarkers a cohort of healthy control children was used. Furthermore, DNA methylation levels of candidate genes were confirmed by pyrosequencing in both the discovery and validation cohorts of patients with mild to moderate symptoms.

Results

A panel of differentially methylated positions (DMPs) distinguishing severe from mild to moderate symptoms in infants was identified. DMPs were determined using a threshold of an adjusted P-value (false discovery rate, FDR) < 0.01 and an absolute difference in DNA methylation (delta beta) > 0.10. Differentially methylated regions (DMRs) were identified in the ZBTB38 (implicated in asthma and pulmonary disease) and the TRIM6-TRM34 gene region (associated with viral infections). The differential DNA methylation of these genes was validated in an independent replication cohort. A weighted correlation network analysis emphasized the pivotal role of a module with RAB11FIP5 as the hub gene, known for its critical function in regulating viral infections.

Conclusions

Oral mucosa methylation may play a role in determining the severity of RSV disease in infants.

Background

Respiratory syncytial virus (RSV), the primary cause of acute lower respiratory infections (ALRIs) in young children across the globe, infects almost all children by the age of two and is responsible for a considerable number of infant hospitalizations [1,2,3,4]. Although most RSV hospitalizations occur in previously healthy children born at term, current prevention strategies still primarily focus on specific at-risk children [5]. Evidence suggests that the severity of RSV bronchiolitis occurring in infancy could be associated with immunological mechanisms [6]. In fact, RSV infections occur more frequently in the first years of life when the immune system is immature.

Understanding the underlying molecular mechanisms of RSV infection and the resulting immune response is important for developing more effective treatments and preventive measures.

While several studies have been carried out on host genome susceptibility to RSV infection [7, 8] and gene expression [9,10,11], only a few of them have been performed in the field of epigenomics [12]. Epigenetic regulation serves as a dynamic interface between the genome and the environment. Thus, in the context of viral or bacterial infections, the regulation of host defense cells by epigenetic mechanisms is directly linked to the development of the disease [13]. Evidence suggests that epigenetic processes play a critical role in modulating the interaction between the host and the causal pathogen [14].

Recently, we have demonstrated the role of blood DNA methylation (DNAm) in the development of respiratory sequelae after RSV infection [15]. One year later, in 2023, Zhaozhong Zhu et al. [16] conducted an epigenome wide association study (EWAS) on blood samples to study the DNAm patterns in infants hospitalized with bronchiolitis, 423 of whom were infected by RSV. While blood is widely used in epigenomic biomedical research and is considered a valuable tissue for biomarker discovery, it is important to recognize the advantages of alternative sample types. Saliva, in particular, has gained attention due to its non-invasive collection process, particularly for research studies involving children. In recent years, saliva and buccal swabs samples have emerged as crucial samples in different areas of medical research, particularly for neurodegenerative diseases [17, 18], due to its ability to reflect brain-related biological processes [19].

Saliva is an attractive, accessible source of cells from which a high quantity and quality of DNA/RNA can be obtained. It is known that both white blood cells originating from mesoderm and epithelial cells from ectoderm are present in the mouth [20]. Despite these data, very few studies exist in the literature on the identification of molecular host factors in saliva in the context of infections [21, 22]. Recent studies have provided evidence that RSV, as other respiratory viruses, is frequently detected in saliva and sputum samples [23]. However, no research have explored the potential of using buccal swab to investigate DNAm signatures associated with RSV disease or used epigenomics in this context.

Methods

Aim of the study

We aimed to investigate the role of host DNA methylation (DNAm) in the severity of RSV infection using buccal swab samples collected from patients with mild/moderate and severe symptoms during the acute phase of infection. To the best of our knowledge, this is the first time that DNAm has been explored to study RSV severity in buccal swab samples.

Study design

Sixteen patients with RSV infection were chosen from a cohort that was prospectively enrolled in an observational study conducted in Spain, facilitated by a national hospital-based research network for pediatric respiratory research called GENDRES (Genetics, Vitamin D and Respiratory Infections Research Network—www.gendres.org). RSV infection was confirmed by polymerase chain reaction (PCR) upon presentation to the hospital. A saliva ORAGENE buccal swab sample was collected from each participant during the acute phase of RSV infection, typically within seven days of symptom onset.

Table 1 presents the clinical history and characteristics of the study subjects. The clinical severity of the disease was determined based on the need for intensive care, invasive respiratory support, and the ReSViNET score [24]. The ReSViNET score, used as a parameter of severity classification, is a clinical scoring system used to assess patients’ severity of acute respiratory infections (ARIs). Seven parameters were considered: feeding intolerance, medical intervention, respiratory difficulty, respiratory frequency, apnea, general condition, and fever. Each parameter is assigned a value ranging from 0 to 3, resulting in a total score of 20 points.

Table 1 Clinical characteristics of the study participants

The RSV cases were obtained from a cohort of more than 3000 patients, over half of whom were infected with RSV. To select severe RSV cases, children with well-defined severity, such as acute bronchiolitis (lower respiratory tract infection [LRTI]), who have the highest ReSViNET score (> 15), who require mechanical ventilation, and who spend the most time in the pediatric intensive care unit (PICU) were chosen. The moderate/mild group was formed by selecting children who only required a few days of hospitalization (less than 5), did not require any kind of respiratory support, and had a ReSViNET score less than 8. This group of children was diagnosed with both LRTI and upper respiratory tract infections (URTIs).

Sample processing

DNA was isolated from saliva samples following the prepIT-L2L (DNAgenotek) protocol. The buccal swabs were incubated at 50 °C in a water incubator for at least one hour. The saliva sample was then transferred to a microcentrifuge tube and mixed with a 1/25 volume of PT-L2P. The mixture was vortexed and incubated on ice. After centrifugation, the supernatant was removed, and ethanol was added to allow DNA precipitation. Subsequently, an ethanol wash was performed, and the DNA pellet was then dissolved in TE solution. After isolation, the DNA sample concentration and quality were checked using common Nanodrop measurements and Qubit. Because DNA in buccal swabs tends to be more degraded than blood DNA, further DNA degradation control was performed to ensure the optimal DNA concentration for the subsequent steps. The whole Illumina methylation EPIC BeadChip experiment was carried out at the Genómica e Investigación Oncológica (GENYO) Center (Granada; Spain). The epigenome data were delivered to our group at Santiago de Compostela for subsequent bioinformatic and statistical analysis.

An appropriate volume of concentrated DNA was subjected to bisulphite conversion, a gold standard technique that allows the deamination of non-methylated cytosines into uracil and prevents methylation of the cytosine. After this treatment, the subsequent amplification recognized uracil as thymine and methylated cytosines as cytosines. The level of DNAm was then calculated by quantitative genotyping of C/T SNPs. After hybridization to the Illumina methylation EPIC BeadChip platform, the fluorescence intensity resulting from the incorporation of a nucleotide in the probe is translated into a level of DNAm for each CpG site. DNAm can be defined as beta- or M-values. Beta values range between 0 and 1, where 0 represents non-methylation, and 1 represents 100% methylation. Beta values are used for graphics due to their easy interpretability. M-values, on the other hand, are logit-transformed beta values; homoscedasticity is generally assumed for these values, and they are preferentially used for linear statistical analysis.

Statistical analysis

Methylation data analysis

Before calculating the beta- and M-values, the raw IDAT data were subjected to different filtering processes and normalization steps, as previously described [15], to obtain high-quality data suitable for downstream analysis.

In the analysis of DNAm data derived from ‘bulk’ tissues such as blood or saliva, it is crucial to estimate the composition of different cell types within the tissue. Since the measured DNAm in ‘bulk’ tissues represents an average of the methylation levels from multiple cell types, there is a risk that the overall DNAm level is predominantly influenced by specific cell types. Moreover, when employing saliva as a biomarker identification tool, it is important to consider its extreme cellular heterogeneity. It has been demonstrated that buccal swabs used to collect saliva samples contain a greater proportion of epithelial cells than saliva, and this proportion is greater in children than in adults. Therefore, estimating the cell composition is an essential step in this type of analysis. Here, we used the Epidish R package to estimate the cell composition of the buccal swab samples. The hepidish function of the Epidish R package, allowing the use of two different DNAm references, was applied to the beta values to specifically estimate epithelial cells, fibroblasts, and total immune cells (neutrophils, eosinophils, monocytes, CD4+ and CD8+ T cells, B cells, and natural killer [NK] cells) [25, 26].

Identification of RSV severity-associated differential methylation

The Limma package [27] was used to identify differentially methylated positions (DMPs), employing a model adjusted for age, sex, and cell estimation. To account for age, we included it as a covariate in the linear model formula. This allowed us to estimate the effect of other variables while controlling for the effect of age. Specifically, we created a design matrix including age as a covariate and fitted a linear model using the lmFit() function. A previous episode of bronchiolitis, which occurred in five out of eight mild/moderate patients, was considered in the statistical analysis to evaluate the effect of this clinical factor. The threshold for identifying DMPs was set at an adjusted P-value (false discovery rate or FDR) of < 0.01 and an absolute difference in DNAm between groups (delta beta) of > 0.10. Receiver operating characteristic (ROC) curve analyses were conducted to evaluate the diagnostic efficacy of the most significant candidate DMPs, and the determined area under the curve (AUC) was used to assess their discriminatory potential.

To identify differentially methylated regions (DMRs), the DMRcate package [28] was utilized with default parameters. The package employs Gaussian kernel smoothing to detect patterns of differential methylation independent of genomic annotation. In this study, a bandwidth (λ) of 1000 base pairs and a scaling factor (C) of 2 were used. Despite working with a sparse set of sites, the aim was to identify regions with multiple probe clusters (minimum of 3 CpGs per region) exhibiting substantial effect sizes and a minimum smoothed FDR < 0.01.

To establish a baseline sample of healthy controls (HC) for comparison, DNAm data from a cohort of healthy children were retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), under accession number GSE252169 [29]. This dataset was generated using the Illumina HumanMethylationEPIC BeadChip (GPL21145 platform). For the analysis, we only retained methylated profiles of buccal DNA samples from healthy offspring of no-smokers pregnant women, collected at three different time points: birth (n = 27), 12 months (n = 14), and 60 months (five years of age; n = 10).

The data were loaded into R using the GEOquery package. The raw IDAT files were preprocessed following the established protocol described previously. Beta and M values, were filtered to remove low-quality probes and normalized before proceeding with downstream analysis.

Pathway analysis

For gene set enrichment analysis, the methylGSA package [30] was applied, which uses the methylglm function.

This function extends GOglm by incorporating gene length as a covariate to adjust length bias in DNAm based on the number of CpGs. Pathway enrichment analysis was customized in two ways: considering all CpGs regardless of their gene group and focusing only on CpGs within the promoters based on the annotation in the “IlluminaHumanMethylationEPICanno.ilm10b4.hg19regions” (TSS200, TSS1500, 5’UTR, 1stExon).

Weighted gene correlation network analysis (WGCNA)

The WGCNA package [31] was used to construct a signed weighted correlation network, using the top 20% of positions with the highest median variance (approximately 120,000). To identify the gene module with the strongest association with RSV severity, a series of steps were performed. First, the soft threshold power was selected based on standard scale-free networks, and a power function was used to calculate all differential CpGs. The adjacency matrix was then transformed into a topological overlap matrix (TOM), and the corresponding dissimilarity matrix (1-TOM) was calculated. The dynamic tree cut method was employed to hierarchically cluster CpGs and identify the module, with a minimum module size of 30, medium sensitivity for cluster splitting, and a dendrogram cut height threshold of 0.25 for module merging as module detection parameters. Unassigned CpG sites were clustered in the “grey” module, which was not considered for further analyses. Gene significance (GS) was computed to detect statistically significant associations between modules and phenotype, and module membership (MM) was also calculated as a measure of intramodular connectivity by correlating (using Pearson correlation) the methylation profile with the eigengene of a given module. The correlation between GS and MM was investigated, and the average absolute gene significance for all CpGs within a module was determined to identify the most important modules. The modules were named according to their most significant hub genes. Furthermore, the hub genes of the RSV severity-associated modules were extracted, and an overrepresentation analysis of the most important module related to the trait was performed using the ClusterProfiler R package [32] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database as a reference to investigate its biological significance.

Pyrosequencing methylation analysis

The DNAm pattern of the identified DMPs was assessed by bisulfite-pyrosequencing in the discovery cohort, along with an independent validation cohort consisting of 14 patients (7 severe, 7 mild/moderate). Genomic DNA was bisulfite-converted using the EZ DNA Methylation Kit (Zymo Research) following the manufacturer’s protocol. Bisulfite-treated DNA was amplified using the PyroMark PCR Kit (Qiagen) in a 25 µL reaction volume. Forward and reverse PCR primers were designed to target the same DMPs identified in the microarray analysis. Following PCR amplification, pyrosequencing assays were performed in duplicate (technical replicates) in sequential runs on the PyroMark Q24 system (Qiagen) using PyroMark Gold Q24 Reagents (Qiagen). The sequencing primers (Primer S) are detailed in Table S1. Pyrosequencing was conducted by the Biomedical Research Institute of Murcia (IMIB).

The R statistical software (v. 4.2.2) was used to perform all the statistical analyses.

Results

Clinical characteristics of patients

Table 1 describes the clinical characteristics of the patients. Most of these characteristics exhibited similarities between the two severity groups. Clinical severity factors, such as the ReSViNET score, PICU admission, and the requirement for mechanical ventilation, were found to be significantly different between the two groups, as expected. Another statistically significant factor was age, with infants in the severe group being younger than those in the mild/moderate group. To investigate the association between the observed severity in the subjects and their “biological age” (i.e. the actual age of cells rather than the chronological one), we employed age predictors. Specifically, we utilized the Pediatric-Buccal-Epigenetic (PedBE) clock, an epigenetic age estimator designed for pediatric buccal swab samples [33]. Through this analysis, we found a strong correlation (rho, ρ = 0.96; P-value < 0.0001) between DNAm age and chronological age (see the scatterplot in Supplementary Fig. 1A). Consequently, we did not observe statistically significant differences in predicted biological age between the mild/moderate and severe groups, as depicted in Supplementary Fig. 1B. To mitigate confounding effects on the results, age was considered a covariate in all analyses.

Identification of differentially methylated positions and regions within the ZBTB38 and TRIM6-TRIM34 genes

The estimation of cell composition did not reveal any significant differences between the two cohorts. Therefore, to minimize potential confounding factors, age was the only significant variable adjusted in the linear model. Among the 803,527 total DNA probes analyzed, we identified 461 DMPs (assuming a nominal P-value < 0.01 and an absolute difference in methylation level [Delta Beta] > 0.10) between infants with severe infection and those with only mild/moderate symptoms (Table S2). Among these positions, the majority (n = 339, 73.5%) exhibited hypomethylation in the severe group compared to the mild/moderate group, while the remaining (n = 122, 26.5%) showed hypermethylation. These DMPs enabled the differentiation of patient methylation profiles into two distinct groups according to principal component analysis (PCA) and a heatmap (Fig. 1A, B). The main difference in the PCA results was from the first principal component (PC1), which accounted for 68.7% of the variation.

Fig. 1

Differentially methylated positions distinguish infants with severe infection from those with mild to moderate symptoms. A PCA and B Heatmap of the significant DMPs (nominal P-value < 0.01 and delta beta > 0.10) between the severe and mild/moderate groups

Notably, we also identified clusters of DMPs within promoter regions belonging to the same gene. Specifically, we found seven positions in the 5’UTR of the ZBTB38 gene (Table 2 and Fig. 2A) and three positions within the TSS200 region of the TRIM6-TRIM34 gene (Table 2 and Fig. 2B).

Table 2 Cluster of differentially methylated positions (DMPs) within the TRIM6-TRIM34 and ZBTB38 genes

Fig. 2

Hypomethylation of the ZBTB38 and TRIM genes in the severe group. Boxplot and ROC curve of the cluster of CpGs observed within A ZBTB38 and B TRIM6-TRIM34 genes

After adjusting the model for both age and bronchiolitis, we observed a decrease in the number of DMPs (n = 129; nominal P-value < 0.01). However, despite this reduction, the top DMPs remained consistent, particularly those located within the ZBTB38 (four positions) and TRIM6-TRIM34 (three positions) genes (Table 2). Additionally, removing patients who had experienced previous episodes of bronchiolitis from the statistical analysis resulted in a slightly different but even larger group of DMPs (n = 465). However, what is more relevant again is the remarkable consistency in the top DMPs and DMRs, pointing to ZBTB38 (with two positions) and TRIM6-TRIM34 (with five positions); (Table 2). This compelling consistency implies that previous episodes of bronchiolitis might not exert a widespread impact on DNAm across the entire mild/moderate cohort, reinforcing the robustness of the findings.

To assess the discriminatory potential of DMPs from these genes (ZBTB38 and TRIM6-TRIM34) as biomarkers of RSV severity, we performed receiver operating characteristic (ROC) curve analysis with area under the curve (AUC) calculations. The total ROC curve for the seven DMPs in the ZBTB38 gene yielded an AUC of 0.984 (0.941–1.000), while the ROC curve for the three DMPs within the TRIM6–TRIM34 region had an AUC of 0.906 (0.756–1.000). Furthermore, we utilized DMRcate to identify DMRs comprising at least three neighboring CpG sites associated with RSV severity in the same direction. This approach revealed eight regions with decreasing methylation and twelve regions with increasing methylation (Table 3). This analysis detected two highly statistically significant regions (Table 3). The most significant one (false discovery rate or FDR = 2.83 × 10^–24) overlaps the TRIM6-TRIM34 genes, and it comprises 11 highly correlated CpGs (Fig. 3A). The second most significant region (FDR of 1.93 × 10^–11) overlaps with the ZBTB38 gene; it contains 15 highly correlated CpGs (Fig. 3B).

Table 3 Differentially methylated regions (DMRs) associated with RSV severity

Fig. 3

Differentially methylated regions distinguish infants with severe infection from those with mild to moderate symptoms. Top panel: DMR plot for the two most significant DMRs within the A ZBTB38 and B TRIM6-TRIM34 genes. DNA methylation patterns in severe and mild/moderate samples differed for most of the CpGs in the DMRs. In the figure, some beta value dots of the CpGs that constitute the DMRs overlap in the two groups with similar values. Bottom panel: pairwise correlation between CpG sites in DMRs in A ZBTB38 and B TRIM6-TRIM34 genes

To verify the basal DNAm levels of the key genes identified in our analysis, we utilized the GEO dataset GSE252169. This dataset was selected for its availability of DNA buccal samples from healthy children at different time points. Specifically, we assessed the basal methylation profiles of the seven CpGs within ZBTB38 and three CpGs within TRIM6-TRIM34, identified through DMPs and DMRs analysis. The GEO subset cohort used for this study included a group of 27 children who were not exposed to maternal smoking in utero, with DNA samples collected at birth, 12 months, and 60 months. Normalized beta values for the 10 CpGs were extracted from this HC cohort. During preliminary analysis, two CpGs within ZBTB38 (cg25395158 and cg23967605) were removed due to data quality issues. The remaining five CpGs within ZBTB38 were analyzed for their methylation levels at the three time points in the HC cohort and these levels were then compared to the methylation profiles of these CpGs in children with RSV.

As illustrated in Fig. 4, the methylation levels of the five CpGs within ZBTB38 at birth in the HC cohort were similar to those observed in children with mild to moderate symptoms. However, as the children aged, methylation levels at these positions in HC cohort increased, demonstrating an opposite trend to that observed in the RSV cohort, where a decrease in DNAm was detected for children with severe infection when compared to those with only mild to moderate symptoms. The analysis of the three CpGs within TRIM6-TRIM34 is presented in Supplementary Fig. 2.

Fig. 4

DNA methylation levels of DMPs within ZBTB38 in RSV and HC cohorts. Boxplots showing DNA methylation levels at five DMPs within ZBTB38 for the RSV cohort compared to the HC cohort from GEO. The HC cohort includes samples collected from children at three-time points: birth, 12 months, and 60 months. The RSV cohort is represented by children with mild to moderate symptoms. Methylation levels at birth are similar to the methylation levels of RSV children with mild to moderate symptoms. For all the CpGs, the HC cohort shows an increase of methylation with age, in contrast to the RSV cohort, where hypomethylation is observed for children severely ill

Pathway analysis

We next investigated gene sets from Reactome, Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways by considering (i) all CpGs and (ii) only CpGs within promoter regions. Remarkably, for both analyses, significant enrichment was observed in three Reactome pathways, three GO biological processes, and one KEGG pathway (Fig. 5). Notably, all three enriched pathways included olfactory signaling pathways, such as transduction and receptor activity.

Fig. 5

Pathway analysis revealed enrichment of pathways associated with olfactory signaling and activation. Bar plot representing the results of the enrichment pathway analysis conducted on all CpGs and specifically on CpGs located within the promoter regions using the methylglm approach. The results were obtained using the KEGG, Reactome and GO databases. The size of each bar along the x-axis indicates the number of genes associated with each pathway. Furthermore, false discovery rate (FDR) P-values are indicated

Technical and biological validation of candidate methylation sites

Seven new patients with severe RSV infection and seven new children with mild to moderate symptoms were selected to evaluate the accuracy of DNAm results obtained from the Illumina EPIC BeadChip assay. Table 4 presents the clinical characteristics of these subjects.

Table 4 The clinical characteristics of mild/moderate and severe patients in the validation cohort

In this new validation cohort, the RESVINET score, need for respiratory support, and PICU admission remained significantly different between severe and mild/moderate cases, consistent with findings from the discovery cohort. To assess the reproducibility of the pyrosequencing assay, technical replicates of the ten DMPs were performed in both the discovery and validation cohorts. Each sample was analyzed in duplicate, and Pearson correlation analysis was conducted to evaluate the consistency between replicates (Fig. 6). The correlation between duplicates was extremely high in both cohorts, with correlation coefficients of 0.9977 and 0.9965 respectively and a highly significant P-value (< 2.16 × 10⁻¹⁶). These results confirm the robustness and reliability of the pyrosequencing assay in accurately quantifying DNAm levels across different cohorts.

Fig. 6

High reproducibility of pyrosequencing technical replicates in discovery and validation cohorts. Correlation plots showing the Pearson correlation between technical replicates for each cohort analyzed using pyrosequencing. Both the Discovery and Validation cohorts demonstrate near-perfect correlation, with Pearson correlation coefficients close to 1 and a highly significant P-value (< 2.2 × 10⁻¹⁶), confirming the high reproducibility of the analysis

The DNAm levels of the 10 DMPs quantified by pyrosequencing in the discovery cohort were found to be significantly differentially methylated between severe and mild/moderate cases, fully aligning with the results from the Infinium assay (Fig. 7A, B). This confirms that the results of the Infinium data were successfully (technically) validated by pyrosequencing. Both technical replicates in the discovery cohort showed highly consistent methylation levels, further supporting the reproducibility of the assay.

Fig. 7

Methylation levels of DMPs in discovery and validation cohorts. Boxplots showing the distribution of DNAm levels for the 10 DMPs analyzed using pyrosequencing in the Discovery and Validation cohorts and in each replicate. The analysis includes seven CpG sites in the ZBTB38 gene and three CpG sites in the TRIM6-TRIM34 gene, confirming the reproducibility of methylation patterns across independent cohorts

In the Validation cohort, six of the seven DMPs within the ZBTB38 gene remained significantly differentially methylated, while the cg23371833 site showed a similar trend of hypomethylation in the severe group, albeit without reaching statistical significance (Fig. 7A). The three DMPs within the TRIM6-TRIM34 locus did not reach statistical significance in the validation cohort (Fig. 7B). However, when both cohorts were analyzed together (Table 5), all 10 DMPs reached statistical significance, reinforcing their role as robust prognostic biomarkers of RSV severity. These findings confirm the significant hypomethylation pattern observed in children with severe phenotype when compared to those with only mild to moderate symptoms, as identified in the discovery stage and in the Infinium assay.

Table 5 Differential methylation analysis by bisulfite pyrosequencing for each cohort of the study

Co-methylated modules were constructed using the WGCNA R package with the top 20% of probes exhibiting the highest mean–variance (n = 114,806). After testing a set of candidate powers, a soft-thresholding power of seven was selected based on the criterion of scale-free topology (Fig. 8A). The analysis revealed 32 modules, each containing at least 30 co-methylated genes, labeled by a color name (Fig. 8B). Applying a cutoff height threshold of 0.25 for module merging resulted in the identification of 21 modules. The module-trait relationships were tested by calculating the correlation between the modules and the severity phenotype (Fig. 8C), resulting in two significantly associated modules: RAB11FIP5 (green) and 2RX3 (magenta) (Fig. 8D). However, after adjusting for multiple testing, only the RAB11FIP5 module survived with a FDR P-value < 0.011. We found a strong correlation between GS and MM for this module, indicating that the CpGs associated with RSV severity are also core elements of that module (Fig. 9A).

Fig. 8

Co-methylated modules associated with RSV severity. A Selection of the soft-thresholding power involved analyzing plots that depict the correlation between the soft-thresholding powers and two metrics: the scale-free fit index (left) and the mean connectivity (right). B Clustering dendrogram of probes, with dissimilarity based on topological overlap, together with assigned module colors. C Bar plot displaying the P-values from correlation tests between the module eigengenes and the RSV severity phenotype. The color of each module is defined by its hub gene. D Heatmap of Pearson correlation analysis of modules and the clinical traits associated with RSV severity. The rows represent the 21 module eigengenes, and the columns represent the phenotypes

A strong negative correlation (R = −0.77) between the RAB11FIP5 module and the trait was observed, indicating hypomethylation in the severely ill group compared to the mild/moderate group. This finding is consistent with the results obtained for DMPs and DMRs, where lower levels of methylation were observed in the severe group than in the mild/moderate group. The heatmap and sample eigengene plots strongly support these results, demonstrating an overall lower methylation value in the cohort of patients with severe symptoms (Fig. 9B, C). Within the RAB11FIP5 module, 153 CpGs of the previously identified 461 DMPs (33%) were included, of which 97 (21%) were considered hub CpGs (GS < −0.7 and MM > 0.8).

Fig. 9

Strong negative correlation between the RAB11FIP5 module and RSV severity. A Plot of correlations between gene significance and module membership for the RAB11FIP5 module. Color coding is equivalent to module names. B Heatmap displaying the beta-values of the CpGs belonging to the RAB11FIP5 module. Additionally, the eigengene values of the samples are also shown. C Boxplot showing the variations in eigengene values of the samples from the RAB11FIP5 module between the severe and mild/moderate groups. D Dot plot showing the top KEGG pathways with an FDR < 0.01 associated with the genes in the RAB11FIP5 module. The size of each dot on the x-axis represents the number of genes associated with each pathway. Additionally, the color of each dot corresponds to the respective FDR P-values linked to the pathways

We identified a total of 1135 hub CpGs (hCpGs) within the RAB11FIP5 module (Table S3). Upon examining the cluster of these hCpGs falling within specific (hub) genes, six hCpGs were observed within the CUX1 and ANK3 genes, four within the ZBTB38, four in the body region of the SSBP3, and four within GLI3, FOXP1, and AFF3. The genes associated with the most significant module hCpGs were primarily involved in pathways related to the PI3K-Akt, MAPK, Rap1, and Ras signaling pathways (Fig. 9D).

Discussion

We present novel evidence suggesting the potential involvement of DNAm in the severity of RSV disease among children younger than 2 years. To the best of our knowledge, this is the first attempt at using non-invasive buccal swab samples to investigate methylation biomarkers associated with varying degrees of severity following RSV infection. Despite the lack of extensive research on DNAm changes specifically related to RSV infection in saliva, the existing scientific literature provides compelling evidence that saliva, which contains both epithelial and immune cells, serves as the primary site for viral infection and subsequent initiation of the immune response during acute respiratory viral infections.

In our study, multiple CpG sites were identified in the ZBTB38 and TRIM6-TRIM34 genes in severe vs mild/moderate patients.

Upon analyzing the clinical history of the patients, we observed that severely ill children displayed more homogenous characteristics than those with mild/moderate symptoms. This pattern is clearly visible in the PCA built on DMPs of donors, where the mild/moderate group displayed greater variability,than the severe group. This difference could be attributed to the challenge of classifying patients as mild/moderate due to the diverse symptoms and age range.

Through DMP and DMR analysis, we identified two genes that were significantly differentially hypomethylated between severely ill children and those with only mild/moderate symptoms: ZBTB38 and TRIM6-TRIM34. Methylation patterns of both genes demonstrated high discriminatory power in distinguishing the two groups of infants, with ZBTB38 having an AUC of 0.984 and TRIM-6-TRIM-34 having an AUC of 0.906. These findings were further validated through pyrosequencing analysis in both the discovery and a new validation cohort, reinforcing the reliability of these methylation markers as prognostic indicators of RSV severity. ZBTB38, which contains seven statistically significant CpG sites, is a protein-coding gene that contains both zinc finger and BTB domains. It functions as a transcription factor and exhibits a strong affinity for methylated sequences in vitro, acting as a methyl-CpG binding protein. Recent research [34] has shown that ZBTB38 is hypomethylated and transcriptionally upregulated in adaptive NK cells. It binds to methylated CpG sites within DNA and plays a role in the negative regulation of apoptosis. This gene has also been implicated in asthma [35] and chronic obstructive pulmonary disease (COPD) [36], two conditions closely related to airway viral infections. TRIM6 is a member of the tripartite motif (TRIM) family, which is a group of proteins that are involved in a wide range of biological processes and are associated with various pathological conditions, including developmental disorders, neurodegenerative diseases, viral infection, and carcinogenesis. The TRIM6 gene plays a critical role in both the type I interferon (IFN-I) production and signaling pathways [37], and its absence inhibits IFN-I signaling leading to increased replication of interferon-sensitive viruses. A recent study on Ebola revealed that the knockout of TRIM6 reduces the replication of the virus, suggesting the gene’s significance as a host cellular factor for viral replication [38]. Similarly, the TRIM34 gene is significantly upregulated in response to IFN-I in macrophages, [39], and may limit HIV replication in a TRIM5-dependent manner [40]. In line with this evidence, the low methylation level observed in the DMPs belonging to this gene in severely ill children could suggest potential overexpression of this gene in this group of patients.

In our study, we observed that the DNAm levels at ZBTB38 gene in healthy controls were comparable to those in children with mild to moderate symptoms of RSV infection. This finding suggests that this gene exhibits a stable baseline methylation pattern under healthy conditions and mild RSV infection. The hypomethylation observed in severe RSV children for those CpGs implies that RSV infection may induce alterations in DNAm levels, particularly in those who develop severe manifestations of the disease. A general hypomethylation of CpGs in the whole blood of patients with bronchiolitis due to RSV or rhinovirus infections was confirmed by the study by Zhu et al. [16]. However, the specific DMPs identified in that study were not observed in our, possibly due to differences in tissue type.

Pathway analysis revealed intriguing enrichment in olfactory signaling, particularly the olfactory signaling pathway (R-HSA-381753), which is linked to genes differentially expressed during RSV infection [41]. Various animal models have demonstrated that certain viruses, such as Japanese encephalitis virus, influenza virus, and Herpes simplex virus, exploit the olfactory pathway to access the olfactory bulb (OB) and spread to other brain regions, including the hypothalamus and cortical areas [42]. While the mechanisms of RSV entry and spread in the CNS are not fully understood, a mouse model study in 2013 [43] demonstrated neuro-invasive ability of RSV through the olfactory pathway. More recently, Bryche and colleagues [44] detected RSV infection and replication in olfactory sensory neurons, leading to inflammation in the OB.

Through weighted correlation network analysis, a set of co-methylated positions with a significant negative correlation with severe RSV infection was identified. This suggests that the CpG blocks associated with these positions are hypomethylated, suggesting potential overactivation of these genes in severely infected children. The most important hub gene of the most significant module was RAB11FIP5, to which the hub position cg13662225 belongs. It has been reported that RAB11FIP plays a crucial role in regulating viral infections. For instance, RAB11FIP3 seems indispensable for the formation of filamentous virions [45]. RAB11FIP1 and RAB11FIP2 regulate the release of RSV [46], while RAB11FIP4 is essential for the movement of components related to cytomegalovirus (HCMV) [47]. The role of RAB11FIP5 in viral infection remained unclear until a recent study which [48] showed its correrlation to NK cell dysfunction, and the productionof broadly neutralizing antibodies (bnAbs) in HIV-1-infected individuals.

Among the other hCpGs (RAB11FIP5 module), four were located within the DMR of the previously described ZBTB38 gene. Additionally, four hub CpGs were found within the SSBP3 gene, whose expression tends to be upregulated in patients with RSV [49]. Another interesting gene is GLI3, which contains four hub CpGs that play essential roles in lung development and in the regulation of innate immune cells [50]. Furthermore, a cluster of hub CpGs was found within the ANK3 gene, encoding proteins that have been detected in various types of epithelial cells, bone marrow macrophages, and neurons. The upregulated expression of the ANK3 protein was reported [50] in alveolar macrophages of individuals with COPD, suggesting its potential involvement in lung function decline and COPD pathogenesis in at-risk smokers. Another gene associated with lung development regulation is FOXP1, which contains four CpGs that may contribute to the development of chronic lung diseases [51].

Further investigation of the RAB11FIP5 module revealed that the genes associated with the hub positions primarily participate in the PI3K-Akt signaling pathway. This pathway has been shown to enhance viral entry into host cells [52], even if, recent studies have proposed that it can serve as a “proviral” kinase upon activation but plays a role in the host’s antiviral response. The second most significant pathway was related to MAPK signaling, which is generally activated by various chemical and physical stimuli, including cytokines, hormones, growth factors, pathogens (including viruses), etc. [53]. The effect of MAPK signaling on virus replication varies depending on the virus in question, as it can either support or downregulate replication [54]. Furthermore, viruses can exploit independently activated MAPK signaling pathways to facilitate their own replication. These two top pathways identified through co-methylation analysis—PI3K-Akt and MAPK—are consistent with findings from a recent RNA-based epigenetic study by Zhu et al. [55], which linked nasal miRNAs during severe bronchiolitis to asthma-related pathways, including PI3K-Akt. These pathways are not only involved in acute viral responses but also play a key role in shaping long-term respiratory outcomes such as asthma and COPD.

In line with these findings, it is important to consider that, in addition to RSV infection, these infants may have a high abundance of pathogenic bacteria in their airways. This interplay between respiratory viral infections and microbiome could significantly impact long-term outcomes such as asthma. As highlighted by this recent study of bronchiolitis in infants (56) understanding how these factors interact will be crucial for future research in elucidating the complex mechanisms leading to chronic respiratory diseases.

Conclusions

The present study represents the first attempt to investigate the role of DNAm in young children with different severities of RSV infection using non-invasive samples. The findings of this study provide compelling evidence indicating that DNAm may indeed contribute to the severity of RSV disease in young children at the local level. The consistency of validation results highlights the robustness of the identified methylation signatures and their potential clinical utility in predicting disease severity. However, further large-scale studies are necessary to confirm these findings.

It is important to note that while our results indicate a correlation between RSV severity and DNAm changes, causation cannot be definitively established from this study alone. It remains uncertain whether the observed methylation changes were already present in the group of children and predisposed them to severe RSV infection or if RSV infection itself directly influenced their methylation levels. Future research should focus on longitudinal studies to elucidate the mechanisms underlying these methylation changes and their role in RSV pathogenesis.

Limitations of the study

One of the most important limitations of the present study is the small sample size of both the mild/moderate and severe groups. This limited size, combined with the extensive number of CpGs examined, increases the likelihood that some DNA probes did not reach the threshold for statistical significance at the minimum FDR level.

Another limitation pertains to age-related differences and their potential implications, as well as disparities between cohorts regarding previous episodes of bronchiolitis and other medical histories. We attempted to address these age variations through appropriate statistical analyses and by adjusting for age in all group comparisons. Future studies may benefit from more nuanced stratification of age cohorts or longitudinal analyses to better elucidate age-specific risk factors and outcomes. Additionally, it would be valuable to explore differential sex differences in disease severity in future studies using larger cohorts of both patients and controls.

Furthermore, the impact of previous episodes of bronchiolitis and other medical history on our study results warrants consideration. In this case, future analysis with a more homogeneous mild/moderate cohort is required to exclude any implications of previous events that could be implicated in DNAm changes.

Another limitation of the study is the unknown exact duration of symptoms before sample collection, which is a common issue in paedriatic research studies.

Despite these limitations, our findings could be validated in an independent cohort, and include promising candidate biomarkers related to RSV infection and severity and shed further light on the local mechanisms that might contribute to the severe phenotype, providing good proof of principle for non-invasive saliva samples.