Methylation Landscape of RUNX3 Promoter Region as a Predictive Marker for Th1/Th2 Imbalance in Bronchiolitis.

BACKGROUND The methylation status of RUNX3 promoter region, its impact on RUNX3 gene expression, and Th1/Th2 imbalance are unknown in bronchiolitis. This study aimed to explore the predictors of bronchiolitis developing into asthma. MATERIAL AND METHODS The methylation status of RUNX3 promoter was assessed using Illumina HiSeq platform method. The relative RUNX3 mRNA levels in PBMCs were measured by qRT-PCR. Serum IL-4 and IFN-γ concentrations were measured by ELISA. RESULTS A series of sites with significantly higher levels of methylation as compared to their corresponding controls were identified, including 24 sites in group Ba vs. group Cn, 13 sites in group Ba vs. group Ca, 7 sites in group Ba vs. group Bn, 16 sites in group Bn vs. group Cn, 11 sites in group Ca vs. group Cn, and 23 sites in group B vs. group C; P<0.05. The relative mRNA levels in group Ba were significantly lower than those in groups Cn, Ca, Bn; P<0.05. The serum IL-4 concentrations in group Ba were significantly higher than those in group Cn; P<0.05. The serum IFN-γ concentrations in group Ba were significantly lower than those in groups Cn, Ca, Bn; P<0.05. Correlation analysis showed that differentially methylated RUNX3 promoter region sites were significantly negatively correlated with levels of relative RUNX3 mRNA and IFN-γ, and were significantly positively correlated with IL-4 levels. CONCLUSIONS The methylation status of RUNX3 promoter region plays a role in Th1/Th2 imbalance by silencing RUNX3 gene expression, which can serve as predictive marker for the development of bronchiolitis into asthma.

Background

Bronchiolitis is a common lower respiratory tract infection disease caused by respiratory syncytial virus (RSV) in early life. Bronchiolitis due to RSV infection is the most common cause of infant’s first wheezing, and it is associated with recurrent wheezing and even asthma [1]. At present, it is generally believed that bronchiolitis and asthma have a similar pathogenesis. Previous studies have shown that both diseases appear to have improper balance and proportion of Th1/Th2 cells, accompanied by increased levels of Th2 cell cytokines and decreased levels of Th1 cell cytokines in the respiratory tract [2-4]. Th1 cells produce cytokines such as IFN-γ and IL-2, which can antagonize inflammation by activating T cells and phagocytes, inducing Th0 cells to transform into Th1 cells and inhibiting Th2 cells [4]. Th2 cells produce cytokines IL-4, IL-5, and IL-13, which can induce inflammation by promoting immunoglobulin E (IgE) secretion and sensitizing mastocytes and eosinophils [5-7]. Th1 cell inactivation and Th2 cell hyperactivation are crucial factors in the development of asthma [8,9]. While bronchiolitis is known as a risk factor for the development of asthma, only a relative some proportion of infants with bronchiolitis develop into asthma, and most infants can be cured after the acute period of treatment. The role and mechanism of Th1/Th2 cells imbalance in infants with bronchiolitis remain unclear.

Human runt-related transcription factor 3 (RUNX3), which is located in the short arm of the human first chromosome (1P36.1) with 2 large and highly conservative CpG islands in the total length of 67 KB, is a transcription factor of the runt domain family [10]. It is the key regulator of lineage-specific gene expression, cell identities, and functions [11]. Repression of the RUNX3 gene expression will influence various aspects of T cell-mediated immune pathologies [12], which can easily lead to human diseases. Recently, some studies reported that RUNX3 influences the differentiation of T cells, modulates the Th1/Th2 balance, and is a vital downstream target of transforming TGF-β [13,14]. In a murine model of asthma study, it was found that upregulation of RUNX3 gene expression attenuates allergic airway inflammation and airway hyperresponsiveness [15]. In contrast, downregulation of RUNX3 gene expression results in the occurrence of asthma with multiple mechanisms involved, such as Th1/Th2 imbalance and DC function abnormalities [13,16]. RUNX3 gene expression is controlled by 2 promoters, a proximal promoter and a distal promoter, located approximately 100 bp upstream of the transcription start site [17]. Studies have shown that methylation of the promoter region is closely related to silencing of RUNX3 gene expression [18,19]. Previous research has indicated that the proximal promoter consists of more CpG islands than does the distal promoter, so the proximal promoter has a relatively higher probability of being methylated than the distal promoter [20]. Methylation of the CpG islands in the RUNX3 gene promoter will result in silencing of the RUNX3 gene expression, which may lead to Th1/Th2 cells imbalance and eventually the development of asthma. Therefore, it is important to study the methylation of the RUNX3 promoter region in infants with bronchiolitis as it may develop into asthma in the future. However, methylation status of RUNX3 gene promoter region and its association with the gene expression levels of RUNX3 and the Th1/Th2 cells imbalance in bronchiolitis remain elusive.

