Virus and cell lines
Monolayers of human bronchial epithelial (16HBE) cells were purchased from Jennino Biological Technology (Guangzhou, China). The cells were seeded at a density of 5 × 105 cells per well in 6-well sterile plastic culture plates and grown in a base media of Roswell Park Memorial Institute (RPMI)-1640 containing 10% (vol/vol) fetal bovine serum (FBS), 100 units/mL penicillin, 100 μg/mL streptomycin and 2 mM L-glutamine at 37°C in a 95% air and 5% carbon dioxide (CO2) incubator. For transcriptome study, 16HBE cells were infected at a multiplicity of infection (MOI) of 0.01 with enterovirus 71 (EV-A71; subgenotype C4, GenBank: EU812515.1) or coxsackievirus A16 (CV-A16; subgenotype B, GenBank: JN590244.1), which were isolated from an epidemic in Fuyang, China in 2008 and an HFMD patient in Guangxi, China in 2010, respectively. Both viruses were incubated for one hour to attain virus attachment, and RPMI-1640 containing 1% FBS was then added. The cells were incubated for a further 6 h, 12 h and 24 h for virus propagation. At the end of each incubation point, cells were collected using RPMI-1640 and centrifuged twice in phosphate-buffered saline (PBS) at 4°C for 5 min at 3000 rpm/min. Finally, the harvested cell pellets were snap-frozen in liquid nitrogen and stored until RNA extraction. For the control sample, the same process was carried out with the exception that 2 mL of sterile PBS was used to replace the virus.
RNA extraction and RNA sequencing (RNA-seq) library construction
Total RNA was extracted using an RNeasy® Mini Kit (Qiagen, USA) as per the methods recommended by the manufacturer. The extracted RNA was cleared of contaminating genomic DNA by RNase-free DNase I (Takara Bio, Japan) treatment. The concentration of RNA was measured with a NanoDrop 2000 (Thermo Scientific, USA), and the RNA quality was determined by an Ultrospec 3000 Pro UV/Visible spectrophotometer (GE Healthcare, UK), where samples with an absorbance ratio (A260/A280) of 1.8 to 2.0 were considered for RNA integrity analysis with an Agilent® 2100 Bioanalyzer. Samples with an RNA integrity number (RIN) greater than 7 were sent for RNA-seq library construction. Briefly, mRNA was enriched by Oligo (dT) beads, and poly(A)-containing mRNA was then purified using Dynabeads (Life Technologies, USA) and further fragmented into smaller pieces with fragmentation buffer using RNase III and an Ion adaptor. Next, the RNA fragments were reverse-transcribed and amplified to form first-stranded complementary DNA (cDNA) with random primers, followed by second-strand cDNA synthesis. The second-strand cDNA was further purified, adenylated at the 3’ ends, and ligated with sequencing adaptor. These fragments (250~300 bp in size) were subjected to PCR amplification with Phusion High-Fidelity DNA polymerase (NEB, Beijing, China), and the products were sequenced on an Illumina HiSeq™ 2000 platform (Illumina, USA).
Computational analysis of RNA-seq data
Transcriptome assembly
The raw RNA-seq paired-end reads were filtered to remove the “dirty” reads, i.e., those containing sequencing adapters, reads with 10% > Q < 20% bases, and reads of low quality (reads with ambiguous bases “N”) using the Fast QC package (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The obtained “clean” reads were mapped against the human reference genome GRCh38 with HISAT2 software. The fragments per kilobase of transcript sequence per million base pair (FPKM) values of each gene were determined by the length of the gene and read count mapped to the gene with Cufflinks (version 1.3.0)[15].
Principal component analysis (PCA)
PCA was employed for quality control, to identify problems with the experimental design, to find mislabeled samples and to visualize variations between the expression analysis samples by using the data from clean reads.
Identification of differentially expressed genes
For gene expression analysis, the expression level of each gene was calculated by quantifying the number of reads and was further normalized by a variation of the FPKM method. To identify differentially expressed genes between the two samples, Cufflinks was applied to calculate the T-statistic and the p-values for each gene. We calculated the expression ratios of 6 h/0 h, 12 h/0 h and 24 h/0 h as fold changes. All P values were adjusted by Benjamini and Hochberg’s approach to control the false discovery rate (FDR). Genes with a 2-fold change and an adjusted P value < 0.05 were regarded to be differentially expressed.
Unsupervised hierarchical clustering analysis of the common differentially expressed genes
To find commonalities between the EV-A71- and CV-A16-infected 16HBE cells, we constructed a Venn diagram using the obtained common differentially expressed genes in each group. The common differentially expressed genes were subject to unsupervised hierarchical clustering using Euclidean distance and average linkage to measure the cluster similarity/dissimilarity. A trend analysis was further utilized to identify commonalities between EV-A71 and CV-A16 infection. The common differentially expressed genes in the EV-A71 and CV-A16 infections that had the same trend were identified as genes with similar changes.
Functional group analysis
The Database for Annotation, Visualization and Integrated Discovery (http://david.abcc.ncifcrf.gov/), which utilizes Gene Ontology (GO) to check the biological process, molecular function and cellular components of common differentially expressed genes, was applied in the current study. In addition, the Biocarta and Reactome database (http://www.genome.jp/kegg/), which uses the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, was applied to identify pathways of common differentially expressed genes in this study. The FDR-corrected P-value threshold was set at 0.05, which denotes the significance of GO term enrichment and pathway correlation.
Coexpression network construction
The coexpression network construction was based on the GenMANIA algorithm by using the common differentially expressed genes. The construction of coexpression networks is conducive to finding potential mechanisms associated with differentially expressed genes.
Validation of differentially expressed genes by quantitative RT-PCR (qRT-PCR) analysis
To validate the transcriptome results, we designed 10 pairs of primers using Primer Premier 5.0 to perform qRT-PCR analysis targeting 7 upregulated (BLM, GBP3, HDAC9, NNMT, PNISR, RNF6 and TAF1) and 3 downregulated genes (ANO1, IRX2 and PCDH7). Briefly, reverse transcription was first performed using total RNA (1 μg) isolated as described above with the Prime Script® RT reagent kit (Takara, Japan) according to the manufacturer’s protocol. Next, qRT-PCR was carried out using 1 μl of cDNA template, 0.2 μl of gene-specific primers (Supplementary Table S1), 5 μl of SYBR® Premix Ex Taq™ (2×), 0.2 μl of ROX Reference Dye II (50×), and water up to 10 μl. The qRT-PCR reaction program was as follows: predenaturation at 95 ℃ for 30 s, followed by 45 cycles of denaturation at 95 ℃ for 10 s, annealing at 55 ℃ for 15 s, and extension at 72 ℃ for 30 s on a 7500 Fast Realtime PCR system (Applied Biosystems, USA). Relative expression levels of the chosen genes were determined by normalizing the data for the target transcripts against the glyceraldehyde-3-phosphate dehydrogenase (GAPDH) transcript data according to the 2−△△Ct method. Three independent biological replicates for each sample and three technical replicates for each biological replicate were analyzed.