Integrative analyses of targeted metabolome and transcriptome of Isatidis Radix autotetraploids highlighted key polyploidization-responsive regulators

Isatidis Radix, the root of Isatis indigotica Fort. (Chinese woad) can produce a variety of efficacious compound with medicinal properties. The tetraploid I. indigotica plants exhibit superior phenotypic traits, such as greater yield, higher bioactive compounds accumulation and enhanced stress tolerance. In this study, a comparative transcriptomic and metabolomic study on Isatidis Radix autotetraploid and its progenitor was performed. Through the targeted metabolic profiling, 283 metabolites were identified in Isatidis Radix, and 70 polyploidization-altered metabolites were obtained. Moreover, the production of lignans was significantly increased post polyploidization, which implied that polyploidization-modulated changes in lignan biosynthesis. Regarding the transcriptomic shift, 2065 differentially expressed genes (DEGs) were identified as being polyploidy-responsive genes, and the polyploidization-altered DEGs were enriched in phenylpropanoid biosynthesis and plant hormone signal transduction. The further integrative analysis of polyploidy-responsive metabolome and transcriptome showed that 1584 DEGs were highly correlated with the 70 polyploidization-altered metabolites, and the transcriptional factors TFs-lignans network highlighted 10 polyploidy-altered TFs and 17 fluctuated phenylpropanoid pathway compounds. These results collectively indicated that polyploidization contributed to the high content of active compounds in autotetraploid roots, and the gene–lignan pathway network analysis highlighted polyploidy–responsive key functional genes and regulators.