In this study, we first used the Illumina HiSeq-based method to examine the methylation status of RUNX3 gene promoter region in PBMCs of infants with bronchiolitis. Then, the expression levels of RUNX3 mRNA were measured by quantitative real-time PCR (qRT-PCR). Finally, we measured serum IL-4 and IFN-γ concentrations by ELISA. We evaluated the significance of the methylation status of RUNX3 promoter region in infants with bronchiolitis through analysis of its correlation with mRNA expression levels of RUNX3 and concentrations of serum IL-4 and IFN-γ. Our findings revealed a strong association of the methylation status of RNUX3 gene promoter region with its mRNA expression and Th1/Th2 cells imbalance, and will provide important clinical evidence for the prediction of whether infants with bronchiolitis will develop into asthma.

Material and Methods

Patients and samples

In this study, 27 infants with bronchiolitis (group B) aged from 3 months to 2 years were recruited from the Lianyungang Maternal and Child Health Hospital from October 2017 to March 2018, including 13 cases with an atopic constitution (group Ba) and 14 cases with a non-atopic constitution (group Bn). Twenty-four age- and gender-matched healthy control infants (group C) were recruited from a physical examination center in the hospital, including 12 cases with an atopic constitution (group Ca) and 12 cases with a non-atopic constitution (group Cn). All infants with bronchiolitis met the diagnostic criteria for bronchiolitis in the 8th edition of ‘Zhu Futang Practice of Pediatrics’. The exclusion criteria were as follows: infants with endocrine, heart, liver, kidney, autoimmune diseases and other serious diseases. The clinical characteristics assessed at enrollment are shown in Table 1. The PBMCs and serum were collected from 51 infants, separated immediately and stored at −80°C.The separated PBMCs were divided into 2 parts, one for DNA extraction and the other for RNA extraction. All protocols used in this study were accepted by the patients and were approved by the Ethics Review Board of Lianyungang Maternal and Child Health Hospital.

Table 1

Study population demographics.

	Ba	Bn	Ca	Cn	P1 (Ba/Bn)	P2 (Ba/Ca)	P3 (Ba/Cn)	P4 (Bn/Cn)	P5 (Ca/Cn)	P6 (B/C)
n	13	14	12	12
Male	10 (76.9%)	10 (71.4%)	9 (69.2%)	6 (50%)	0.75	0.91	0.16	0.26	0.21	0.37
Female	3 (23.1%)	4 (28.6%)	3 (31.8%)	6 (50%)
Age (months)	12.46±8.35	8.64±4.38	10.75±6.98	15±13.87	0.15	0.59	0.58	0.12	0.35	0.29
Weight (kg)	2.97±0.50	3.32±0.51	3.21±0.55	3.17±0.46	0.08	0.26	0.32	0.43	0.84	0.81
Caesarean	4 (30.7%)	9 (64.3%)	5 (41.7%)	5 (41.7%)	0.08	0.57	0.57	0.25	1	0.64
Breast Feeding	6 (46.2%)	6 (42.9%)	6 (50%)	8 (66.7%)	0.86	0.85	0.30	0.23	0.41	0.32
Premature	1 (7.7%)	1 (7.1%)	0 (0)	0 (0)	0.96	0.33	0.33	0.35	–	0.17
Passive smoking	8 (61.5%)	8 (57.1%)	8 (66.7%)	7 (58.3%)	0.82	0.79	0.87	0.95	0.67	0.81
Paroxysmal continuous cough	4 (30.8%)	6 (42.9%)	–	–	0.52	–	–	–	–	–
Paroxysmal solitary cough	9 (69.2%)	8 (57.1%)	–	–
No fever	5 (38.5%)	5 (35.7%)	–	–	0.09	–	–	–	–	–
Moderate fever	6 (46.2%)	2 (14.3%)	–	–
Ardent fever	2 (15.4%)	7 (50%)	–	–
Disease course (<7 days)	8 (61.5%)	9 (64.3%)	–	–	0.88	–	–	–	–	–
Disease course (>7 days)	5 (38.5%)	5 (35.7%)	–	–

B – represents bronchiolitis; Ba – represents atopic constitution bronchiolitis; Bn – represents non-atopic constitution bronchiolitis; C – represents heathy control; Ca – represents atopic constitution heathy control; Cn – represents non-atopic constitution heathy control.

Extraction of genomic DNA and total RNA from PBMCs

Genomic DNA from PBMCs was isolated using a DNA Cell Mini Kit (Tiangen, Dalian, China) according to the manufacturer’s instructions and was then used for methylation detection. Total RNA from PBMCs was isolated using TRIzol reagent (Invitrogen, California, USA); after reverse transcription, it was then used for qPCR. The DNA and RNA were quantified by a NanoDrop™ ND-2000 spectrophotometer (NanoDrop, DE, USA).

Bisulfite modification and targeted bisulfite pyrosequencing assay

