A total of 6 female BALB/c mice (8 weeks old and weight 20–25 g) were obtained from Daehan BioLink Co., Ltd., and were housed under standard laboratory conditions (22±2°C, 60% relative humidity, 12/12 h light/dark cycle, and food and water ad libitum) for at least 1 week before the experiment began. Ethical and scientific management procedures for all animals were approved by The Institutional Animal Care and Use Committee of the Korea Institute of Oriental Medicine (approval no. 17-073).

Mice were randomly divided into two groups: Vehicle-treated (n=3) and OVA-treated (n=3). Asthma was induced using OVA according to the previously described experimental procedures (19). To boost immune responses, mice were injected intraperitoneally with OVA (50 µg) emulsified with 200 µl aluminum sulfate (InvivoGen) twice, 7 days apart (on day 0 and 7). The mice were administered intranasally with 50 µl OVA (25 µg) from day 12 to 15 under anesthesia with 2% isoflurane (Piramal Critical Care, Ltd.).

Pathophysiological evidence of OVA-induced asthma was indicated by airway inflammation, eosinophil and interleukin (IL)-4 production (Fig. S1). The right lungs of the mice were fixed with 10% neutral buffered formalin at 25°C for 24 h, embedded in paraffin, sectioned into 4-µm thick slices, and stained with hematoxylin for 5 min and eosin 3 min (Sigma-Aldrich; Merck KGaA) at 25°C to evaluate inflammatory responses in lung tissues, as well as with periodic acid-Schiff (IMEB Inc.) to evaluate mucus production. A tracheostomy was performed in mice anesthetized with 85 mg/kg alfaxan (Alfaxalone^®; Jurox Pty Ltd.) to obtain bronchoalveolar lavage fluid (BALF). The lungs were lavaged twice with cold PBS (total volume 1.4 ml) and the obtained BALF was centrifuged at 380 × g for 5 min at 4°C. The supernatant was used for an ELISA assay (cat. no. M400B; R&D Systems, Inc.) in order to analyze level of IL-4, according to the manufacturer's protocol. The pellets of BALF were used to count inflammatory cells in the BALF. The eosinophil count was determined using a cell counting machine (Countess II™ Automated Cell Counter, Thermo Fisher Scientific, Inc.).

Statistical analysis was performed using Student's t-test with Prism 8 software (GraphPad Software, Inc.). P<0.05 was considered to indicate a statistically significant difference.

A schematic diagram for identifying target genes involved in asthma development is presented in Fig. 1. Gene expression profile data was used from a previous study by Baek et al (18), which were deposited in Gene Expression Omnibus (GSE114587; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE114587). A total of six chips were available for further analysis, including samples from three vehicle- and three OVA-treated animals. Detailed information is available from the GEO database (18).

