2.4. RNA sequencing and bioinformatic analyses
RNA integrity number (RIN) was measured using Agilent 2100 bioanalyzer (Palo Alto, CA, USA) and samples with RIN of 7.5 or above were sent to CnK genomics for sequencing. Libraries were prepared using the Illumina TruSeq Stranded mRNA sample preparation kit (San Diego, CA, USA), and a total of six libraries were constructed using an Illumina Nexseq 500 platform (San Diego, CA, USA). Raw reads were trimmed by filtering out adaptors with a minimum length of 75 bp using Trimmomatic v0.36 (http://www.usadellab.org/cms/index.php?page=trimmomatic) (Bolger et al., 2014). The quality score was evaluated using fastqc (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) both before and after trimming. Trimmed reads were aligned to apple reference genome GDDH13 v1.1 (http://www.iris.angers.inra.fr/gddh13) using HISAT2 software (Kim et al., 2015). Across all the libraries, gene counts were calculated for each predicted coding DNA sequence using FeatureCounts v1.5.2. Read counts in the range of 21,543,441 to 24,876,113 were obtained with Q30 ratio > 0.9. Reads per kilobase per million (RPKM) values were counted from BAM files. Differentially expressed genes (DEGs) with a false discovery rate (FDR) value of < .05 and |log2 fold change| > 1 were selected using EdgeR Bioconductor software (Empirical analysis of digital gene expression data in R) (Robinson et al., 2010). Functional enrichment analysis of DEGs was performed using InterPro2 (http://www.ebi.ac.uk/interpro) and SwissProt (www.ebi.ac.uk/swissprot). All sequencing data were deposited in National Center for Biotechnology Information Sequence Read Archive database bearing the BioProject ID PRJNA613278.