RNA Extraction and RNA-Seq Analysis

Total RNA was extracted from 15 shoot apices after removing all visible leaves under a binocular for each of the three independent biological replicates using an RNeasy Plant Mini Kit (Qiagen) and treated with DNase (Ambion) to remove residual genomic DNA. Library for sequencing was prepared using a TruSeq Library Preparation Kit (Illumina) according to the manufacturer’s protocol. Sequencing was performed using the HiSeq3000 platform (Illumina) in 150-bp single reads. For each sample, ∼15,000,000 reads were generated. The software tool FastQC was used to assess quality control parameters (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). To estimate expression levels, the RNA-seq reads were mapped to the Arabidopsis TAIR10 (Lamesch et al., 2012) reference genome using the software TopHat2 under default settings (Kim et al., 2013), except that only a single alignment was permitted per read and the coverage-based junction search was disabled (settings: −g 1 –no-coverage-search). The program Samtools was used to sort and index BAM alignment files and to calculate BAM file statistics (Li et al., 2009). The software HTSeq was used to tabulate the number of reads mapping to each genomic feature, with counts tabulated only for genes that completely overlapped a given feature (Anders et al., 2015). We used the Wald test implemented in the program DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) to detect DEGs for pair-wised comparison. To visualize the expression levels of candidate genes, the expression level for each gene was calculated as transcripts per million (TPM).