Genomic DNA conversion via bisulfite modification using the EZ DNA Methylation™ Kit (ZYMO, CA, USA) was performed on 2 μg of extracted DNA. Based on the positions of the candidate CpG sites for RUNX3 (NM_004350), we designed 16 sequential pairs of primers by Primer 5.0. The primers and primer composition and concentration of the Multiplex PCR Panel were optimized by capillary electrophoresis. The primer sequence information is shown in Table 2, and the gene sequence amplification strategy is shown in Figure 1. Multiplex PCR was performed to amplify the targeted DNA sequence. Then, the targeted DNA fragments were sequenced by Illumina HiSeq 2000 (San Diego, CA, USA). BS-Seeker2 was used for mapping bisulfite-treated reads as well as for methylation calling and for analyzing the bisulfite sequencing results [21].

Table 2

The methylation primer sequence informations.

Number	Forward primer (5′→3′)	Reverse primer (5′→3′)	Size
1	GTYGTGGTTAGGTAGGGTTTTG	CAAAAATCCCTCATTCTCTAAAAACT	215 bp
2	AAAGGAAGGAYGGTAAAGGTTTTT	AACRAATCATAATCCCRCAACCT	263 bp
3	GTTTTAGGGGYGAAGGATTTTG	AAATAAAAACRCCCAAACTAAATCACTC	252 bp
4	TTTTGGAGAYGGGGTTAGATT	AAATCCCRAAATAAATAAAAACCAACTC	254 bp
5	AGTAGATTTAGGGGTYGGTTATAGTTG	CRCTAAACAAAAATCAAATCCATAATAAAATC	267 bp
6	TTAAAATAAGTTTTAGGATTTTGTAGGTT	CATAAAAAAAAAAAATTCTCTCCTCTAAA	237 bp
7	GGGGTTTTAGGGTYGGTGA	AAAACTCCCTTCCRCCTATCC	222 bp
8	TGTGGGGTGGAGGTTTTT	CAACCCCTTCRCTCCTCCT	200 bp
9	GTAGTTTTAGAATAAATTTTTTAGAATTAAGTGG	CCCTAAACTATAACACTAAAACCTCCTC	246 bp
10	GGGTTTTAGGGYGTAGGG	ACCCCRCCACTTAATTCTAA	290 bp
11	GGGGAGGGAGGTGTGAAG	AACRAACCCTAACRAATAATAACC	228 bp
12	GTTTYGGGTTTYGTATTTATTTTGAAG	ACCTCCCTCCCCRACCTT	261 bp
13	GGGGTAATTTYGTTAAAAGTTATTTTTGTAG	CCCTCRCACTAACRCTACAACAAA	263 bp
14	TTYGTTATTYGTAGGGTTGTATTTGAG	AACCAAAACCTCRAAAACAATAAAA	269 bp
15	GTTYGAGGGTGTGAGTAGGAAGT	AAATTTTACRCCCCCCTAAC	184 bp
16	ATGGTTGGGGGTTGTTAGG	CCTCCAACACCCCTCCTC	191 bp

Figure 1

RUNX3 promoter region CpG sites amplification map.

Quantification of RUNX3 mRNA level

For determination of relative RUNX3 mRNA levels, we used a relative quantitative method. Total RNA was used for reverse-transcription reactions by the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China) to synthesize cDNA. The obtained cDNA was amplified by qRT-PCR using SYBR Premix Ex Taq™ II (TaKaRa, Dalian, China). To normalize the results, OAZ1 was used as a reference gene for RUNX3 detection. The reaction conditions for qRT-PCR were 95°C for 30 s, followed by 40 cycles of 95°C for 5 s, optimal annealing temperature 60°C for 31 s, and 72°C for 10 s. The relative RUNX3 mRNA levels were calculated by the 2−ΔΔCT method. The primer sequence information for RUNX3 and OAZ1 is shown in Table 3.

Table 3

The primer sequence informations of RUNX3 and OAZ1.

Number	Forward primer (5′→3′)	Reverse primer (5′→3′)	Size
RUNX3	GTCGCCTTCAAGGTGGTGGC	TAGGTCGCCACTTGGGTGGG	227 bp
OAZ1	AGCAAGGACAGCTTTGCAGTTCTC	GATGCCCCGGTCTCACAATC	154 bp

Measurement of serum cytokines

IFN-γ and IL-4 are representative cytokines of Th1 and Th2 cells, respectively, which can reflect the Th1/Th2 balance. The serum IL-4 and IFN-γ concentrations in the bronchiolitis (Ba and Bn) and healthy control (Ca and Cn) groups were detected using a human IL-4 ELISA kit and a human IFN-γ ELISA kit according to the manufacturer’s protocol (Neobioscience, Shenzhen, China). The OD450 nm value was detected by a microplate reader (Labsystems Dragon, USA). The serum IL-4 and IFN-γ concentrations were calculated according to standard curves.

Statistical analysis