Introduction
Isatis indigotica Fort. (Chinese woad, 2n = 14) belongs to Isatideae tribe of the Brassicaceae family. The root of I. indigotica called Isatidis Radix, which can produce a variety of chemicals with medicinal properties, can be used in clinical treatment of regular seasonal influenza and plays an immune regulatory role in vitro and in vivo [1], while the leaves of I. indigotica called Isatidis Folium composed of isatin, tryptanthrin, indirubin and so on [2,3]. Several categories of metabolites including alkaloids, phenylpropanol, organic acids and polysaccharides identified from Isatidis Radix, were demonstrated to achieve the antiviral and antioxidant effects [4][5][6]. Therefore, increasing the abundance of the active compounds is critical for improving the quality of Isatidis Radix [7].
Polyploidy is widespread in plants, and nearly 70 % of angiosperms are polyploids including many important crops [8]. Polyploidization, also known as whole genome duplication (WGD), plays a pivotal role in promoting the evolution of plant morphological, physiological and reproductive diversity [8][9][10][11][12]. Compared with their diploid progenitors, polyploid plants often exhibit superior phenotypic traits, such as stronger tolerance, higher content of active compounds, and enlarged organs together with increased vigor [8,[13][14][15][16]. The most conspicuous features of polyploidy are the increased cell size, slowed cell division and tissue development, and increased organ size at maturity, which is referred to as the 'gigas effect' [11,15]. The tetraploid I. indigotica accumulate more lignans than diploid, including lariciresinol and its derivatives, which present effective antiviral ingredients of I. indigotica [17]. The giant organs and enhanced concentrations of secondary metabolites realized by autopolyploidy are attractive for breeding the respective medicinal and agricultural plants.
However, there has been no report on the metabolomic and transcriptomic changes post polyploidization of Isatidis Radix until now. In the past several years, the research in the field of polyploidy is mainly focused on the transcriptional level, using RNAseq-based transcriptomic analysis to reveal the relationship between polyploidization and gene expression [16]. At present, only two reports of the Chinese Woad leaf (Isatidis Folium) transcriptomic changes induced by autotetraploidization were available [18,19]. However, the root (Isatidis Radix) differed from leaf (Isatidis Folium) whether in biological function or in medicinal usage. And there are three monographs of I. indigotica included in the Chinese Pharmacopoeia, namely Isatidis Folium, Indigo Naturalis and Isatidis Radix [20], so these three kinds of Chinese herbal medicine preparations are somewhat different. Moreover, the metabolome is closer to phenotype than transcriptome.
Given that the metabolic activity was altered by the fluctuated gene expression, which led to the change of the concentration of secondary metabolites, we carried out a comparative transcriptomic and metabolomic study on Isatidis Radix autotetraploid and its progenitor. Through the integrative analysis of Isatidis Radix transcriptome and metabolome, the differentially expressed genes affecting the metabolic pathway of active components such as lignan were identified. As a result, the gene-lignan pathway network analysis highlighted polyploidy-responsive key functional genes and regulators.

Plant materials and sampling
Appropriate permissions for collection and use of seed of Isatis indigotica Fort. (2n = 2x = 14) was obtained from Jiangsu Germplasm Repository Center. I. indigotica (2n = 2x = 14) used as diploid donor. Autotetraploid I. indigotica was artificially synthesized by colchicinemediated polyploidy induction in vitro as described previously [16,19]. Briefly, adventitious buds induced from diploid planet were subjected to 0.20 % colchicine treatment for 12 h, and transferred to MS medium without colchicine for 2 weeks. Then, the synthesized autotetraploid plantlets were transferred to 1/2 MS medium for rooting. The root tips (0.5-1 cm long) were excised and pretreated with 2 mmol·L − 1 8-hydroxyquinoline solution for 4 h, and fixed with Carnoy's solution at 4°C for 24 h. Samples were then hydrolyzed using 1 mol·L − 1 HCl at 60°C for 10 min. The hydrolyzed root tips were soaked in a drop of Carbol fuchsin for 10 min and squashed on the microscopic slide to observe the metaphase chromosomes. Finally, seedlings with roots were transplanted into nutritional soil. The diploid and autotetraploid I. indigotica seedlings were planted in the experimental fields in our campus for 1-year. Then, their roots were sampled for the subsequent transcriptomic and metabolomic analysis with three repeats (Fig.S1). The sampled fresh roots of I. indigotica were frozen with liquid nitrogen, transported and stored at -80 ℃.

Targeted metabolomic analysis of Isatidis Radix metabolites
The Isatidis Radix samples were freeze-dried and ground into fine powder for metabonomic analysis. The widely targeted metabolic profile and quantitative detection of metabolites were performed by MetWare Biotechnology Co.,Ltd (Wuhan, China) (www.metware.cn). The quantification of metabolites was carried out using a predetermined multi-reaction monitoring method [21].
The elemental composition and mass fragmentation were compared to those registered inaccessible databases of NIST as well as the standards in a database compiled by MetWare Biotechnology Co.,Ltd [22].

RNAseq libraries preparation and sequencing
Total RNA for RNAseq was extracted from seedling roots and about 1 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) and index codes were added to attribute sequences to each sample. For high-throughput sequencing, the library preparations were sequenced on an Illumina Hiseq X Ten platform and 150 bp pairedend reads were generated [19]. After the adaptor and low-quality sequences were trimmed, a total of 38.71 Gb clean data from 6 cDNA libraries were retained (Table S1).
The RNAseq reads for each sample were mapped to the reference genome using HISAT2, and the output SAM files were sorted and converted to BAM files using SAMtools (version 0. 1.19). Then the sorted alignments were assembled into transcripts and the expression levels of all genes and transcripts were estimated using StringTie.

Analysis of the differentially expressed genes (DEGs)
The expression values were represented by fragments per kilobase transcript per million reads mapped (FPKM), and the differential expression analysis of genes and transcripts across two conditions was performed using the Cuffdiff utility. Foldchange ≥ 2 and FDR ≤ 0.05 was set as the threshold to determine the DEGs between the compared samples. The KEGG pathway enrichment analysis of DEGs was conducted by Path_finder software with Q-value ≤ 0.1 [19].

