Study subjects and sample collection
Five patients with AERD and another 5 patients with ATA were enrolled at Ajou University Hospital in Suwon, South Korea. Asthma was diagnosed according to the Global Initiative for Asthma guideline (GINA) 2021 by the allergy specialists. AERD was defined by a recurrent clinical history (exacerbation of upper or lower respiratory reactions after ingestion of aspirin/NSAIDs) and/or a positive response to the lysine-aspirin bronchial provocation test (Lys-ASA BPT). The Lys-ASA BPT was performed with increasing doses of Lys-ASA solution up to 300 mg/mL using the method previously reported.9 The positive result of the Lys-ASA BPT was defined when FEV1% was decreased by more than 20% after the challenge. Asthmatics with negative results to the Lys-ASA BPT or denied any changes in upper or lower respiratory tract symptoms after ingestion of aspirin/NSAIDs were defined as ATA. Severe asthma was defined as the international ERS/ATS guidelines.10 Exclusion criteria for enrollment were as follows: 1) asthmatics who had been treated with asthma type 2 biologics, including omalizumab, mepolizumab, reslizumab, benralizumab, and dupilumab, within 130 days of enrollment; 2) current smokers or ex-smokers who quit smoking within 30 days of enrollment; and 3) asthmatics who used any of intranasal corticosteroids, intranasal antihistamines, and intranasal anticholinergics within 7 days of enrollment.
Nasal scraping was performed with a pencil-shaped disposable nasal curette (Rhino-probe®,Arlington Scientific, Inc., Springville, UT, USA) for 3 times from the middle portion of the each side of the inferior turbinate after 2 nasal lavages (10-mL saline each) to remove overlying mucus. Peripheral eosinophil counts, sputum eosinophil counts, total immunoglobulin E (IgE), fractional exhaled nitric oxide (FeNO) levels, pulmonary function test results, and atopy status were evaluated on regular asthma medication including inhaled corticosteroids with long-acting beta2-agonists at enrollment. All subjects gave written informed consent at the time of enrolment, and the study was approved by the Institutional Review Board of Ajou University Hospital (AJIRB-BMR-SUR-17-182).
RNA sequencing and data analysis
Total RNA was also extracted from nasal mucosa tissues using a Nucleo spin RNA kit (MACHREY-NAGEL, Dueren, Germany). The quantity and quality of total RNA were evaluated using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and a RNA quality indicator (RIN). Library preparation was conducted with an Illumina TruSeq Stranded mRNA sample preparation kit (Illumina, San Diego, CA, USA) according to the manufacturer′s instructions. Indexing adapters were ligated to the ends of the cDNA fragments using ligation mix reagent at 30°C for 10 min. After twice washing with sample purification beads, PCR was performed to enrich those cDNA fragments that have adapter molecules on both ends. The quality and band size of the libraries were assessed using an Agilent 2100 Bioanalyzer. Libraries were quantified by qPCR using the CFX96 Real Time System (Biorad, Hercules, CA, USA). Sequencing of the prepared library was conducted on the Hiseq2500 platform (Illumina) with 100 bp paired-end reads.
Reads were trimmed to remove adapters and low-quality reads (per-base quality < 20) and thereby improve paired-end mapping. High-quality sequence reads were mapped to the Human genome (hg19) using STAR,11 and gene expression levels were quantified with the RSEM (v. 2.12.0).12 Differentially expressed genes (DEGs) between the 2 groups were evaluated by the edgeR package (v. 3.0.8),13 which is based on negative binomial models for RNA-seq count data. DEGs were screened with a cutoff threshold of log2 (FC) ≥ |1| and P value < 0.05.
Methylation assay and data analysis
Genomic DNA was also extracted from nasal mucosal tissues in the same sample used for RNA-seq using Illumina Infinium Methylation EPIC Bead Chip kits (Illumina, Inc.), and arrays were scanned by the Standard Illumina procedures using an Illumina iScan scanner. Each methylation data point is represented by fluorescent signals from the M (methylated) and U (unmethylated) alleles. Background intensity computed from a set of negative controls was subtracted from each analytical data point. The ratio of fluorescent signals was then computed from the 2 alleles ß = (max (M, 0))/(|U| + |M| + 100). The ß-value reflects the methylation level of each CpG site. A ß-value of 0-1 was reported signifying percent methylation, from 0% to 100%, for each CpG site.
After methylation assay, array data export processing and analysis were performed using Illumina GenomeStudio v2011.1 (Methylatioin Module v1.9.0) and R 3.0.2 (http://www.r-project.org). Differentially methylated probes (DMP) and regions were detected through the in house python scripts using the numpy and scipy package14 and ChAMP15 with a cutoff threshold of P value < 0.001.
Correlation analysis between RNA expression and DNA methylation
Correlation analysis was performed to identify genes regulated by DNA methylation. To perform analysis with as many combinations between DEG and DMR as possible, the candidate gene for DEG and DMR were selected based on P value rather than the adjusted P value. To identify differentially correlated genes (DCG) which show an inverse relationship between expression and methylation, Spearman rank correlations were used to assess the relationship between CpG methylation and gene expression. R-package was used for statistical analysis
Pathway and network analysis
Ingenuity pathway analysis (IPA, Redwood City, CA, USA) was used to understand the functional characteristics of DEGs of AERD.16 IPA canonical pathway analysis was performed with DEGs to identify biological pathways that significantly affect mRNA expression in AERD. The gene symbols of probe, which significantly up-regulated or down-regulated in AERD (log2 (FC) ≥ |1| and P value < 0.05) compared to ATA were selected as significantly differentiated genes and used for IPA analysis. The significance of the association between the data set and the canonical pathways was measured as the number of molecules in each pathway that meets cutoff criteria (false discovery rate (FDR)-adjusted P values < 0.05).
Network analysis was performed using asthma-related DEGs to understand and predict interactions among these genes. The biological functions that were significantly associated with the genes in the core networks were identified by functional analysis based on the Ingenuity’s knowledge base. As our study is focused on asthma population, asthma-related terms were selected and used in further analysis. We determined asthma-related terms such as refractory asthma, severe asthma, uncontrolled asthma, eosinophilia of lung, pulmonary eosinophilia, T helper 2 (Th2)-mediated asthma, and aspirin-induced asthma. The DEGs were mapped onto the networks to explore the dynamic changes during the asthma process.
Isolation of total RNA from the cells and qRT-PCR
Results from RNA and methylation sequencing were verified by RT-qPCRin second cohort of AERD (n=12) and ATA (n=12). Total RNA was isolated from peripheral blood mononuclear cells from AERD and ATA patients using a Pure-Link RNA mini kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The RNA was reverse transcribed into cDNA with the ReverTraAce qPCR RT Kit (TOOBO, Osaka, Japan). The following primers (Bioneer, Daejeon, South Korea) were used: human STX-2 (forward: 5′-CGA GCT AGA AGA GAT GCT GG-3′; reverse: 5′-AAC TCT CGG ATG CTG GTC TC-3′), human SMPD3 (forward: 5′-TCA TGG ACG TGG CCT ATC AC -3′; reverse: 5′-TGC ACC TTG AGA AAC AGA GC -3′), human retinoic acid receptor–related orphan receptor γt (RORγt) (forward: 5′-TCT CAA AGC AGG AGC AAT GG-3′; reverse: 5′-ATG CCA CCG TAT TTG CCT TC -3′), and human FcER1A (forward: 5′-GCA GCT GGA CTA TGA GTC TG -3′; reverse: 5′-GAA TCA CCA CCA ACA ATG GG -3′). Real-time PCR was conducted with THUNDERBIRD SYBR Green Mix (TOOBO) on an ABI PRISM 7300 Sequence Detection System (Applied Biosystems, Carlsbad, CA, USA). Relative target gene expression was calculated by normalization to GAPDH with the 2−△△Ct method, and the results are presented as fold changes relative to the controls.