We applied the Wilcoxon rank-sum test for the differential methylation status in each CpG site. The FDR (false discovery rate) was applied for multiple test correction. The differential analysis of relative RUNX3 mRNA levels, DNA methylation levels, and the comparative analysis of serum IL-4 and IFN-γ concentrations were assessed using ANOVA and the SNK-q test. Correlations were performed using Pearson’s correlation test. All statistical analyses were conducted with SPSS 22.0 software and R version 3.2.1. P<0.05 was considered significantly different.

Results

Methylation status of the RUNX3 promoter region in PBMCs

An Illumina HiSeq 2000 was used to analyze the methylation state of the RUNX3 gene promoter region in PBMCs of infants with bronchiolitis (group B) and healthy control (group C) infants. It was found that for 358 CpG sites in the RUNX3 gene promoter region, the hypermethylation status occurred mainly in the RUNX3 gene upstream promoter, and the hypomethylation status occurred mainly in the RUNX3 gene downstream promoter. The methylation heatmap for the 4 groups is shown in Figure 2. The Wilcoxon rank-sum test was used to detect the difference at each CpG site among the 4 groups. We found that 24 methylation sites in group Ba exhibited significantly higher mean methylation levels than those in group Cn, 13 methylation sites in group Ba had higher mean methylation levels compared to those in group Ca, 7 methylation sites in group Ba had higher mean methylation levels compared to those in group Bn, 16 methylation sites in group Bn had higher mean methylation levels compared to those in group Cn, 11 methylation sites in group Ca had higher mean methylation levels compared to those in group Cn, and 23 methylation sites in group B had higher mean methylation levels compared to those in group C. All P values were less than 0.05 (P<0.05), and the results are shown in Supplementary Tables 1–6.

Figure 2

RUNX3 promoter region methylation heatmap.

Analysis of RUNX3 relative mRNA levels in PBMCs

In the 4 groups of PBMCs, relative RUNX3 mRNA levels were measured by the 2−ΔΔCT method, and the OAZ1 gene was used as a reference gene. We found that the relative RUNX3 mRNA levels in PBMCs of infants with bronchiolitis (group B) were significantly lower than those in healthy control infants (group C), (P<0.05); in particular, those in the atopic infants with bronchiolitis (group Ba) were significantly lower than those in atopic healthy control infants (group Ca), non-atopic healthy control infants (group Cn) and non-atopic infants with bronchiolitis (groups Bn), (P<0.05). These results suggest that the relative RUNX3 mRNA levels were significantly decreased in PBMCs from infants with bronchiolitis, especially for those with atopic constitution. The results are shown in Figure 3.

Figure 3

RUNX3 relative mRNA levels in PBMCs, with OAZ1 as reference gene (* P<0.05).

Analysis of serum IFN-γ and IL-4 concentrations

Two cytokine concentrations in the 4 groups were measured by ELISA. The serum IL-4 and IFN-γ concentrations in infants with bronchiolitis (group B) were significantly different (P<0.05) from those in healthy control infants (group C). The serum IL-4 concentration in infants with bronchiolitis (group B) was significantly higher than that in healthy control infants (group C) (P<0.05); in particular, serum IL-4 concentration in atopic infants with bronchiolitis (group Ba) was significantly higher than that in non-atopic healthy control infants (group Cn) (P<0.05), as shown in Figure 4A. The serum IFN-γ concentration in infants with bronchiolitis (group B) was significantly lower than that in healthy control infants (group C) (P<0.05); in particular, serum IFN-γ concentration in atopic bronchiolitis infants (group Ba) was significantly lower than that in atopic and non-atopic healthy control infants and non-atopic infants with bronchiolitis (groups Ca, Cn and Bn) (P<0.05), as shown in Figure 4B. These results suggest that the serum IL-4 concentration was significantly increased in PBMCs of infants with bronchiolitis, especially for those with atopic constitution, while the IFN-γ concentration showed the opposite result.

Figure 4

(A) Serum IL-4 concentration in each group. (B) Serum IFN-γ concentration in each group (* P<0.05).

The correlation between relative RUNX3 mRNA levels and DNA methylation

It is widely accepted that promoter methylation can strongly influence gene expression. To evaluate the correlation between differential methylation of RUNX3 gene promoter region sites and relative RUNX3 mRNA levels, Pearson’s correlation test was applied. For comparison of groups Ba/Cn, the average methylation levels of 24 differentially methylated CpG sites were significantly negatively correlated with relative RUNX3 mRNA levels (Figure 5A; P<0.05). For comparison of groups Ba/Ca, the average methylation levels of 13 differentially methylated CpG sites were negatively correlated with relative RUNX3 mRNA levels, but did not reach significance (Figure 5B). For comparison of groups Ba/Bn, the average methylation levels of 7 differentially methylated CpG sites were significantly negatively correlated with relative RUNX3 mRNA levels (Figure 5C; P<0.05). For comparison of groups Bn/Cn, the average methylation levels of 16 differentially methylated CpG sites were significantly negatively correlated with relative RUNX3 mRNA levels (Figure 5D; P<0.01). For comparison of groups Ca/Cn, the average methylation levels of 11 differentially methylated CpG sites were significantly negatively correlated with relative RUNX3 mRNA levels (Figure 5E; P<0.05). For comparison of groups B/C, the average methylation levels of 23 differentially methylated CpG sites were significantly negatively correlated with relative RUNX3 mRNA levels (Figure 5F; P<0.01). These results suggest that the methylation levels of the RUNX3 promoter region are negatively associated with its gene expression.

