*4.4. mRNA Transcriptome Data Analysis*

The reads with adaptors were removed firstly. The reads in which unknown bases comprised more than 5% of the total and low quality reads (the percentage of the low quality bases of quality value ≤ 15 is more than 20% in a read) were also removed. The clean reads were aligned to *Phalaenopsis equestris* genome accessed from (https://www.ncbi.nlm.nih.gov/bioproject/192198) by HISAT2 (http://www.ccb.jhu.edu/software/hisat/index.shtml) [37].

Gene alignment rate was counted with Bowtie2 v2.2.5 (http://bowtie-bio.sourceforge.net/bowtie2/ index.shtml), gene and transcript expression levels were calculated with RSEM v1.2.12 (http://deweylab. biostat.wisc.edu/RSEM). The differentially expressed genes (DEGs) between spot and non-spot tissues were identified and filtered with the R package DESeq2 based on |log2 (foldchange)| ≥ 1 and Padj < 0.05 [38]. The heatmap displays of the Trimmed Mean of M-values (TMM) normalized against the Fragment Per Kilobase of transcripts per Million (FPKM) mapped reads were performed using the R package pheatmap (https://cran.r-project.org/src/contrib/Archive/pheatmap/). The DGEs were imported into Blast2GO software v2.5 [39] and in-house perl scripts for gene ontology (GO) term analysis, while the KEGG pathways were assigned to the assembled genes using the online KEGG Automatic Annotation Server (KAAS, http://www.genome.jp/kegg/), and then the phyper function in R software was used for enrichment analysis with *FDR* ≤ 0.01. To identify transcription factors, the assembled transcriptome were searched against the Plant Transcription Factor Database PlnTFDB (http://plntfdb.bio.uni-potsdam.de/v3.0/downloads.php) using hmmseach v3.0 (http://hmmer.org).
