BACKGROUND
autotoxicity indicates that one plant releases allelochemicals into the adjacent environment, which directly or indirectly inhibit their own normal growth or that of their relatives in continuous cropping or short rotation .
autotoxicity commonly occurs both in nature and cultivation and affects the yield and quality of plants .
the phenomenon was thought to be caused by soil-borne pathogens that gradually accumulated during continuous cropping at first.
however, evidence showed that secondary micromolecular metabolites released by plants are important factors  <cit> .
autotoxicity has been found in many plants, such as rice  <cit> , wheat  <cit> , sweet potato  <cit> , apple  <cit>  and alfalfa  <cit> .
recently, medicinal plants were also reported to exhibit autotoxic phenomena .
autotoxicity has been reported in more than  <dig> plant species thus far  <cit> .
potential autotoxins have been explored for different compounds.
for example, water-soluble organic acids, aliphatic aldehydes, lactones, long-chain fatty acids, naphthoquinones, anthraquinones, coumarins, tannins, steroids, alkaloids, cyanohydrins, sulfides, oil glycosides and purines have been isolated and identified as allelochemicals  <cit> .
other types of autotoxins including phenolic compounds, flavones and terpenoids, and phenolic compounds are the main types of autotoxins  <cit> .
the main phenolic acids are ferulic acid, o-hydroxy phenyl acetic acid, p-coumaric acid, ferulic acid, iso-ferulic acid, malic, critic, fumatic, and caffeic acids .
microarrays have been used to investigate the gene expression profiling in plants in response to allelochemicals such as 2-benzoxazolinone  <cit> , gallic acid  <cit> , 3--l-alanine  <cit> , juglone  <cit>  and ferulic acid  <cit> .
however, dna microarrays have limitations including background hybridization and known probe information.
the recently developed next generation high throughput rna sequencing  technique directly sequences transcripts to obtain a large dynamic range of expression levels and may provide comprehensive clues to solve this problem  <cit> .
furthermore, rna-seq can be used directly for transcriptome profiling for species with  <cit>  or without  <cit>  available genomes.
therefore, rna-seq has potential to replace microarray technology  <cit> .
panax ginseng is a highly valuable perennial herb with medicinal properties that is native to china and korea  <cit>  .
however, continuous cropping of p. ginseng results in autotoxicity and a decline in biomass.
to date, transcriptome studies of p. ginseng have focused on ginsenoside biosynthetic genes .
transcriptome studies of p. ginseng after autotoxin treatment have not been reported.
benzoic acid is one of the major autotoxins identified in p. ginseng root exudates and rhizosphere soil.
it significantly inhibited seed germination and growth  <cit> .
although the autotoxins in p. ginseng have been isolated and identified, their mechanisms remain unknown.
in this study, we constructed rna-seq libraries using rna extracted from roots, stems, and leaves treated with benzoic acid at six different time points and revealed enriched functional terms in response to this stress.
interestingly, several transcript factors were up-regulated in roots, suggesting the importance of transcription factors in response to benzoic acid.
moreover, peroxidase  and superoxide dismutase  response to benzoic acid induction were identified.
several key contigs involved in the flavonoid and ginsenoside biosynthesis pathways were repressed.
these results provide a comprehensive understanding of the p. ginseng response to benzoic acid stress and lay a foundation for improving the resistance or endurance of p. ginseng to autotoxins in the environment.
RESULTS
rna sequencing and de novo assembly
rna samples were collected from roots, stems and leaves at 0 days post treatment ,  <dig> dpt,  <dig> dpt,  <dig> dpt,  <dig> dpt and  <dig> dpt after benzoic acid treatment.
then, rna-seq was performed to investigate degs or pathway responses to benzoic acid.
in total,  <dig> rna-seq libraries generated approximately  <dig>   <dig> clean reads of 100 nt in length .
all the sequencing reads are deposited in the ncbi short read archive  under the accession number srp <dig>  the reads were pooled together and assembled into reference sequences, yielding  <dig>  contigs representing  <dig> , <dig> total assembled bases for the p. ginseng transcriptome.
the gc content was  <dig>  %.
the n <dig> and average length of assembled sequences were 1794 bp and 1428 bp, respectively .
this transcriptome assembly project has been deposited in ddbj/embl/genbank under the accession gdqw <dig> table  <dig> summary of illumina sequencing and transcriptome assemblies for rna-seq libraries
annotation of the ginseng transcriptome
the assembled transcripts were searched against the sequences in the nt database using the blastn algorithm  for functional annotations.
to obtain comprehensive annotation, coding regions were extracted using trinity software  <cit> .
these protein sequences were searched against the uniprot and nr databases using the blastp algorithm .
complete annotation information was listed in additional file  <dig>  a total of  <dig> ,  <dig> , and  <dig>  genes were annotated in the uniprot, nr and nt databases, respectively.
finally,  <dig>  genes were annotated, and the nt database had the largest match, followed by the nr and uniprot databases.
blast was also performed against the kegg database to annotate the metabolic pathways for each gene.
a total of  <dig> genes in  <dig> pathways were identified according to the kegg database .
calculation of gene expression
gene expression was measured by mapping rna-seq reads from  <dig> libraries to the assembled sequences.
on average, 85 % of total reads were successfully mapped to the assembled transcript sequences using bowtie <dig>  <dig> .
<dig>  <cit> .
subsequently, fpkm was adopted to quantify the expression of  <dig>  contigs.
both pairs of pe reads with unique location were retained to calculate fpkm value in the following analysis to detect the degs associated with benzoic acid stress.
there were  <dig> , <dig>  reads that were mapped back in pairs with unique locations on the assembled genes and used for the downstream fpkm calculation.
the average fpkm value for all samples examined was  <dig>  additional file  <dig> shows the distributions of fpkm values followed by log <dig> transformation.
in total,  <dig>  genes with the minimum expression threshold of fpkm >  <dig> were identified in at least one rna-seq library.
identification of degs and dynamic profile
each deg was evaluated by the ma-plot-based method with random sampling model  <cit>  to calculate the p value by the pair-wise comparison of control samples and benzoic acid induced samples.
the stringent cut-offs for degs included a three-fold or higher change and p <  <dig> .
a total of  <dig> and  <dig> genes were identified with significantly up- and down-regulated expression, respectively, compared with control .
all the differential expressed genes could be found at additional file  <dig>  to validate the rna-seq results, qrt-pcr was carried out on  <dig> randomly chosen degs at  <dig> dpt,  <dig> dpt,  <dig> dpt,  <dig> dpt,  <dig> dpt and  <dig> dpt.
the expression profiles based on rna-seq could be found at additional file  <dig>  the qrt-pcr results were highly consistent with those from the rna-seq .
the pearson correlation coefficient was  <dig> , indicating that the qrt-pcr results were positively linearly related with the rna-seq results .
the rna-seq profiles of  <dig> randomly chosen degs and all the genes that described in the result parts were validated by qrt-pcr .
genes with similar expression patterns might be functionally correlated  <cit> .
in order to reveal the expression change after benzoic acid treatment, six time points and three tissues across  <dig> rna-seq libraries were used to identify statistically significant expression profiles and the genes associated with the profiles.
to investigate the dynamic patterns of these degs after benzoic acid treatment, gene expression cluster analysis was performed using the stem algorithms  <cit> .
in brief, the stem algorithm first selects a set of representative model profiles.
then the correlation coefficient for each gene was used to identify the most closely model profile.
then the profiles with statistically significant higher number of genes were selected using permutation test.
stem algorithm analyses of the expressed genes generated three significant clusters  based on the expression profile similarity.
cluster  <dig> included  <dig> genes which increased sharply at early time points  after benzoic acid treatment in roots .
genes in this cluster were predominantly enriched in response to toxin metabolic process  and response to stress  based on go analysis.
however, for cluster  <dig> in the stems and cluster  <dig> in the leaves, the maximum abundance shifted to later time points  .
the gene numbers for these two clusters were  <dig> and  <dig>  respectively.
go terms enriched in cluster  <dig> also included response to stress  and response to stimulus , suggesting the importance of the genes that response to stress after benzoic acid treatment.
within contrast to these two clusters, genes in cluster  <dig> were enriched in nuclease activity  and dna metabolic process , respectively.table  <dig> the number of up- and down-regulated degs based on pair-wise comparison with 0 day control
fig.
<dig> expression patterns of  <dig> randomly selected degs by qrt-pcr and the correlation between qrt-pcr and rna-seq.
fold changes of the transcript levels at different time points are shown.
the average expression level in  <dig> dpt was set to  <dig>  error bars represent standard error.
l, s and r represent leaves, stems and roots, respectively.
the x-axis represents log <dig> fold change of the qrt-pcr.
the y-axis indicates the log <dig> value of the fold change from rna-seq.
the size of each point is proportional to the log <dig>  in  <dig> dpt.
leaves, stems and roots were represented by points with different color
fig.
<dig> representative time-course profile clusters.
representative profile clusters for roots, stems and leaves.
for each time point and each gene, the log <dig>  is shown for the mean expression profile for each cluster
fig.
<dig> venn diagram illustrating the overlap of the degs
degs overlap at six different time points
degs that overlapped at all time points were genes induced by benzoic acid across the entire stress process.
in total, seven contigs were induced across all time points  in roots .
among these, the eg45-like domain containing protein  was significantly up-regulated in roots.
the qrt-pcr validation for the gene was consistent with the result of rna-seq .
in this study,  <dig> contigs including five transcription factors were up-regulated across all time points in stems .
the function of these transcription factors was deduced from the homology protein of uniprot database in other species.
two myb transcription factors  can respond to ethylene, jasmonic acid and abscisic acid.
two no apical meristem  transcription factors  functioned in multidimensional cell growth and the negative regulation of the abscisic acid activated signaling pathway.
transcription factor wrky  played a critical role in the response to chitin.
the expression profiles based rna-seq and qrt-pcr validation for the above-mentioned transcript factors could be found at additional file  <dig>
in the leaves, seven genes were up-regulated across all the time points .
the plant invertase/pectin methylesterase inhibitor has been reported as a stress-response gene  <cit> .
in this study, we also found that one plant invertase/pectin methylesterase inhibitor  gene was up-regulated in the leaves .
we also observed down-regulated genes across all the time points.
for example, one unigene that coded for zinc-binding dehydrogenase  was down-regulated in roots, and another unigene that coded for calcium-binding protein with ef-hand domain  was down-regulated in stems .
these results are consistent with previous reports that the lipophilicity of benzoic and cinnamic acids could impact the ion uptake and root elongation in cucumber  <cit> .
all the differentially expressed genes at all time points could be found at additional file  <dig>
go functional enrichment analyses
the function of the  <dig> up-regulated degs was determined by go enrichment analysis.
the results showed that the up-regulated genes were enriched for ‘transcription factor activity’, ‘carboxylesterase activity’, and ‘pectinesterase activity’ .
up-regulated genes were also involved in ‘centrosome cycle and duplication’ and ‘mitotic spindle elongation’.
other enriched terms included ‘response to water deprivation’ , consistent with previous findings that exposure of plant roots to ferulic acid can affect water usage  <cit> .
we also found that ‘polysaccharide metabolic process’  and ‘lignin metabolic process’  were enriched.fig.
<dig> most significantly enriched go terms  of the degs.
the most significant go terms  are presented graphically.
the x-axis represents the log <dig> of the enrichment p value.
the y-axis indicates the number of degs in log <dig> value.
the size of each point is proportional to the percentage .
a and b presents the enriched go term for up-regulation genes based on biological procession and molecular function, respectively.
c and d presents the enriched go term for down-regulation genes based on biological procession and molecular function, respectively
among the down-regulated genes, go terms associated with ‘cell wall organization or biogenesis’, ‘stress response’, ‘secondary metabolic processes’ and ‘ion binding’ were enriched .
down-regulated genes involved in ‘polysaccharide metabolic process’  and ‘lignin metabolic process’  showed significant responses to benzoic acid.
autotoxin usually led to a 20-50 % decrease in growth rate  <cit> .
the photosystem is an important factor that influences plant growth.
contigs involved in ‘photosystem’  were significantly repressed, which can partly explain the reduction of production due to autotoxicity.
reactive oxygen species
reactive oxygen species  produced from various metabolic pathways can cause oxidative damage  <cit> .
for example, ferulic acid can produce oxidative stress in rice  <cit> .
plants can produce reactive oxygen species  during abiotic stress as byproducts of aerobic metabolism.
at the same time plant must have the ability to remove the toxin.
benzoic acid treatment can be regarded as abiotic stress, thus the investigation of the generation and removal of ros can help us to understand the possible function of ros upon benzoic acid treatment.
in this study,  <dig> ros-related genes were annotated as ‘response to reactive oxygen species’ .
among them, three contigs  were significantly up-regulated, suggesting that these genes might take part in the oxidative damage response upon benzoic acid induction.
in leaves, c14304_g <dig> was up-regulated at an early stage .
however, c42911_g <dig> and c34265_g <dig> were up-regulated at a later stage in roots .
the rna-seq profiles and qrt-pcr validations for the above three ros-related genes could be found at additional file  <dig>
we also investigate the differential expression levels of contigs coding pod and sod.
in total, we identified  <dig> contigs with peroxidase activity.
two pods  were up-regulated at  <dig> dpt in roots.
we also identified nine contigs with sod activity.
among them, c41673_g <dig> was up-regulated in roots  and leaves .
the rna-seq profiles and qrt-pcr validations for two pods genes  and one sod genes  could be found at additional file  <dig>
measurement of physiological and biochemical indexes of proline and malondialdehyde
the impacts of benzoic acid over a concentration gradient of  <dig> -250 mg · l− <dig> on the accumulation of proline and malondialdehyde  was investigated.
mda is the product of lipid membrane peroxidation and is a marker of stress-induced damage  <cit> .
in general, the content of both free proline and mda increased as the concentration of benzoic acid increased compared with control .
the accumulation of proline at  <dig> dpt was higher than at  <dig>   <dig>   <dig>   <dig> and  <dig> dpt.
in the proline biosynthetic pathway, pyrroline-5-carboxylate synthetase  and pyrroline-5-carboxylate reductase  were two important enzymes with important roles in turning glutamic acid into proline  <cit> .
the expression level of p5cs  showed a similar profile as proline.
however, the expression level of p5cr  was unchanged during benzoic acid treatment.
for mda , the content at  <dig> dpt was higher than at  <dig>   <dig>   <dig>   <dig> and  <dig> dpt, suggesting that benzoic acid treatment promoted lipid peroxidation at later stage.
ascorbate is catalyzed and converted into mda by ascorbate peroxidase   <cit> .
according to the rna-seq results, the expression of apx  increased at the early stage and decreased at the later stage, consistent with the change in mda content.
the rna-seq profiles for p5cs, p5cr and apx, together with the qrt-pcr validations were shown at additional file  <dig> fig.
<dig> measurement of physiological and biochemical indexes of proline  and malondialdehyde
effect of autotoxicity on root hair development
autotoxins produced by plants inhibit their own growth directly or indirectly.
for example, ferulic acid inhibited root elongation in rice  <cit> .
therefore, we also observed the root morphology after benzoic acid treatment.
the results showed that benzoic acid caused permanent morphological reductions in root hair density by observation of root morphology using stereomicroscope  <dig>  the adverse effect of autotoxins on root development increased as the concentration of benzoic acid increased .
the root hair growth was almost completely inhibited when the concentration of benzoic acid reached 250 mg · l− <dig>  suggesting that benzoic acid repressed the root hair in a concentration-dependent manner.
the possible mechanism of the effect of benzoic acid on root hair formation was investigated.
in total,  <dig> root hair related contigs were identified from the homology protein of uniprot database in other species, and the contigs  were highly expressed in roots compared with leaves and stems .
among these root hair related contigs, four up-regulated  and two down-regulated  contigs might contribute to the low density of root hairs upon benzoic acid treatment.
c35599_g <dig> was up-regulated at  <dig> dpt in roots; both c44448_g <dig> and c47561_g <dig> were up-regulated in stems  and c54588_g <dig> was up-regulated in roots .
two down-regulated contigs were observed in stems.
the qrt-pcr validations with rna-seq profiles for the above nine genes could be found at additional file  <dig> fig.
<dig> plot of the  <dig> root hair related genes according to fpkm values
flavonoid and ginsenoside biosynthetic pathways upon benzoic acid treatment
previous studies revealed that flavonoids were also involved in autotoxicity and resulted in root growth inhibition and cell division reduction in root meristematic regions  <cit> .
chalcone synthase , chalcone isomerase , flavanone 3-hydroxylase , flavonoid 3’ 5’-hydroxylase , flavonoid 3’-hydroxylase , anthocyanidin synthase , and udp-glucose:anthocyanidin 3-o-glucosyltransferase  genes are key enzymes in the flavonoid biosynthesis pathway  <cit> .
in this study, we identified one chs unigene , two chi contigs , one f3h contigs  and two a3gt contigs .
interestingly, chi  was repressed in stems at the later stage  after benzoic acid treatment.
f3h  was also repressed in roots at the later stage .
the rna-seq profiles and qrt-pcr validations for chi and f3h could be found at additional file  <dig>
ginsenosides are the major bioactive components of p. ginseng.
due to their important pharmacological effects, the ginsenoside biosynthetic pathway has been widely investigated .
cytochrome p <dig> involved in the ginsenoside biosynthetic pathway might represent the largest family and catalyze most of the oxidation steps during the biogenesis of plant secondary metabolism.
in this study,  <dig> potential cytochrome p <dig> contigs were down-regulated and  <dig> were up-regulated.
the remaining contigs involved in the ginsenoside biosynthetic pathway were down-regulated upon benzoic acid treatment.
the down-regulated degs included three farnesyl diphosphates , two-geranylgeranyl pyrophosphate synthase , one hydroxymethylglutaryl-coa synthase  and two ugt .
the rna-seq profiles and qrt-pcr validations for the two fps , two gps , one hmgs  and two ugt  for chi and f3h could be found at additional file  <dig>
the expression of transcription factors upon benzoic acid treatment
tfs function in the abiotic stress tolerance response, and overexpression of these tfs generated transgenic plants with significant tolerance to abiotic stresses .
a total of  <dig> tfs were significantly up-regulated by benzoic acid, and tf expression levels were higher in roots than in other tissues .
in total, 10 % of tfs in the tf family were up-regulated in roots .
our results revealed that the expression of three ap2/erf , four myb , and nine wrky genes  were up-regulated in response to benzoic acid, suggesting that responses to benzoic acid might require the up-regulated tfs in roots to regulate downstream response genes.
the qrt-pcr validations and rna-seq profiles for these tfs could be found at additional file  <dig> fig.
<dig> hierarchical clustering of differentially expressed tfs from  <dig> libraries.
hierarchical clustering of the degs among roots, stems and leaves at  <dig> dpt,  <dig> dpt,  <dig> dpt,  <dig> dpt,  <dig> dpt and  <dig> dpt using fpkm expression values from rna-seq.
each column represents an rna-seq library and each row represents a deg.
different colors shown different expression levels, with green representing low expression levels and red representing high expression levels
alternative splicing of the transcriptome upon benzoic acid induction
in eukaryotes, alternative splicing  increases the diversity of the transcriptome and the proteome and takes part in biological process such as plant-virus interactions  <cit>  and stress-related regulation  <cit> .
however, the as events upon benzoic acid treatment have never been reported.
we used previous published method to identify as events in a genome-wide way  <cit> .
in brief, we used blat to compare trinity assembled sequences with the following options: −tilesize =  <dig> -oneoff =  <dig> -minidentity =  <dig> -out = sim <dig> -maxintron =  <dig>  then the candidates were revealed according to the pair-wise comparisons to search for indels.
finally, indels that were flanked by two sequentially matching sections were regarded as candidate as regions.
in total,  <dig>  as events in  <dig>  as genes were found upon benzoic acid induction.
as events from  <dig> randomly selected genes were validated by rt-pcr.
then the expected band sizes can be predicted according to the rna-seq results and were listed at additional file  <dig>  the real band sizes from rt-pcr that matched the expected band size from rna-seq were arrowed in fig.
<dig>  for genes c15452_g <dig> and c11325_g <dig>  we only observed one clear expected band, however, the other weak bands might be caused by low abundance.
the other  <dig> genes can be detected all expected bands and the obtained production sizes of rt-pcr were consistent with that of the rna-seq results.
splicing factors are extensively alternatively spliced.
in this study, we found  <dig> splicing factors.
among these splicing factors,  <dig> splicing factors were regulated by as.
in total, the as splicing factor percentage was 54 %, higher than the average as ratio .
among the contigs involved in the ginsenoside biosynthetic pathway, one aact, one as,  <dig> fps,  <dig> gps,  <dig> hmgr,  <dig> hmgs,  <dig> se,  <dig> ss and  <dig> ugt contigs were regulated by as.fig.
<dig> validation of randomly selected as events by rt-pcr
discussion
autotoxicity reduced both yield and quality of crops in a continuous cropping pattern  <cit> .
p. ginseng, one of the most important traditional chinese herbal medicines, exhibited a severe autotoxicity problem.
in our previous study, secondary metabolites including benzoic acid were identified with strong autotoxic and growth-inhibiting activity on p. ginseng  <cit> .
however, the molecular mechanism by which p. ginseng responds to benzoic acid remains unknown.
in this study, we investigated the effect of benzoic acid on gene expression profiles using rna-seq technology and identified the differentially expressed genes that directly respond to autotoxins and the enriched go terms that were involved in multiple biological processes.
moreover, we identified three up-regulated ros-related genes.
ros can damage cell membranes and other components.
the ros-scavenging system involves three major enzymes: sod, pod and cat  <cit> .
in this study, we found that sod and pod contigs were up- regulated at later stages, which might enhance the anti-oxidation mechanisms at the later stage of benzoic acid induction.
the physiological and biochemical indexes for proline and mda also support this conclusion.
consistent with the gene expression profiles of sod and pod, the proline and mda content decreased at the later stage of benzoic acid treatment.
both free proline and mda were restored to the lower level at a later stage , implying that p. ginseng partly became resistant to benzoic acid stress at the later stage.
we also found that two important contigs that coded for the enzymes chi and f3h in the flavonoid biosynthesis pathway were repressed at the later stage.
flavonoid production might be repressed at the later stage by the down-regulation of chi and f3h upon benzoic acid treatment.
in this way, flavonoid could be presented at a relatively low concentration, and the exudation of flavonoid into soil could be reduced.
the reduced concentrations of flavonoid in the soil can the release some of the adverse effects of autotoxin.
in the future, the actual flavonoid concentrations upon benzoic acid treatment will be measured to obtain accurate information about the regulation of flavonoid biogenesis and transporters in the response to autotoxins.
autotoxicity and pathogens are two major causes of continuous cropping  <cit> .
thus, several groups have investigated the interplay between autotoxins and pathogens in soil.
autotoxins can be degraded by pathogens  <cit>  and repress the prevalence of soil pathogens  <cit> .
however, another study showed that autotoxins exuded from the root of ginseng saponins stimulated the growth of soil pathogens  <cit> .
baicalin released from scutellaria baicalensis georgi increases the activity of soil pathogens  <cit> .
the suppression or enhancement effect of autotoxins on pathogens can be explained by the diversity of responses of plants to different autotoxins and pathogens.
the accumulation of c. destructans, one of the most destructive pathogens, results in serious root rot disease and prevents the continuous cultivation of p. ginseng  <cit> .
autotoxins in the rhizosphere in p. ginseng affect and interact with pathogens in soils.
however, a study on this topic has never been reported in p. ginseng.
we still do not know the suppression or enhancement effects of autotoxins on the most destructive pathogens.
in this study, we found that ‘stress response genes’ was enriched in differentially expressed genes after exposure to benzoic acid.
these genes might also affect disease resistance upon cylindrocarpon destructans infection.
understanding the interaction underlying autotoxins and pathogens is necessary to effectively control both p. ginseng root diseases and the adverse effects of autotoxicity.
previous studies showed that the autotoxin potential of many plants is species-dependent  <cit> .
the autotoxin was toxic to cucumber but not to fig leaf gourd plants  <cit> .
in this study, we only uncovered the transcriptional response induced by benzoic acid.
thus, in the future it will be interesting to investigate different p. ginseng autotoxins to reveal the similarities and differences in the molecular mechanisms of autotoxin-induced responses.
CONCLUSIONS
p. ginseng is a suitable plant for the study of autotoxicity because it is a highly valuable perennial herb and its biomass is affected by continuous cropping.
based on  <dig> rna-seq libraries from roots, stems, and leaves at different time points, we obtained transcriptome profiles and identified the degs upon benzoic acid stress.
more importantly, changes in transcription factors upon benzoic acid induction were identified.
additionally, the study analyzed the reactive oxygen species genes and measured the physiological and biochemical indices of proline and mda, which will help to reveal the molecular mechanisms of autotoxicity.
moreover, several differentially expressed root hair related genes that coded for key enzymes in the flavonoid biosynthesis pathway were identified, which might contribute to the dose-dependent inhibitory effect of benzoic acid on root hair development.
finally, the down-regulated contigs in the ginsenoside biosynthetic pathway were identified, and the repression of these contigs can reduce the biosynthesis of ginsenoside and partly release the adverse effects of autotoxins.
this study provided a comprehensive characterization of the transcriptional response induced by autotoxins.

 autotoxicity indicates that one plant releases allelochemicals into the adjacent environment, which directly or indirectly inhibit their own normal growth or that of their relatives in continuous cropping or short rotation. autotoxicity commonly occurs both in nature and cultivation and affects the yield and quality of plants. the phenomenon was thought to be caused by soil-borne pathogens that gradually accumulated during continuous cropping at first. however, evidence showed that secondary micromolecular metabolites released by plants are important factors. autotoxicity has been found in many plants, such as rice, wheat, sweet potato, apple and alfalfa. recently, medicinal plants were also reported to exhibit autotoxic phenomena. autotoxicity has been reported in more than plant species thus far. potential autotoxins have been explored for different compounds. for example, water-soluble organic acids, aliphatic aldehydes, lactones, long-chain fatty acids, naphthoquinones, anthraquinones, coumarins, tannins, steroids, alkaloids, cyanohydrins, sulfides, oil glycosides and purines have been isolated and identified as allelochemicals. other types of autotoxins including phenolic compounds, flavones and terpenoids, and phenolic compounds are the main types of autotoxins. the main phenolic acids are ferulic acid, o-hydroxy phenyl acetic acid, p-coumaric acid, ferulic acid, iso-ferulic acid, malic, critic, fumatic, and caffeic acids. microarrays have been used to investigate the gene expression profiling in plants in response to allelochemicals such as 2-benzoxazolinone, gallic acid, 3--l-alanine, juglone and ferulic acid. however, dna microarrays have limitations including background hybridization and known probe information. the recently developed next generation high throughput rna sequencing technique directly sequences transcripts to obtain a large dynamic range of expression levels and may provide comprehensive clues to solve this problem. furthermore, rna-seq can be used directly for transcriptome profiling for species with or without available genomes. therefore, rna-seq has potential to replace microarray technology. panax ginseng is a highly valuable perennial herb with medicinal properties that is native to china and korea. however, continuous cropping of p. ginseng results in autotoxicity and a decline in biomass. to date, transcriptome studies of p. ginseng have focused on ginsenoside biosynthetic genes. transcriptome studies of p. ginseng after autotoxin treatment have not been reported. benzoic acid is one of the major autotoxins identified in p. ginseng root exudates and rhizosphere soil. it significantly inhibited seed germination and growth. although the autotoxins in p. ginseng have been isolated and identified, their mechanisms remain unknown. in this study, we constructed rna-seq libraries using rna extracted from roots, stems, and leaves treated with benzoic acid at six different time points and revealed enriched functional terms in response to this stress. interestingly, several transcript factors were up-regulated in roots, suggesting the importance of transcription factors in response to benzoic acid. moreover, peroxidase and superoxide dismutase response to benzoic acid induction were identified. several key contigs involved in the flavonoid and ginsenoside biosynthesis pathways were repressed. these results provide a comprehensive understanding of the p. ginseng response to benzoic acid stress and lay a foundation for improving the resistance or endurance of p. ginseng to autotoxins in the environment.  rna sequencing and de novo assembly rna samples were collected from roots, stems and leaves at 0 days post treatment, dpt, dpt, dpt, dpt and dpt after benzoic acid treatment. then, rna-seq was performed to investigate degs or pathway responses to benzoic acid. in total, rna-seq libraries generated approximately clean reads of 100 nt in length. all the sequencing reads are deposited in the ncbi short read archive under the accession number srp the reads were pooled together and assembled into reference sequences, yielding contigs representing, total assembled bases for the p. ginseng transcriptome. the gc content was  %. the n and average length of assembled sequences were 1794 bp and 1428 bp, respectively. this transcriptome assembly project has been deposited in ddbj/embl/genbank under the accession gdqw table summary of illumina sequencing and transcriptome assemblies for rna-seq libraries annotation of the ginseng transcriptome the assembled transcripts were searched against the sequences in the nt database using the blastn algorithm for functional annotations. to obtain comprehensive annotation, coding regions were extracted using trinity software. these protein sequences were searched against the uniprot and nr databases using the blastp algorithm. complete annotation information was listed in additional file a total of,, and genes were annotated in the uniprot, nr and nt databases, respectively. finally, genes were annotated, and the nt database had the largest match, followed by the nr and uniprot databases. blast was also performed against the kegg database to annotate the metabolic pathways for each gene. a total of genes in pathways were identified according to the kegg database. calculation of gene expression gene expression was measured by mapping rna-seq reads from libraries to the assembled sequences. on average, 85 % of total reads were successfully mapped to the assembled transcript sequences using bowtie.. subsequently, fpkm was adopted to quantify the expression of contigs. both pairs of pe reads with unique location were retained to calculate fpkm value in the following analysis to detect the degs associated with benzoic acid stress. there were, reads that were mapped back in pairs with unique locations on the assembled genes and used for the downstream fpkm calculation. the average fpkm value for all samples examined was additional file shows the distributions of fpkm values followed by log transformation. in total, genes with the minimum expression threshold of fpkm >  were identified in at least one rna-seq library. identification of degs and dynamic profile each deg was evaluated by the ma-plot-based method with random sampling model to calculate the p value by the pair-wise comparison of control samples and benzoic acid induced samples. the stringent cut-offs for degs included a three-fold or higher change and p < . a total of and genes were identified with significantly up- and down-regulated expression, respectively, compared with control. all the differential expressed genes could be found at additional file to validate the rna-seq results, qrt-pcr was carried out on randomly chosen degs at dpt, dpt, dpt, dpt, dpt and dpt. the expression profiles based on rna-seq could be found at additional file the qrt-pcr results were highly consistent with those from the rna-seq. the pearson correlation coefficient was, indicating that the qrt-pcr results were positively linearly related with the rna-seq results. the rna-seq profiles of randomly chosen degs and all the genes that described in the result parts were validated by qrt-pcr. genes with similar expression patterns might be functionally correlated. in order to reveal the expression change after benzoic acid treatment, six time points and three tissues across rna-seq libraries were used to identify statistically significant expression profiles and the genes associated with the profiles. to investigate the dynamic patterns of these degs after benzoic acid treatment, gene expression cluster analysis was performed using the stem algorithms. in brief, the stem algorithm first selects a set of representative model profiles. then the correlation coefficient for each gene was used to identify the most closely model profile. then the profiles with statistically significant higher number of genes were selected using permutation test. stem algorithm analyses of the expressed genes generated three significant clusters based on the expression profile similarity. cluster included genes which increased sharply at early time points after benzoic acid treatment in roots. genes in this cluster were predominantly enriched in response to toxin metabolic process and response to stress based on go analysis. however, for cluster in the stems and cluster in the leaves, the maximum abundance shifted to later time points. the gene numbers for these two clusters were and respectively. go terms enriched in cluster also included response to stress and response to stimulus, suggesting the importance of the genes that response to stress after benzoic acid treatment. within contrast to these two clusters, genes in cluster were enriched in nuclease activity and dna metabolic process, respectively.table the number of up- and down-regulated degs based on pair-wise comparison with 0 day control fig. expression patterns of randomly selected degs by qrt-pcr and the correlation between qrt-pcr and rna-seq. fold changes of the transcript levels at different time points are shown. the average expression level in dpt was set to error bars represent standard error. l, s and r represent leaves, stems and roots, respectively. the x-axis represents log fold change of the qrt-pcr. the y-axis indicates the log value of the fold change from rna-seq. the size of each point is proportional to the log in dpt. leaves, stems and roots were represented by points with different color fig. representative time-course profile clusters. representative profile clusters for roots, stems and leaves. for each time point and each gene, the log is shown for the mean expression profile for each cluster fig. venn diagram illustrating the overlap of the degs degs overlap at six different time points degs that overlapped at all time points were genes induced by benzoic acid across the entire stress process. in total, seven contigs were induced across all time points in roots. among these, the eg45-like domain containing protein was significantly up-regulated in roots. the qrt-pcr validation for the gene was consistent with the result of rna-seq. in this study, contigs including five transcription factors were up-regulated across all time points in stems. the function of these transcription factors was deduced from the homology protein of uniprot database in other species. two myb transcription factors can respond to ethylene, jasmonic acid and abscisic acid. two no apical meristem transcription factors functioned in multidimensional cell growth and the negative regulation of the abscisic acid activated signaling pathway. transcription factor wrky played a critical role in the response to chitin. the expression profiles based rna-seq and qrt-pcr validation for the above-mentioned transcript factors could be found at additional file in the leaves, seven genes were up-regulated across all the time points. the plant invertase/pectin methylesterase inhibitor has been reported as a stress-response gene. in this study, we also found that one plant invertase/pectin methylesterase inhibitor gene was up-regulated in the leaves. we also observed down-regulated genes across all the time points. for example, one unigene that coded for zinc-binding dehydrogenase was down-regulated in roots, and another unigene that coded for calcium-binding protein with ef-hand domain was down-regulated in stems. these results are consistent with previous reports that the lipophilicity of benzoic and cinnamic acids could impact the ion uptake and root elongation in cucumber. all the differentially expressed genes at all time points could be found at additional file go functional enrichment analyses the function of the up-regulated degs was determined by go enrichment analysis. the results showed that the up-regulated genes were enriched for ‘transcription factor activity’, ‘carboxylesterase activity’, and ‘pectinesterase activity’. up-regulated genes were also involved in ‘centrosome cycle and duplication’ and ‘mitotic spindle elongation’. other enriched terms included ‘response to water deprivation’, consistent with previous findings that exposure of plant roots to ferulic acid can affect water usage. we also found that ‘polysaccharide metabolic process’ and ‘lignin metabolic process’ were enriched.fig. most significantly enriched go terms of the degs. the most significant go terms are presented graphically. the x-axis represents the log of the enrichment p value. the y-axis indicates the number of degs in log value. the size of each point is proportional to the percentage. a and b presents the enriched go term for up-regulation genes based on biological procession and molecular function, respectively. c and d presents the enriched go term for down-regulation genes based on biological procession and molecular function, respectively among the down-regulated genes, go terms associated with ‘cell wall organization or biogenesis’, ‘stress response’, ‘secondary metabolic processes’ and ‘ion binding’ were enriched. down-regulated genes involved in ‘polysaccharide metabolic process’ and ‘lignin metabolic process’ showed significant responses to benzoic acid. autotoxin usually led to a 20-50 % decrease in growth rate. the photosystem is an important factor that influences plant growth. contigs involved in ‘photosystem’ were significantly repressed, which can partly explain the reduction of production due to autotoxicity. reactive oxygen species reactive oxygen species produced from various metabolic pathways can cause oxidative damage. for example, ferulic acid can produce oxidative stress in rice. plants can produce reactive oxygen species during abiotic stress as byproducts of aerobic metabolism. at the same time plant must have the ability to remove the toxin. benzoic acid treatment can be regarded as abiotic stress, thus the investigation of the generation and removal of ros can help us to understand the possible function of ros upon benzoic acid treatment. in this study, ros-related genes were annotated as ‘response to reactive oxygen species’. among them, three contigs were significantly up-regulated, suggesting that these genes might take part in the oxidative damage response upon benzoic acid induction. in leaves, c14304_g was up-regulated at an early stage. however, c42911_g and c34265_g were up-regulated at a later stage in roots. the rna-seq profiles and qrt-pcr validations for the above three ros-related genes could be found at additional file we also investigate the differential expression levels of contigs coding pod and sod. in total, we identified contigs with peroxidase activity. two pods were up-regulated at dpt in roots. we also identified nine contigs with sod activity. among them, c41673_g was up-regulated in roots and leaves. the rna-seq profiles and qrt-pcr validations for two pods genes and one sod genes could be found at additional file measurement of physiological and biochemical indexes of proline and malondialdehyde the impacts of benzoic acid over a concentration gradient of -250 mg · l− on the accumulation of proline and malondialdehyde was investigated. mda is the product of lipid membrane peroxidation and is a marker of stress-induced damage. in general, the content of both free proline and mda increased as the concentration of benzoic acid increased compared with control. the accumulation of proline at dpt was higher than at and dpt. in the proline biosynthetic pathway, pyrroline-5-carboxylate synthetase and pyrroline-5-carboxylate reductase were two important enzymes with important roles in turning glutamic acid into proline. the expression level of p5cs showed a similar profile as proline. however, the expression level of p5cr was unchanged during benzoic acid treatment. for mda, the content at dpt was higher than at and dpt, suggesting that benzoic acid treatment promoted lipid peroxidation at later stage. ascorbate is catalyzed and converted into mda by ascorbate peroxidase. according to the rna-seq results, the expression of apx increased at the early stage and decreased at the later stage, consistent with the change in mda content. the rna-seq profiles for p5cs, p5cr and apx, together with the qrt-pcr validations were shown at additional file fig. measurement of physiological and biochemical indexes of proline  and malondialdehyde  effect of autotoxicity on root hair development autotoxins produced by plants inhibit their own growth directly or indirectly. for example, ferulic acid inhibited root elongation in rice. therefore, we also observed the root morphology after benzoic acid treatment. the results showed that benzoic acid caused permanent morphological reductions in root hair density by observation of root morphology using stereomicroscope the adverse effect of autotoxins on root development increased as the concentration of benzoic acid increased. the root hair growth was almost completely inhibited when the concentration of benzoic acid reached 250 mg · l− suggesting that benzoic acid repressed the root hair in a concentration-dependent manner. the possible mechanism of the effect of benzoic acid on root hair formation was investigated. in total, root hair related contigs were identified from the homology protein of uniprot database in other species, and the contigs were highly expressed in roots compared with leaves and stems. among these root hair related contigs, four up-regulated and two down-regulated contigs might contribute to the low density of root hairs upon benzoic acid treatment. c35599_g was up-regulated at dpt in roots; both c44448_g and c47561_g were up-regulated in stems and c54588_g was up-regulated in roots. two down-regulated contigs were observed in stems. the qrt-pcr validations with rna-seq profiles for the above nine genes could be found at additional file fig. plot of the root hair related genes according to fpkm values flavonoid and ginsenoside biosynthetic pathways upon benzoic acid treatment previous studies revealed that flavonoids were also involved in autotoxicity and resulted in root growth inhibition and cell division reduction in root meristematic regions. chalcone synthase, chalcone isomerase, flavanone 3-hydroxylase, flavonoid 3’ 5’-hydroxylase, flavonoid 3’-hydroxylase, anthocyanidin synthase, and udp-glucose:anthocyanidin 3-o-glucosyltransferase genes are key enzymes in the flavonoid biosynthesis pathway. in this study, we identified one chs unigene, two chi contigs, one f3h contigs and two a3gt contigs. interestingly, chi was repressed in stems at the later stage after benzoic acid treatment. f3h was also repressed in roots at the later stage. the rna-seq profiles and qrt-pcr validations for chi and f3h could be found at additional file ginsenosides are the major bioactive components of p. ginseng. due to their important pharmacological effects, the ginsenoside biosynthetic pathway has been widely investigated. cytochrome p involved in the ginsenoside biosynthetic pathway might represent the largest family and catalyze most of the oxidation steps during the biogenesis of plant secondary metabolism. in this study, potential cytochrome p contigs were down-regulated and were up-regulated. the remaining contigs involved in the ginsenoside biosynthetic pathway were down-regulated upon benzoic acid treatment. the down-regulated degs included three farnesyl diphosphates, two-geranylgeranyl pyrophosphate synthase, one hydroxymethylglutaryl-coa synthase and two ugt. the rna-seq profiles and qrt-pcr validations for the two fps, two gps, one hmgs and two ugt for chi and f3h could be found at additional file the expression of transcription factors upon benzoic acid treatment tfs function in the abiotic stress tolerance response, and overexpression of these tfs generated transgenic plants with significant tolerance to abiotic stresses. a total of tfs were significantly up-regulated by benzoic acid, and tf expression levels were higher in roots than in other tissues. in total, 10 % of tfs in the tf family were up-regulated in roots. our results revealed that the expression of three ap2/erf, four myb, and nine wrky genes were up-regulated in response to benzoic acid, suggesting that responses to benzoic acid might require the up-regulated tfs in roots to regulate downstream response genes. the qrt-pcr validations and rna-seq profiles for these tfs could be found at additional file fig. hierarchical clustering of differentially expressed tfs from libraries. hierarchical clustering of the degs among roots, stems and leaves at dpt, dpt, dpt, dpt, dpt and dpt using fpkm expression values from rna-seq. each column represents an rna-seq library and each row represents a deg. different colors shown different expression levels, with green representing low expression levels and red representing high expression levels alternative splicing of the transcriptome upon benzoic acid induction in eukaryotes, alternative splicing increases the diversity of the transcriptome and the proteome and takes part in biological process such as plant-virus interactions and stress-related regulation. however, the as events upon benzoic acid treatment have never been reported. we used previous published method to identify as events in a genome-wide way. in brief, we used blat to compare trinity assembled sequences with the following options: −tilesize =  -oneoff =  -minidentity =  -out = sim -maxintron =  then the candidates were revealed according to the pair-wise comparisons to search for indels. finally, indels that were flanked by two sequentially matching sections were regarded as candidate as regions. in total, as events in as genes were found upon benzoic acid induction. as events from randomly selected genes were validated by rt-pcr. then the expected band sizes can be predicted according to the rna-seq results and were listed at additional file the real band sizes from rt-pcr that matched the expected band size from rna-seq were arrowed in fig.  for genes c15452_g and c11325_g we only observed one clear expected band, however, the other weak bands might be caused by low abundance. the other genes can be detected all expected bands and the obtained production sizes of rt-pcr were consistent with that of the rna-seq results. splicing factors are extensively alternatively spliced. in this study, we found splicing factors. among these splicing factors, splicing factors were regulated by as. in total, the as splicing factor percentage was 54 %, higher than the average as ratio. among the contigs involved in the ginsenoside biosynthetic pathway, one aact, one as, fps, gps, hmgr, hmgs, se, ss and ugt contigs were regulated by as.fig. validation of randomly selected as events by rt-pcr discussion autotoxicity reduced both yield and quality of crops in a continuous cropping pattern. p. ginseng, one of the most important traditional chinese herbal medicines, exhibited a severe autotoxicity problem. in our previous study, secondary metabolites including benzoic acid were identified with strong autotoxic and growth-inhibiting activity on p. ginseng. however, the molecular mechanism by which p. ginseng responds to benzoic acid remains unknown. in this study, we investigated the effect of benzoic acid on gene expression profiles using rna-seq technology and identified the differentially expressed genes that directly respond to autotoxins and the enriched go terms that were involved in multiple biological processes. moreover, we identified three up-regulated ros-related genes. ros can damage cell membranes and other components. the ros-scavenging system involves three major enzymes: sod, pod and cat. in this study, we found that sod and pod contigs were up- regulated at later stages, which might enhance the anti-oxidation mechanisms at the later stage of benzoic acid induction. the physiological and biochemical indexes for proline and mda also support this conclusion. consistent with the gene expression profiles of sod and pod, the proline and mda content decreased at the later stage of benzoic acid treatment. both free proline and mda were restored to the lower level at a later stage, implying that p. ginseng partly became resistant to benzoic acid stress at the later stage. we also found that two important contigs that coded for the enzymes chi and f3h in the flavonoid biosynthesis pathway were repressed at the later stage. flavonoid production might be repressed at the later stage by the down-regulation of chi and f3h upon benzoic acid treatment. in this way, flavonoid could be presented at a relatively low concentration, and the exudation of flavonoid into soil could be reduced. the reduced concentrations of flavonoid in the soil can the release some of the adverse effects of autotoxin. in the future, the actual flavonoid concentrations upon benzoic acid treatment will be measured to obtain accurate information about the regulation of flavonoid biogenesis and transporters in the response to autotoxins. autotoxicity and pathogens are two major causes of continuous cropping. thus, several groups have investigated the interplay between autotoxins and pathogens in soil. autotoxins can be degraded by pathogens and repress the prevalence of soil pathogens. however, another study showed that autotoxins exuded from the root of ginseng saponins stimulated the growth of soil pathogens. baicalin released from scutellaria baicalensis georgi increases the activity of soil pathogens. the suppression or enhancement effect of autotoxins on pathogens can be explained by the diversity of responses of plants to different autotoxins and pathogens. the accumulation of c. destructans, one of the most destructive pathogens, results in serious root rot disease and prevents the continuous cultivation of p. ginseng. autotoxins in the rhizosphere in p. ginseng affect and interact with pathogens in soils. however, a study on this topic has never been reported in p. ginseng. we still do not know the suppression or enhancement effects of autotoxins on the most destructive pathogens. in this study, we found that ‘stress response genes’ was enriched in differentially expressed genes after exposure to benzoic acid. these genes might also affect disease resistance upon cylindrocarpon destructans infection. understanding the interaction underlying autotoxins and pathogens is necessary to effectively control both p. ginseng root diseases and the adverse effects of autotoxicity. previous studies showed that the autotoxin potential of many plants is species-dependent. the autotoxin was toxic to cucumber but not to fig leaf gourd plants. in this study, we only uncovered the transcriptional response induced by benzoic acid. thus, in the future it will be interesting to investigate different p. ginseng autotoxins to reveal the similarities and differences in the molecular mechanisms of autotoxin-induced responses.  p. ginseng is a suitable plant for the study of autotoxicity because it is a highly valuable perennial herb and its biomass is affected by continuous cropping. based on rna-seq libraries from roots, stems, and leaves at different time points, we obtained transcriptome profiles and identified the degs upon benzoic acid stress. more importantly, changes in transcription factors upon benzoic acid induction were identified. additionally, the study analyzed the reactive oxygen species genes and measured the physiological and biochemical indices of proline and mda, which will help to reveal the molecular mechanisms of autotoxicity. moreover, several differentially expressed root hair related genes that coded for key enzymes in the flavonoid biosynthesis pathway were identified, which might contribute to the dose-dependent inhibitory effect of benzoic acid on root hair development. finally, the down-regulated contigs in the ginsenoside biosynthetic pathway were identified, and the repression of these contigs can reduce the biosynthesis of ginsenoside and partly release the adverse effects of autotoxins. this study provided a comprehensive characterization of the transcriptional response induced by autotoxins. 