Figure 5

The correlation between RUNX3 promoter region methylation levels and relative mRNA levels. (A) Represents the result of comparison between group Ba and group Cn. (B) Represents the result of comparison between group Ba and group Ca. (C) Represents the result of comparison between group Ba and group Bn. (D) Represents the result of comparison between group Bn and group Cn. (E) Represents the result of comparison between group Ca and group Cn. (F) represents the result of comparison between group B and group C.

The correlation between the methylation of RUNX3 promoter region and serum cytokine levels

In addition, we also evaluated the correlation between differential methylation of RUNX3 gene promoter region sites and serum cytokine concentrations, including IL-4 and IFN-γ, in the serum of 27 infants with bronchiolitis and 24 control infants by ELISA. Pearson’s correlation test was conducted to assess the correlation between serum IL-4 and IFN-γ levels and differential methylation of the CpG sites. For comparison of groups Ba/Cn, the average methylation levels of 24 differentially methylated CpG sites were significantly positively correlated with serum IL-4 concentration (Figure 6A; P<0.05) and significantly negatively correlated with serum IFN-γ concentration (Figure 6A; P<0.05). For comparison of groups Ba/Bn, the average methylation levels of 7 differentially methylated CpG sites were significantly positively correlated with serum IL-4 concentration (Figure 6B; P<0.05) and significantly negatively correlated with serum IFN-γ concentration (Figure 6B; P<0.05). For comparison of groups Ba/Ca, the average methylation levels of 13 differentially methylated CpG sites were positively correlated with serum IL-4 concentration and negatively correlated with serum IFN-γ concentration, but the correlations were not significant (Figure 6C). For comparison of groups Bn/Cn, the average methylation levels of 16 differentially methylated CpG sites were significantly positively correlated with serum IL-4 concentration (Figure 6D; P<0.05) and significantly negatively correlated with serum IFN-γ concentration (Figure 6D; P<0.05). For comparison of groups Ca/Cn, the average methylation levels of 11 differentially methylated CpG sites were significantly positively correlated with serum IL-4 concentration (Figure 6E; P<0.05) and significantly negatively correlated with serum IFN-γ concentration (Figure 6E; P<0.01). For comparison of groups B/C, the average methylation levels of 23 differentially methylated CpG sites were significantly positively correlated with serum IL-4 concentration (Figure 6F; P<0.01) and significantly negatively correlated with serum IFN-γ concentration (Figure 6F; P<0.01). These results indicate that the methylation states of RUNX3 promoter region affected the Th1/Th2 balance in vivo.

Figure 6

The correlation between RUNX3 promoter region methylation levels and relative mRNA levels. (A) Represents the result of comparison between group Ba and group Cn. (B) Represents the result of comparison between group Ba and group Ca. (C) Represents the result of comparison between group Ba and group Bn. (D) Represents the result of comparison between group Bn and group Cn. (E) Represents the result of comparison between group Ca and group Cn. (F) Represents the result of comparison between group B and group C.

Discussion

Epigenetic modification refers to hereditary alteration of gene expression without changes in DNA sequence, which can be transmitted steadily during development and cell proliferation and includes DNA methylation, histone acetylation, and RNA-associated silencing [22]. Among these modifications, DNA methylation is the most common and important mechanism influencing gene expression, and it is an enzyme-mediated chemical modification that occurs after birth with the methylation of the 5′ end of cytosine [23,24]. It is important for normal development of mammalian genes and embryos, during which gene promoter methylation inhibits gene expression [25]. Research has shown that methylation of the 5′-CpG islands in the RUNX3 gene promoter is a major gene-silencing mechanism [17]. Therefore, our study was based on whether RUNX3 gene promoter methylation influences its gene expression in bronchiolitis. For DNA methylation detection, we used the high-accuracy Illumina HiSeq platform method to ensure the reliability of the detection. The data showed that some CpG sites in the RUNX3 gene promoter region were significantly methylated in PBMCs from subjects with bronchiolitis, especially in subjects with an atopic constitution (group Ba, high-risk asthma infants), and the methylation levels of CpG sites were significantly higher in subjects with bronchiolitis than in the 2 control groups (Ca and Cn). These data indicate that RUNX3 gene promoter region methylation plays an important role in bronchiolitis. To our knowledge, this is the first study revealing the methylation status of the RUNX3 promoter region in different constitutions of bronchiolitis.

