Transcriptome Analysis of Monoterpenoid Biosynthesis

For RNA library construction, poly-T oligo-attached magnetic beads were used to purify the RNA from total RNA (approximately 5μg) with twice purification. Next, the mRNA was fragmented into small pieces via divalent cations at high temperature. Then the final cDNA library was established by reverse-transcribed of the cleavage RNA pieces in accordance with the guidance for the mRNASeq Sample Preparation Kit RS-122-2103 (Illumina, San Diego, USA). The average insert size for the libraries was 300 bp (±50 bp). The prepared libraries were then sequenced on an IlluminaHiseq 2000 (LC Sciences, San Diego, USA) platform according to the vendor’s recommended protocol, and 150bp paired-end reads were generated. After filtering the primers, adapters and reads with low-quality and >5% bases using Cutadapt software [38 ], de novo assembling of the transcriptome was performed by Trinity 2.4.0 [39 ] with default parameters on the foundation of an overall 60.61GB RNA-seq data. Benchmarking Universal Single-Copy Orthologs (BUSCO) [40 (link)] were used to assess transcriptome assembly and annotation completeness. Unigenes were then queried against the Swiss-Prot, Non-redundant (Nr), Protein family (Pfam), Kyoto Encyclopedia of Genes and Genomes (KEGG), eukaryotic Orthologous Groups (KOG), Gene Ontology (GO) public databases by BLASTx with an E-value <10⁻⁵ to obtain annotation. RPKM (Reads Per Kilobase per Million mapped reads) were employed to evaluate genetic expressing levels. Differential expression analysis were performed using edgeR [41 (link)] to obtain p value, and Benjamini-Hochberg (BH) algorithm was used to adjust p value for controlling false discovery rate (FDR). To define differentially expressed genes, we used fold change (FC) ≥1.5 and FDR ≤0.05 as cutoff. Subsequently, GO and KEGG enrichment analyses and the heatmaps of differential expressed genes (DEGs) among four tissues were performed using the OmicStudio tools at https://www.omicstudio.cn/tool.
GeneMANIA app in cytoscape version 3.7 was employed to discover genes modulated by identified transcription factors (TFs) [42 (link)]. Genes screened and annotated as monoterpenoid biosynthesis related genes (TPs) along with TFs activities were adopted for the purpose of drawing a co-expressing network. Arabidopsis co-expression network was employed as a reference to query [43 (link)], genes with the greatest normalized weight (high associated genes) were forecasted to be related genes correlated with a specific role. The normalized weight was calculated by default weighting method of GeneMANIA app.