Total RNA was extracted from the lung tissues of saline-treated (normal) and OVA-induced mice. RNA isolation was performed using the PureLink™ RNA Mini kit (cat. no. 12183018A; Thermo Fisher Scientific, Inc.). RNA-seq libraries were prepared using the TruSeq RNA Sample Prep kit (Illumina, Inc.) and sequenced with a NextSeq^® 500 System Whole-Genome Sequencing Solution (Illumina, Inc.) using the 76-bp paired-end reads. The raw FASTQ reads were trimmed to remove adapters and low-quality reads (per-base quality <20) using cutadapt (v.1.10, http://cutadapt.readthedocs.io/en/stable/) (20), then the high-quality sequence reads were aligned to the Mus musculus genome (mm10) using HISAT2 (v.2.1.0) and StringTie (v.1.3.4, http://ccb.jhu.edu/software/stringtie/) (21). Gene expression was quantified using the ballgown package in R (v.2.6.0) (22). Differentially expressed genes (DEGs) between the normal and OVA-induced mice groups were identified using the edgeR package (v.3.16.5) (23) based on negative binomial models of RNA-seq count data. Candidate DEGs were filtered using log₂(FC) ≥|1| and P<0.05 as thresholds, where FC indicates the fold change. DAVID was used to functionally annotate the genes (Tables SI and SII) (24,25).

DNA methylation profile data from the previous study by Baek et al (18) was used. DNA libraries were prepared according to the SureSelectXT Methyl-Seq Target Enrichment System protocol (Agilent Technologies, Inc.) and sequenced using an Illumina HiSeq 2500 platform (Illumina, Inc.) to generate 101-bp paired-end reads. After sequencing, the raw FASTQ reads were trimmed to remove adapters and low-quality reads (per-base quality <20) using cutadapt (v.1.10) and mapped onto the bisulfite converted Mus musculus genome (mm10) using Bismark v.0.19.0 (https://www.bioinformatics.babraham.ac.uk/projects/bismark/). The methylKit R package (v.1.6.0) (26) was used to analyze DNA methylation. Differentially methylated regions (DMRs) are presented in Fig. 2. The base-pair methylation profile was summarized over a tiling window (window size: 1,000 bp and step-size: 500 bp). Statistically significant DMRs were selected based on their Q-value and the percentage of methylation difference (CpGs, mean methylation difference in groups >10%; Q<0.01). The genomation R package (27) was used to annotate the genomic features of DMRs, which were classified into six types (intergenic, 5′UTR, 3′UTR, promoter, CDS and intron) based on University of California Santa Cruz Genome Browser (28).

DNA methylation and gene expression are known to be closely related. To identify the genes whose RNA expression levels were affected by DNA methylation (i.e. inverse correlation between RNA expression and DNA methylation) in asthmatic mice, the expression level was adjusted to between 0 and 1 (as methylation β values range between 0 and 1). Spearman's rank correlation coefficient was used to evaluate the correlation between CpG methylation and gene expression. To understand the biological mechanisms and pathways related to these genes in asthma, Gene Ontology (GO) enrichment analysis (29,30) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (31) mapping were conducted using DAVID (24).

Ingenuity Pathway Analysis (IPA; v.1.13, Qiagen, Inc.) software was used to predict putative sub-networks containing the four candidate genes, whose expression altered due to changes in DNA methylation in asthma. Briefly, the Molecule Activity Predictor tool in IPA was used to interrogate sub-networks using the methylation-expression candidate genes, then integrated the sub-network by overlaying the genes involved in asthma. To examine how each gene affected the sub-network, the fold change of gene expression and DNA methylation β values were included.

The data from the present study were deposited in the GEO under accession number GSE114587 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE114587; secure token for reviewer: itytgkuoptclnej) (18).

The DNA methylation of samples from the same animals were profiled in order to identify DMRs between the vehicle and the OVA-induced asthma samples (Fig. 2). After filtering the methylation profile, a total of 35,401 DMRs (base pair-level) were identified between the OVA-induced asthma and the saline-treated mice. Of these, 3,060 were located in promoter regions and 370 of the genes containing these DMRs demonstrated an inverse correlation between methylation and gene expression. Additionally, the number of hyper-methylated regions was twice as high as the number of hypo-methylated regions. Nevertheless, the proportion of functional regions, such as promoters and exons, was similar between the hyper- and hypo-methylated regions (hyper-methylated: 20.1%; hypo-methylated: 20.7%). The regional distribution of DMRs was assessed based on their distance from the nearest CpG island (CpGi); shores were defined as the regions flanking the CpGi (0–2 kb range). The majority of DMRs identified in mice with OVA-induced asthma were located in other regions (89% non-CpGi and non-shore); this was followed by shores (10%) and CpGi (2%). The basepair-level DMRs were combined with the differentially methylated genes to perform an integrated analysis of methylation and gene expression. The number of hyper-methylated genes was also higher than the number of hypo-methylated genes. These DMRs, which were verified as important genes, were used to infer the mechanisms underlying the development of asthma. DAVID was used to functionally annotate the genes (Tables I and II). KEGG pathway analysis identified that 1,235 genes were more highly expressed in the OVA-induced asthma samples, than in the vehicle samples, and were mainly involved in ‘chemokine signalling pathway’, ‘cytokine-cytokine receptor interaction’, ‘MAPK signalling pathway’, ‘neuroactive ligand receptor interaction’, ‘olfactory transduction’ and ‘oxidative phosphorylation’ pathways (Table I). Conversely, 170 genes were expressed at lower levels in the OVA-induced asthma samples than in the vehicle samples, and were enriched in the ‘cytokine-cytokine receptor interaction’, MAPK signalling pathway’, ‘neuroactive ligand receptor interaction’ and ‘olfactory transduction’ pathways (Table II).

To verify the genes modulated by DNA methylation in the development of asthma, the RNA-seq data were analyzed for correlations between DNA methylation and gene expression (Fig. 3). A total of 95 genes were upregulated, whilst 275 genes were downregulated in the OVA-induced asthma samples. DAVID was used to functionally annotate the genes. KEGG pathway analysis identified that 368 genes were upregulated or downregulated in the OVA-induced asthma samples, and were mainly involved in ‘acute myeloid leukemia’, ‘chemokine signalling pathway’, ‘focal adhesion’, ‘leukocyte transendothelial migration’, ‘vascular smooth muscle contraction’ and ‘neurotrophin signalling pathway’ (Table III).

An integrated network analysis of DNA methylation and gene expression was conducted to understand the interactions between the genes related to asthma development. The predicted pathway is presented in Fig. 4. There were four hub genes, consisting of three upregulated genes [forkhead box O1 (FOXO1), SP1 transcription factor (SP1) and amyloid β precursor protein (APP)] and one downregulated gene [RUNX family transcription factor 1 (RUNX1)], which demonstrated the association between DNA methylation and gene expression. These genes were closely interconnected nodes in the IPA module, and were functionally significant. Therefore, FOXO1, SP1, APP and RUNX1 were identified as key integrated network regulators in the development of asthma.