The RUNX3 gene has 2 promoters. The distal promoter does not contain a CpG island, and the proximal promoter has a CpG island with binding sites for Sp1, and Sp1 elements located upstream and downstream of the transcription start site that have been shown to protect CpG islands from methylation [20,26]. Our results showed that the hypermethylation sites were mainly located in the upstream of the promoter, whereas the downstream sites were mostly maintained in the hypomethylation states. In accordance, most of the differentially methylated sites identified in comparison between sub-groups were located in the upstream hypermethylation region. These results are essentially consistent with previous studies on the methylation of RUNX3 promoter region [17,26].

Since the relative RUNX3 mRNA levels in PBMCs of infants with bronchiolitis were significantly reduced, we wondered whether the significantly methylated CpG sites in RUNX3 gene promoter influenced its gene expression in these infants, as in other diseases [27-29]. Therefore, we performed a correlation analysis between the RUNX3 promoter region methylation levels and the gene expression levels. RUNX3 promoter region methylation levels were negatively correlated with relative RUNX3 mRNA levels. In the pathogenesis of allergic respiratory diseases such as asthma, it has been confirmed that the decrease in RUNX3 gene expression is caused by RUNX3 DNA methylation [30]. RUNX3 promoter region methylation, as an epigenetic regulation mechanism, resulted in the decrease in RUNX3 expression in bronchiolitis.

Many studies have confirmed that DNA methylation plays an important role in the mechanism of infant asthma by influencing gene expression [31,32]. Our results confirmed a negative correlation between RUNX3 methylation and gene expression in bronchiolitis. Therefore, we evaluated the relationship between RUNX3 promoter region methylation and serum cytokine levels. We found that RUNX3 promoter region methylation was positively correlated with serum IL-4 levels and negatively correlated with IFN-γ levels. This result further supported previous studies of the relationship between serum cytokines and RUNX3 gene expression. Therefore, we conclude that when the RUNX3 gene promoter region is significantly methylated, its mRNA expression will decrease, which may influence the secretion and serum levels of IL-4 and IFN-γ and T lymphocyte differentiation (Th1/Th2 imbalances), and ultimately lead to bronchiolitis.

Based on our present results, we conclude that RUNX3 promoter region methylation can play a role in the development of bronchiolitis into asthma. However, there are still some limitations in our study. The research was carried out in PBMCs. Although PBMCs can be used as a representative of methylation changes in many diseases [33-35], respiratory epithelial cells may be better to confirm our findings. In addition, due to the small number of samples included in this study, the results may not be exactly accurate. Also, more scientific methods to predict biomarkers should be employed, such as logit-lapnet method or network-regularized cox model [36,37]. The inclusion of respiratory epithelial cells, larger sample sizes and more scientific methods in later studies will yield more reliable results.

Conclusions

RUNX3 is a promising specific marker for the development of bronchiolitis into asthma, especially in diagnosis and disease prognosis. The methylation landscape and significance of the RUNX3 promoter regions that can be used in the risk assessment of bronchiolitis and provide useful information for genetic diagnosis and treatment. However, this study had a small sample size and was not used respiratory epithelial cells. Therefore, further study is need in the future.

Supplementary Table 1

The mean methylation level analysis of methylation sites in group Ba and group Cn.

	Chromosomal site	Ba	Cn	P (Ba/Cn)	FDR
1	Chr1: 25258221	0.977±0.011	0.961±0.015	0.006000	1.029×10−2*
2	Chr1: 25258264	0.756±0.052	0.689±0.037	0.001331	3.194×10−3**
3	Chr1: 25258257	0.841±0.030	0.788±0.034	0.000611	1.833×10−3**
4	Chr1: 25258289	0.802±0.033	0.763±0.029	0.004678	9.356×10−3**
5	Chr1: 25258384	0.912±0.012	0.854±0.034	0.000113	7.000×10−4**
6	Chr1: 25258426	0.809±0.030	0.741±0.026	0.000090	7.000×10−4**
7	Chr1: 25258402	0.915±0.014	0.872±0.021	0.000140	7.000×10−4**
8	Chr1: 25258404	0.753±0.023	0.691±0.031	0.000175	7.000×10−4**
9	Chr1: 25258406	0.965±0.007	0.938±0.017	0.000907	2.419×10−3**
10	Chr1: 25258412	0.840±0.019	0.785±0.033	0.000268	9.189×10−4**
11	Chr1: 25258416	0.946±0.011	0.909±0.022	0.000175	7.000×10−4**
12	Chr1: 25258418	0.883±0.015	0.819±0.030	0.000072	7.000×10−4**
13	Chr1: 25258479	0.947±0.009	0.930±0.018	0.019341	2.443×10−2*
14	Chr1: 25258485	0.872±0.024	0.838±0.037	0.025740	3.089×10−2*
15	Chr1: 25258546	0.707±0.059	0.629±0.075	0.012347	1.976×10−2*
16	Chr1: 25258531	0.917±0.021	0.880±0.042	0.029577	3.380×10−2*
17	Chr1: 25258557	0.482±0.057	0.420±0.066	0.038471	4.043×10−2*
18	Chr1: 25258613	0.191±0.027	0.151±0.030	0.005528	1.021×10−2*
19	Chr1: 25258623	0.340±0.044	0.282±0.043	0.004678	9.356×10−3**
20	Chr1: 25258633	0.169±0.028	0.135±0.028	0.015074	2.128×10−2*
21	Chr1: 25258654	0.728±0.058	0.677±0.069	0.014200	2.128×10−2*
22	Chr1: 25258673	0.811±0.021	0.788±0.027	0.044164	4.416×10−2*
23	Chr1: 25258679	0.782±0.042	0.730±0.068	0.038741	4.043×10−2*
24	Chr1: 25257392	0.472±0.017	0.441±0.027	0.016698	2.226×10−2*