several key contigs in the flavonoid and ginsenoside biosynthesis pathways were repressed. physiological and biochemical indexes showed that the proline and malondialdehyde content were restored to lower levels at a later stage after benzoic acid treatment. transcription factors including ap2s/erfs, mybs, and wrkys were up-regulated in roots after benzoic acid treatment. moreover, reactive oxygen species, peroxidases and superoxide dismutase contigs were up-regulated in roots after benzoic acid treatment. benzoic acid inhibited the root hair development in a dose-dependent manner, and several differential expressed genes potentially involved in hair development were identified. finally, alternative splicing events from genes were found after benzoic acid treatment. this study revealed the stress-response molecular mechanisms in p. ginseng induced by benzoic acid. rna-seqpanax ginsengautotoxicitybenzoic acidissue-copyright-statement© the author 2015 as a valuable medicinal plant, the yield of panax ginseng is seriously affected by autotoxicity, which is a common phenomenon due to continuous cropping. however, the mechanism of autotoxicity in p. ginseng is still unknown. interestingly, contigs in the ginsenoside biosynthetic pathway underwent as, providing useful information about post-transcriptional regulation in p. ginseng. compared with control, and genes were significantly up- and down-regulated across different tissues and time points, respectively.
 as a valuable medicinal plant, the yield of panax ginseng is seriously affected by autotoxicity, which is a common phenomenon due to continuous cropping. however, the mechanism of autotoxicity in p. ginseng is still unknown.  in total, high throughput sequencing of rna-seq libraries produced 100-nt reads that were assembled into contigs. compared with control, and genes were significantly up- and down-regulated across different tissues and time points, respectively. gene ontology enrichment analysis showed that ‘enzyme inhibitor activity’, ‘carboxylesterase activity’, ‘pectinesterase activity’, ‘centrosome cycle and duplication’ and ‘mitotic spindle elongation’ were enriched for the up-regulated genes. transcription factors including ap2s/erfs, mybs, and wrkys were up-regulated in roots after benzoic acid treatment. moreover, reactive oxygen species, peroxidases and superoxide dismutase contigs were up-regulated in roots after benzoic acid treatment. physiological and biochemical indexes showed that the proline and malondialdehyde content were restored to lower levels at a later stage after benzoic acid treatment. benzoic acid inhibited the root hair development in a dose-dependent manner, and several differential expressed genes potentially involved in hair development were identified. several key contigs in the flavonoid and ginsenoside biosynthesis pathways were repressed. finally, alternative splicing events from genes were found after benzoic acid treatment. interestingly, contigs in the ginsenoside biosynthetic pathway underwent as, providing useful information about post-transcriptional regulation in p. ginseng.  this study revealed the stress-response molecular mechanisms in p. ginseng induced by benzoic acid. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords rna-seqpanax ginsengautotoxicitybenzoic acidissue-copyright-statement© the author 2015 