Integrative targeted metabolomic and transcriptomic profiling analysis
The data of metabolites profiling were normalized and exported to Simca-P software (12.0, http://www. umetrics.com/simca) employing partial least-squares discriminant analysis (PLS-DA) model. The differentially expressed metabolites were discriminated according to a threshold of variable importance in the projection (VIP) values (VIP > 1) after PLS-DA processing using the previously published protocols [17].
The correlation among lignan biosynthetic genes and lignans was constructed using the Pearson correlation coefficient according to the co-occurrence principle. The correlation network was generated using Cytoscape [17].

qRT-PCR Analysis
In order to verify the differentially expressed genes (DEGs), the total RNA of 3 individuals of each genotype which were used for the aforementioned metabolite profiling was extracted by Total RNA Kit II (Qiagen). Then, DNaseI treatment, RNA concentration measurement and cDNA synthesis were carried out. According to the RNA-seq data, Primer5 software was used to design primer pairs for randomly selected DEGs (primers were listed in the supplementary Table S2). The housekeeping gene UBIQUITIN1 was used as the reference gene to calculate the relative expression of genes using the comparative Ct method [11].

Metabolomic alterations in Isatidis Radix following autopolyploidization
In order to assess the impact of polyploidization on the metabolomic shifts, the extracts from Isatidis Radix autotetraploids and their diploid parents were subjected to the targeted metabolic profiling by UPLC-TOF/MS. Totally, 283 annotated metabolites were identified in Isatidis Radix, the roots of I. indigotica.
It is noteworthy that the biosynthesis of Lariciresinol glucosides was enhanced in autopolyploid ( Fig. 3, Table  S3). Regarding the flavonoid, the accumulation of Phlorizin was enhanced in autotetraploid in comparison to diploid (Table S3). These compounds can be regarded as the indictor components of polyploidization based on metabolomics data.

Polyploidization-modulated changes in lignan biosynthesis
Metabolic analysis revealed that the production of lignans was significantly increased post polyploidization, since Coniferyl alcohol and Lariciresinol glucosides accumulated more in autotetraploid than in diploid. Moreover, Coniferyl aldehyde and the subsequent Coniferyl alcohol, the critical precursor for lignan biosynthesis, were all enhanced in autotetraploid.
However, ferulic acid, the precursor of Coniferyl alcohol as well as Sinapyl alcohol, was less accumulated in autotetraploid than in diploid. Pinoresinol and its deriva- displayed less accumulation in autotetraploid than in diploid, as was the case for Secoisolariciresinol glucoside or Lariciresinol-4,4'-di-O-glucoside (also named clemastanin B) (Table S3), which implied that polyploidization-modulated changes in lignan biosynthesis.

Polyploidization-responsive genes in Isatidis Radix
RNAseq-based transcriptomic profiling was performed to investigate the polyploidization imposed profound impacts on gene expression and the subsequent metabolic pathways in Isatidis Radix. Using a stringent cutoff (Foldchange > 2 and FDR ≤ 0.05), a total of 2065 differentially expressed genes (DEGs) were identified as being polyploidization-responsive genes, of which 1251 were polyploidy-induced and 814 were polyploidy-repressed in I. indigotica seedlings roots (Table S4). To further evaluate the functions and the biological pathways represented by the DEGs, we compared these genes with that included in the KEGG database [25]. The annotation and classification of root DEGs indicated that the polyploidization-altered root genes were enriched in phenylpropanoid biosynthesis and plant hormone signal transduction (Fig. 4).
To gain insights into the functionality of the 2065 DEGs that are likely to be associated with the process of polyploidization, all of these polyploidy-responsive transcripts were functionally grouped (Table 1). Among DEGs mainly involved in stress response, L-ascorbate peroxidase S, trehalose-phosphate synthase TPS7 and Senescence/dehydration-associated protein were upregulated. Regarding the upregulated DEGs involved in growth and development, Light-regulated protein 1, GIGANTEA, Glycine-rich protein 3 and two HAIKU members were of particular interest. In the kinase and signaling category, three kinase (Receptor-like protein kinase FERONIA, Wall-associated receptor kinase WAK14 and Mitogen-activated protein kinase MPK19) and Gibberellin receptor GID1C were upregulated by polyploidization. Regarding the transporters, Glutathione S-transferase GSTZ1, Sulfate transporter AST12 and two ABC transporter C family members (ABCG14/36) were upregulated by polyploidization (Table 1).

Systematic transcriptomic and metabolomic shift post polyploidization
To integrate the analysis of polyploidy-responsive metabolome and transcriptome, a canonical correlation analysis using Pearson's correlation coefficient was performed to display the dynamic variation over the polyploidization course. This integrative analysis showed that 1584 DEGs were highly correlated with the 70 polyploidization-altered metabolites, with |PCC|>0.917 (Table S5).
A TF-metabolite correlation network was built that consisted of 15 polyploidy-altered TFs and 67 fluctuated compounds to characterize TFs involved in polyploidyinduced alteration in roots metabolome and transcriptome (Fig. 5, Table S6).

Integrated metabolomic and transcriptomic analysis of lignan metabolism modulated by polyploidization
To have a systematic view on the polyploidy-responsive variation of lignan biosynthesis, the transcripts involved in the general phenylpropanoid pathway, lignan biosynthesis and the corresponding metabolites were subjected to construct lignan biosynthesis pathway.
In the polyploidy-altered TFs-lignans network, there were 10 polyploidy-altered TFs and 17 fluctuated phenylpropanoid pathway compounds, which indicated the transcriptomic and metabolic shifts in lignan metabolism as a result of polyploidization-mediated transcriptional regulation (Fig. 7, Table S7). Among the 10 polyploidy-altered transcriptional factors, polyploidy-inhibited bHLH44 and polyploidyinduced bHLH129 were highly correlated with all the fluctuated phenylpropanoid pathway compounds (7 up and 10 down) except Syringaresinol, but with the reverse trend.

Discussion
Polyploidization contributed to the 'gigas effect' and high content of active compounds in I. indigotica autotetraploid roots Polyploidization, also known as whole-genome duplication (WGD), results in the "gigas effect" that includes increased cell size, enlarged vegetative or reproductive organs and prolonged vegetative growth [11,15,26]. Compared to their diploid progenitors, the autotetraploid I. indigotica plants exhibit bigger robustness and larger leaves with deeper color, which was in accordance with the "gigas effect".
Among DEGs mainly involved in stress response, Lascorbate peroxidase S, trehalose-phosphate synthase (TPS7) and senescence/dehydration-associated protein were upregulated by polyploidization (Table 1). In Arabidopsis, TPS1 catalyzes the synthesis of the sucrosesignaling metabolite trehalose 6-phosphate which acts as a potent regulator of post-embryonic growth and development [27]. Moreover, rice OsTPS1 may improve the abiotic stress tolerance by increasing the accumulation of trehalose and proline, and modulating the expression of stress-related genes [28]. Regarding the senescence/ dehydration-associated protein, Arabidopsis ERD7 and its homologs play essential roles in plant stress responses and development and are associated with modification of membrane lipid composition [29]. Therefore, the role of these polyploidy-enhanced genes in 'gigas effect' and stress tolerance of autotetraploid needs further establishment.
Regarding the kinase genes upregulated by polyploidization, the receptor-like protein kinase FERONIA was of particular interest (Table 1). In Arabidopsis, the couple of extracellular peptide RAPID ALKALINIZATION FACTOR1 (RALF1) and FERONIA (FER) acted as a central hub between the cell surface and downstream signaling events, and the RALF-FER pathway functioned as an essential regulator of plant stress responses [30]. Furthermore, the RALF1-FER-GRP7 module provided a paradigm for regulatory mechanisms of RNA splicing to regulate plant fitness and flowering time [30,31]. Glycine-rich proteins (GRPs) were demonstrated to participate in cold stress responses, plant defense, cell elongation and fertility. Moreover, rice glycine-rich protein OsDG2 plays important roles in chloroplast development during early seedling stage [32]. In this study, glycine-rich protein (Iin13953) was one of the polyploidy-upregulated DEGs which were involved in growth and development (Table 1). Hence, it is interesting to investigate the contribution of polyploidy-induced FER together with glycine-rich protein Iin13953 to the modulation of cell growth and stress responses.
Arabidopsis root hair defective 6-like 4 (RSL4), a bHLH transcription factor, triggers the expression of hundreds of root hair genes which promote ectopic root hair growth, and the autocrine regulation of root hair size by the RALF-FERONIA-RSL4 signaling pathway has been revealed [33]. In this study, polyploidy-induced bHLH129 and bHLH63 together with FER were identified (Tables 1 and 2), but whether they acted as central hubs orchestrating complex intracellular and extracellular signals required further elucidation.
One of the ideal expectations for the medicinal autopolyploid was that the organ giantism was accompanied by the higher content of some chemical compositions, especially the active compounds [8,[13][14][15][16]19]. I. indigotica, like I. tinctoria of the Brassicaceae family, represents a valuable source of bioactive compounds such as alkaloids, phenolic compounds, phenylpropanoids and terpenoids [3,19]. In this study, the content of active compounds in the roots of autotetraploid I. indigotica was higher than that in diploid roots, and some new compounds including Phlorizin, Tannins and Solatuberenol A were also isolated in Isatidis Radix (Table S3). Moreover, Isatidis Radix autotetraploid accumulated more Lariciresinol glucosides than the diploid counterparts, consistent with previous report [17]. Known indolic alkaloids called Indirubin and indicant (iso) which were reported in the dried I. tinctoria leaves [2,3] were also identified to be polyploidy-upregulated metabolite in this study (Table S3). These collectively implied the potentiality of enhancing active compounds accumulation through polyploidization.
The penetration-resistance gene PEN3/ABCG36/PDR8 and PDR12 function redundantly to mediate the secretion of camalexin, and they have multiple functions in Arabidopsis immunity via transport of distinct Trp metabolic products [34]. PEN2 encodes a myrosinase that catalyzes the degradation of indole glucosinolates, and the catalyzed products of PEN2 are postulated to be transported to the apoplast by PEN3. Moreover, the indole compound 4-methoxyindole-3-methanol that is a substrate for PEN3 stimulates bacterial flg22-induced callose deposition [35]. In this study, one myrosinase (Iin26136) and two myrosinase-binding protein genes were induced by polyploidization (Table S4). Nonetheless, whether the coordinated function mechanism of PEN2-PEN3 play its role in synthesis and export of the active metabolites (including but not limited to indole compounds) in Isatidis Radix merits further investigation.
Gene-lignan pathway network analysis highlighted polyploidy-responsive key functional genes and regulators Phenylpropanoid is the major group of secondary metabolites, which metabolism generate diverse metabolites including lignans and flavonoids, and lignans are identified to be the pharmacologically active compounds [7], therefore the correlations between the identified DEGs and phenylpropanoid pathway compounds were inferred based on the co-occurrence principle between the transcript and metabolite levels. Several key metabolites involved in general phenylpropanoid pathway (e.g. coniferyl aldehyde and coniferyl alcohol) and lignan compound (e.g. Lariciresinol glucoside) were markedly increased post polyploidization. Moreover, various catalytic genes (e.g. C4H, 4CL, COMT and F5H) showed similar up-regulated patterns in correspondence with the increased metabolites (Fig. 6), suggesting that transcriptomic and metabolomic profile of lignan biosynthesis pathway was modulated by polyploidization.
Transcriptional factors were predicted to act as key regulators of lignan synthesis in I. indigotica. Among the 10 polyploidy-altered TFs, polyploidy-inhibited bHLH44 and polyploidy-induced bHLH129 were highly correlated with all the fluctuated phenylpropanoid pathway compounds (7 up and 10 down) except Syringaresinol, but with the reverse trend (Fig. 7). It was also intriguing to note that the polyploidy-induced NAC54 was positively correlated with polyploidization-enhanced Sinapyl alcohol (Fig. 7). Further studies of the regulatory mechanism of polyploidy-induced bHLH129 and NAC54 may provide fruitful means to reveal the beneath mechanism for polyploidy vigor and lignan biosynthesis of I. indigotica.
It was reported that IiWRKY34 significantly contributed to the polyploidy vigor of I. indigotica, and IiWRKY34 positively contributed to the yield, lignan biosynthesis and stress tolerance in I indigotica hairy roots, however, this key regulator was not identified here using the genuine roots namely Isatidis Radix. One possible explanation is that the expression pattern of genes in the induced hairy roots of tetraploid I. indigotica greatly differed from that in its original root [17].
Not surprisingly, different TFs may play distinct roles in lignan biosynthesis and allow the autotetraploid roots to prioritize toward a more efficient lignan biosynthesis. Therefore, whether these highlighted TFs regulate DEGs for lignan biosynthesis is the most important issue to elucidate the genuine regulators of lignan biosynthesis in Isatidis Radix.