Ba – represents atopic constitution bronchiolitis; Cn – represents non-atopic constitution heathy control.

P<0.05,

P<0.01.

Supplementary Table 2

The mean methylation level analysis of methylation sites in group Ba and group Ca.

	Chromosomal site	Ba	Ca	P (Ba/Ca)	FDR
1	Chr1: 25258136	0.985±0.012	0.973±0.009	0.007681	1.248×10−2*
2	Chr1: 25258264	0.756±0.052	0.721±0.018	0.025740	3.042×10−2*
3	Chr1: 25258257	0.841±0.030	0.813±0.013	0.005537	1.028×10−2*
4	Chr1: 25258384	0.912±0.012	0.874±0.016	0.000036	4.680×10−4**
5	Chr1: 25258426	0.809±0.030	0.769±0.026	0.004670	1.012×10−2*
6	Chr1: 25258402	0.915±0.014	0.892±0.016	0.001606	5.220×10−3**
7	Chr1: 25258404	0.753±0.023	0.716±0.021	0.000611	3.972×10−3**
8	Chr1: 25258406	0.965±0.007	0.957±0.008	0.016698	2.171×10−2*
9	Chr1: 25258412	0.840±0.019	0.810±0.020	0.002319	6.029×10−3**
10	Chr1: 25258416	0.946±0.011	0.930±0.017	0.033895	3.390×10−2*
11	Chr1: 25258418	0.883±0.015	0.852±0.021	0.001100	4.767×10−3**
12	Chr1: 25258613	0.191±0.027	0.166±0.019	0.033861	3.390×10−2*
13	Chr1: 25257629	0.834±0.030	0.804±0.026	0.016698	2.171×10−2*

Ba – represents atopic constitution bronchiolitis; Ca – represents atopic constitution heathy control.

P<0.05,

P<0.01

Supplementary Table 3

The mean methylation level analysis of methylation sites in group Ba and group Bn.

	Chromosomal site	Ba	Bn	P (Ba/Bn)	FDR
1	Chr1: 25258384	0.912±0.012	0.882±0.023	0.000262	1.386×10−3**
2	Chr1: 25258416	0.951±0.012	0.934±0.013	0.022564	2.256×10−2*
3	Chr1: 25258426	0.809±0.030	0.776±0.031	0.002620	6.113×10−3**
4	Chr1: 25258404	0.753±0.023	0.714±0.027	0.000396	1.386×10−3**
5	Chr1: 25258406	0.965±0.007	0.953±0.014	0.017417	2.032×10−2*
6	Chr1: 25258412	0.840±0.019	0.813±0.022	0.004885	8.549×10−3**
7	Chr1: 25258418	0.883±0.015	0.861±0.022	0.013330	1.866×10−2*

Ba – represents atopic constitution bronchiolitis; Bn – represents atopic constitution heathy control.

P<0.05,

P<0.01.

Supplementary Table 4

The mean methylation level analysis of methylation sites in group Bn and group Cn.

	Chromosomal site	Bn	Cn	P (Bn/Cn)	FDR
1	Chr1: 25258264	0.729±0.049	0.688±0.037	0.042153	4.215×10−2*
2	Chr1: 25258257	0.822±0.030	0.788±0.034	0.023629	3.926×10−2*
3	Chr1: 25258402	0.904±0.019	0.872±0.021	0.001704	1.284×10−2*
4	Chr1: 25258404	0.714±0.027	0.691±0.031	0.039650	4.215×10−2*
5	Chr1: 25258406	0.953±0.014	0.938±0.017	0.034961	4.215×10−2*
6	Chr1: 25258412	0.812±0.022	0.785±0.033	0.017982	3.926×10−2*
7	Chr1: 25258416	0.934±0.013	0.909±0.022	0.002408	1.284×10−2*
8	Chr1: 25258418	0.861±0.022	0.819±0.030	0.002028	1.284×10−2*
9	Chr1: 25258426	0.776±0.031	0.741±0.026	0.003973	1.589×10−2*
10	Chr1: 25258479	0.971±0.011	0.930±0.018	0.034961	4.215×10−2*
11	Chr1: 25258537	0.960±0.012	0.946±0.015	0.023629	3.926×10−2*
12	Chr1: 25258568	0.825±0.037	0.781±0.060	0.026989	3.926×10−2*
13	Chr1: 25258633	0.162±0.030	0.135±0.028	0.039650	4.215×10−2*
14	Chr1: 25257392	0.467±0.025	0.441±0.027	0.015631	3.926×10−2*
15	Chr1: 25257229	0.127±0.017	0.114±0.015	0.026936	3.926×10−2*
16	Chr1: 25257217	0.117±0.016	0.104±0.014	0.026936	3.926×10−2*

