Comparative transcriptomic of Stevia rebaudiana provides insight into rebaudioside D and rebaudioside M biosynthesis

Rebaudioside D (Reb D) and rebaudioside M (Reb M) are commercially important low/no-calorie natural sweeteners. However, they are present in a minor proportion of all steviol glycosides (SGs) in Stevia rebaudiana Bertoni (S. rebaudiana). Strain-dependent deviation in Reb D and Reb M biosynthesis is one key breach for breeding of S. rebaudiana, which has not been studied at the transcriptional level. Herein, five different S. rebaudiana varieties with distinct SGs contents, one cultivar having high stevioside content (HST), one cultivar having high Reb A content (HRA) and three cultivars having high Reb D and Reb M content (HDM₁, HDM₂, HDM₃), were selected for RNA-seq analysis. In total, 131,655 de novo assembled unigenes were found in the RNA-seq data. According to Reb D and Reb M content divergence of S. rebaudiana accessions, 2186 differentially expressed genes (DEGs) were selected as potential genes related to Reb D and Reb M biosynthesis. Weighted Gene Co-expression Network Analysis (WGCNA) was used to explore the genes associated with the Reb D and Reb M biosynthesis. The unigenes from the positively associated turquoise module formed a layered co-expression network. There are 7 UDP-dependent glycosyltransferases (UGT) and 76 transcription factors (TFs) distributing at different regions which represented varying coherence of Reb D and Reb M biosynthesis. Particularly, two TFs having a strong correlation with two UGTs in the network were also discovered. The present study provided a comprehensive insight into networks for regulation of Reb D and Reb M contents in S. rebaudiana.