Bn – represents non-atopic constitution bronchiolitis; Cn – represents non-atopic constitution heathy control.

P<0.05.

Supplementary Table 5

The mean methylation level analysis of methylation sites in group Ca and group Cn.

	Chromosomal site	Ca	Cn	P (Ca/Cn)	FDR
1	Chr1: 25258264	0.721±0.018	0.689±0.037	0.011074	3.045×10−2*
2	Chr1: 25258257	0.813±0.013	0.788±0.034	0.024343	3.347×10−2*
3	Chr1: 25258289	0.786±0.027	0.763±0.029	0.043308	4.331×10−2*
4	Chr1: 25258390	0.982±0.003	0.978±0.005	0.017901	3.282×10−2*
5	Chr1: 25258426	0.769±0.026	0.741±0.026	0.007912	2.901×10−2*
6	Chr1: 25258402	0.892±0.016	0.872±0.021	0.015314	3.282×10−2*
7	Chr1: 25258404	0.716±0.021	0.691±0.031	0.020921	3.288×10−2*
8	Chr1: 25258406	0.957±0.008	0.938±0.017	0.003235	2.901×10−2*
9	Chr1: 25258412	0.810±0.020	0.785±0.033	0.037667	4.143×10−2*
10	Chr1: 25258416	0.930±0.017	0.909±0.022	0.037667	4.143×10−2*
11	Chr1: 25258418	0.852±0.021	0.819±0.030	0.007912	2.901×10−2*

Ca – represents atopic constitution heathy control; Cn – represents non-atopic constitution heathy control.

P<0.05.

Supplementary Table 6

The mean methylation level analysis of methylation sites in group B and group C.

	Chromosomal site	B	C	P (B/C)	FDR
1	Chr1: 25258264	0.742±0.051	0.705±0.032	0.007577	1.793×10−2*
2	Chr1: 25258257	0.831±0.031	0.801±0.028	0.000958	2.947×10−3**
3	Chr1: 25258384	0.897±0.023	0.864±0.028	0.000068	5.213×10−4**
4	Chr1: 25258426	0.792±0.034	0.755±0.029	0.000217	1.248×10−3**
5	Chr1: 25258402	0.909±0.018	0.882±0.021	0.000028	5.213×10−4**
6	Chr1: 25258404	0.733±0.032	0.704±0.029	0.000958	2.947×10−3**
7	Chr1: 25258406	0.959±0.013	0.948±0.016	0.009208	1.925×10−2*
8	Chr1: 25258412	0.826±0.024	0.797±0.030	0.001025	2.947×10−3**
9	Chr1: 25258416	0.940±0.013	0.920±0.022	0.000553	2.544×10−3**
10	Chr1: 25258418	0.872±0.022	0.835±0.030	0.000050	5.213×10−4**
11	Chr1: 25258479	0.945±0.009	0.935±0.017	0.018329	3.011×10−2*
12	Chr1: 25258485	0.866±0.026	0.846±0.036	0.028559	3.285×10−2*
13	Chr1: 25258546	0.687±0.066	0.648±0.074	0.041540	4.154×10−2*
14	Chr1: 25258531	0.914±0.022	0.889±0.038	0.024724	3.285×10−2*
15	Chr1: 25258537	0.958±0.012	0.949±0.015	0.023540	3.285×10−2*
16	Chr1: 25258557	0.473±0.057	0.437±0.068	0.037909	3.963×10−2*
17	Chr1: 25258568	0.823±0.040	0.794±0.057	0.029992	3.285×10−2*
18	Chr1: 25258613	0.181±0.028	0.159±0.026	0.013429	2.574×10−2*
19	Chr1: 25258623	0.327±0.043	0.297±0.048	0.028592	3.285×10−2*
20	Chr1: 25258633	0.165±0.029	0.144±0.033	0.023450	3.285×10−2*
21	Chr1: 25258673	0.809±0.017	0.792±0.023	0.017417	3.011×10−2*
22	Chr1: 25257629	0.827±0.034	0.812±0.026	0.029992	3.285×10−2*
23	Chr1: 25257392	0.469±0.030	0.445±0.032	0.007794	1.793×10−2*

B – represents bronchiolitis; C – represents heathy control.

P<0.05,

P<0.01.