BACKGROUND
soybean is a paleopolyploid with two duplication events, approximately  <dig> and  <dig> million years ago  <cit> .
the majority of phenylpropanoid and specialized metabolism enzymes that have been characterized belong to multi-gene families .
functional divergence within these families have been conjectured to be partially responsible for the array of plant specialized metabolites  <cit> .
in the phenylpropanoid pathway chalcone synthase  produces the chalcone nucleus that is then utilized in downstream metabolite biosynthesis.
there are at least nine chs genes in soybean  that share  <dig>  to  <dig> % sequence identity at the amino acid level  <cit> .
members of the chs gene family are differentially expressed, respond to different stimuli, and have shown functional divergence.
gmchs <dig> and gmchs <dig> are critical for isoflavonoid biosynthesis and accumulation in soybean seeds  <cit> .
these two enzymes have shown differential localization; gmchs <dig> was localized to the cytoplasm, while gmchs <dig> was localized to the nucleus and cytoplasm  <cit> .
the level of differentiation in the expression of almost identical genes in the chs family, and the putative functional specialization of their cognate proteins underlines the complexity associated with multi-gene families in large genomes such as soybean.
genetic and functional variation have also been evidenced in other closely linked enzyme families such as chalcone isomerase   <cit> , chalcone reductase   and isoflavonoid transferases  <cit> .
the characterization of multi-gene families is important in the study of overlapping branches in the phenylpropanoid pathway.
specialization can explain mechanisms that mediate competition between shared enzymes and metabolites.
further, it can describe the evolutionary path that leads to new legume or species-specific metabolites such as isoflavonoids.
the competition between the flavonoid and isoflavonoid branches of the phenylpropanoid pathway has been described in the shared usage of flavanone substrate, naringenin, and enzymes, such as chs and chi  <cit> .
the evolution of enzymatic capacity is apparent in the evolution of the chi fold from protist homologs to the legume-specific catalysis of isoliquiritigenin and naringenin chalcone to liquiritigenin and naringenin, respectively  <cit> .
isoflavonoid biosynthesis is a legume-specific branch of the diverse phenylpropanoid pathway .
this class of specialized metabolites is involved in pathogen inhibition and nitrogen-fixing symbiosis .
as part of human diet isoflavonoids are linked with a reduction in the risk of cardiovascular disease and hormone-dependent cancers .
the isoflavonoid pathway is in direct competition with the concurrent flavonoid pathways for flavanone substrates.
efforts to manipulate isoflavonoid biosynthesis or to engineer the pathway in non-legumes have underlined this effective bottleneck to metabolite accumulation .fig.
<dig> phenylalanine recruited into a diverse network of metabolism including the production of isoflavones and other specialized metabolites.
competition for the common substrate naringenin is highlighted, with downstream metabolites in yellow.
chalcone synthase  produces the chalcone, naringenin chalcone, and in legumes chs coactively with chalcone reductase  forms isoliquiritigenin chalcone.
chalcone isomerase  converts the two chalcones to their flavanone derivatives: liquiritigenin and naringenin; the latter, as stated above, can be utilized by isoflavones synthase  or a whole host of other enzymes, to produce a variety of metabolites including the isoflavones.
the multiple arrows indicate two or more steps in the pathway.
other enzymes labeled: pal, phenylalanine ammonia lyase; c4h, cinnamate 4-hydroxylase; 4cl, 4-coumarate:coa ligase; 2hid, 2-hydroxyisoflavanone dehydratase; iomt, isoflavone o-methyltransferase; f6h, flavonoid 6-hydroxylase
the present study applied a transcriptomic approach to investigate the factors that underlie isoflavonoid content variation.
the root was chosen as the organ of study; due to the dual importance of isoflavonoids in the root as a) signaling molecules for nod gene induction in nitrogen-fixing symbiotic bacteria  <cit> , and b) as phytoalexins and anti-microbial agents warding off pathogen infection  <cit> .
four soybean cultivars were chosen for this study: two with increased resistance to stem and root rot disease caused by phytophthora sojae, conrad and ac colombe; two with susceptibility to p. sojae, ac glengarry and pagoda .
multiple sources have cited increased resistance in conrad soybeans to a broad range of pathogen infection   <cit> .
here we report that the root isoflavonoid contents of conrad and ac colombe are significantly higher than that of ac glengarry and pagoda.
the variation in isoflavonoid levels in these four soybean cultivars was used as a basis for a transcriptomic approach to find the underlying mechanisms of metabolic variation.
five differentially expressed  genes encoding flavonoid and isoflavonoid metabolic enzymes were found.
this list was further refined by quantitative pcr  to two genes encoding flavonoid-metabolizing enzymes: flavonoid 3′-hydroxylase  and dihydroflavonol 4-reductase .
these genes were significantly downregulated in high isoflavonoid cultivars.
transcriptomic downregulation of genes from overlapping phenylpropanoid pathways could divert common substrates, such as naringenin, toward isoflavonoid production.
further de genes could play a role in transcriptional regulation, metabolism, transport, and plant stress response.
the lattermost could be responsible for the favorable resistance traits of conrad and ac colombe.
methods
plant material
soybean seeds from the cultivars, conrad, ac colombe, ac glengarry and pagoda, were planted in pro-mix bx mycorrhizae™ soil  in a growth chamber with a 16 h light cycle at 25 °c, and an 8 h dark cycle at 20 °c, with 60-70% relative humidity, with light intensity 100–150 μmol m−2s− <dig>
rna and isoflavonoid extraction; hplc analysis
two week-old soybean roots were harvested, flash-frozen with liquid nitrogen and ground to a fine powder.
three biological replicates were prepared for each of the four cultivars.
rna was extracted from the twelve samples according to the instructions in the rneasy mini kit .
rna integrity and concentration were recorded using the agilent  <dig> .
the same root samples were used for isoflavonoids extraction following the method described previously for seed isoflavonoid extraction and analyzed by high-performance liquid chromatography   <cit>  .
quantitative rt-pcr analysis
the same rna samples used for rna sequencing were also used in a reverse transcription reaction using the thermoscript™ rt-pcr system .
primer sequences for qpcr are listed in additional file 1: table s <dig>  all reactions were performed in three technical triplicates, and fold expression was normalized to the reference genes cons <dig> and cons <dig>  <cit> .
illumina sequencing, data quality and mapping
root tissue total rna from the four cultivars, three replicates per cultivar, was sequenced on an illumina hiseq <dig> at the dna technologies unit of the plant biotechnology institute  using 100 bp paired-end runs.
contaminating adaptor sequences were removed using a custom perl script, reads subjected to 3′ end trimming   <cit> .
reads from each library were mapped against the glycine max transcriptome v <dig> , wm <dig> a <dig> v <dig>  <cit>  using the burrows-wheeler aligner   <cit> .
pcr duplicates and reads aligned with a low mapping quality  were dropped, and uniquely mapped reads per transcript counted using samtools  <cit> .
differential gene expression analysis
counts for reads mapping to each transcript were imported into the r statistical environment for differential expression analysis with the deseq package  <cit> .
transcripts were normalized by counts across biological replicates and cultivars and the bottom 10% of low expressing transcripts dropped.
dispersion estimates were calculated for each gene using deseq functions to estimate transcript expression variation across replicates for a given cultivar  <cit> .
for the purposes of gene discovery two sets of de genes were produced based on ‘baseline’ p <  <dig>  or ‘high’ p <  <dig>  significance, representing moderate and high stringency.
p-values were adjusted using benjamini-hochberg multiple-testing correction and all diagnostic and analysis plots generated using deseq or r functions.
heatmaps were generated using the r programming package deseq, and the heatmap.
<dig> function in gplots, accessed from the cran library   <cit> .
raw gene expression data were normalized for replicate and library differences in read counts and coverage by variance stabilized transformation, using the deseq package, prior to heatmap generation.
functional annotation of differentially expressed transcripts
de genes were annotated using the soybean transcriptome annotation .
arabidopsis orthologs were used to associate ‘go’ annotations   <cit>  to soybean genes using the tair   <cit> .
over-representation analysis was performed using the ‘panther over-representation test’   <cit> .
pathway enrichment analysis was carried out using the phytomine tool ‘pathway enrichment’,   <cit> .
data from kegg  <cit>  and plantcyc , www.plantcyc.org) resources were used to conduct the analysis, with the glycine max database selected as the reference.
the benjamini-hochberg statistical test was used to adjust p values for multiple testing.
RESULTS
stem and root rot resistant soybean cultivars contain high isoflavonoid content in roots
four soybean cultivars were assayed for isoflavonoid content: two with reported resistance, conrad and ac colombe, and two with susceptibility to p. sojae, ac glengarry and pagoda.
total isoflavonoid content includes daidzein, genistein, glycitein and their corresponding glycosides as measured by hplc analysis.
as shown in fig.
<dig>  conrad and ac colombe showed higher total isoflavonoid content in the roots, as compared to ac glengarry and pagoda.fig.
<dig> root isoflavonoid content of four soybean cultivars.
total isoflavonoid content, including the six isoflavonoids, daidzein, genistein, glycitein, daidzin, genistin and glycitin, as measured by hplc.
isoflavonoids were extracted from the roots of four soybean cultivars: conrad, ac colombe, ac glengarry and pagoda.
error bars indicate sem of three independent experiments
data quality and coverage of soybean transcriptome in four soybean cultivars
the root transcriptomes of four soybean cultivars were studied, including high, conrad and ac colombe, and low, ac glengarry and pagoda, root isoflavonoid content.
root mrna from three replicates for each of the four cultivars was sequenced, using 100 bp paired-end reads.
the resulting data were then mapped against the soybean transcriptome v <dig>    <cit> .
mapping and quality information including total reads per rna sequencing library, base quality ≥  <dig> , mapped reads, uniquely mapped reads and gene coverage  are summarized in table  <dig>  coverage of the  <dig>  protein-coding loci ranged from  <dig>  to  <dig> %.
biological variation for all the gene models in each cultivar was displayed using dispersion graphs .
dispersion plots have, overall ‘tight’ dispersion, which is indicative of good quality rna sequencing libraries for the three biological replicates per cultivar.table  <dig> rna sequencing quality and coverage of soybean transcriptome
coverage of the soybean transcriptome  in soybean roots of conrad, ac colombe, ac glengarry and pagoda biological replicates ranged from  <dig>  to  <dig> %.
library replicate: the soybean cultivar replicate sequenced; total reads: the total number of sequence reads obtained; q <dig> bases : the percentage of bases with >  <dig> quality after 3′ end trimming; mapped reads: the number of reads that mapped to one of the soybean transcript models; uniquely mapped reads: number of mapped reads mapping uniquely to a transcript model with a mapping quality ≥ 20; genes hit: number of gene models with one or more reads mapped to its transcript as a percentage of the total number of gene models
differential expression analyses between high and low root isoflavonoid cultivars
alignments of rna sequencing libraries to the soybean transcriptome were imported into r and assessed for differential gene expression using deseq  <cit> .
the benjamini-hochberg adjusted p-value distribution for all genes assessed for differential expression between high and low cultivars showed a large number of significantly differentially expressed genes with a uniform dispersion .
to categorize and rank the level of expression variation among these genes, the log p-value distribution for all the genes was binned and plotted in terms of significance levels: , , and  .
an abundance of gene models being significantly, and highly significantly, differentially expressed in each of the comparisons is representative of the genetic differences between the cultivars.
four lists were generated from pair-wise differential expression analyses of low against high root isoflavonoid cultivars: ac glengarry with ac colombe; ac glengarry with conrad; pagoda with ac colombe; pagoda with conrad.
these four lists were filtered to include genes significantly differentially expressed  and are compiled in additional file 6: table s <dig>  additional file 7: table s <dig>  additional file 8: table s <dig>  and additional file 9: table s <dig>  the lists generated were assessed to find gene models that were consistently up- or down-regulated in all four differential expression analyses.
the overlap studies are illustrated in four, four-way venn diagrams for up- and down-regulation at both baseline and high p-values .
the numbers in the ovals represent the number of gene models in the four lists, and their corresponding overlap is also enumerated.
the core overlaps of  <dig>  and  <dig>  are the number of genes consistently up- and down-regulated, respectively  .
these lists of genes can be studied for their correlation with isoflavonoid biosynthesis and a putative role in accumulation.
the list of up- and down-regulated genes will hereafter be referred to as ‘high isoflavonoid’ and ‘low isoflavonoid’ genes  .
this does not denote a role in the pathway, but an association with the isoflavonoid content of the cultivars.fig.
<dig> overlap study of up- and down-regulated genes in high root isoflavonoid cultivars.
venn diagrams depicting the number of  upregulated and  downregulated genes and their overlap in four pair-wise differential expression studies: ac glengarry with ac colombe, ac glengarry with conrad, pagoda with ac colombe, pagoda with conrad.
venn diagrams are based on differential expression analysis with ‘baseline’ p <  <dig>   or ‘high’ p <  <dig>   significance.
therefore, a core of  <dig> and  <dig> genes are consistently up- and down-regulated, respectively; the higher p-value reduces this list to  <dig> and  <dig> highly up- and down-regulated genes, respectively
genes that fall outside of this core overlap would include genes up- or down-regulated in one or more of the comparisons but not all four.
these genes might represent cultivar-specific differences in the root transcriptome that are not consistently differentially expressed, within the parameters of the study.
heatmaps were generated using raw read counts for the set of  <dig> ‘high isoflavonoid’ genes, in pair-wise comparisons between low and high root isoflavonoid cultivars .
the heatmaps display the higher expression  of the ‘high isoflavonoid’ genes in the cultivars, ac colombe and conrad, as compared with ac glengarry and pagoda.
functional and structural annotation of differentially expressed genes
arabidopsis genes orthologous to the soybean differentially expressed genes were used to find corresponding ‘go’ annotations from the tair.
figure  <dig> depicts the putative molecular  and biological function , and cellular compartmentalization  of the protein products for these genes.
figures 4a, b, and c are representative of p <  <dig> , and figs.
4d, e, and f are representative of p <  <dig>  de genes.fig.
<dig> gene product annotations of upregulated and downregulated genes in high root isoflavonoid cultivars.
a, d cellular compartment; b, e biological processes; c, f molecular function annotation, using the gene ontology database annotations of tair identifiers that are homologous to differentially expressed soybean genes.
a-c are upregulated and d-f downregulated genes in high root isoflavonoid cultivars.
dark blue and light blue bars denote differentially expressed genes at a p-value of < <dig>  and < <dig> , respectively
aside from cytoplasmic proteins, the majority  of ‘high isoflavonoid’ gene products were predicted to localize to the nucleus , while the plasma membrane is the prominent cellular compartment  with ‘low isoflavonoids’ .
twelve genes were annotated as transcription factors and three as nucleic acid binding proteins, with the sum of dna-dependent transcriptional elements coming to seventeen ‘high isoflavonoid’ genes.
an aprataxin-like / bhlh protein was the only ‘low isoflavonoid’ gene product associated with transcription factor activity, putatively involved in dna metabolism or repair.
five ‘high’ and fifteen ‘low isoflavonoid’ gene products were associated with the er and golgi, as part of the secretory pathway.
among the ‘high isoflavonoid’ genes,  <dig> were annotated as being responsive to stress.
these genes included member encoding the dirigent family of proteins.
several nucleotide binding disease resistance proteins were identified, containing c-terminal leucine-rich repeat  domains fused to central nucleotide-binding  domain  .
fifteen ‘high’ and sixteen ‘low isoflavonoid’ gene products upregulated in high isoflavonoid cultivars were annotated for transferase activity.
furthermore, eight ‘high’ and nine ‘low isoflavonoid’ gene products with putative transporter activity were identified.table  <dig> differentially expressed genes that are upregulated  in high isoflavonoid cultivars: de gene models were inspected manually using annotations from the soybean database , and the ‘go’ annotations associated with their arabidopsis homologs ; they were compiled into families for consideration based on function, and potential involvement in the phenylpropanoid pathway and as underlying factors for isoflavonoid content
the log <dig> fold changes for the gene models in each of the four pair-wise comparisons between cultivars are detailed in the table
the lists of de genes were manually inspected for genes putatively involved in specialized metabolism and in particular the isoflavonoid pathway.
several members of oxidoreductase families, including fad-dependent oxidoreductases , cyps, and 2-oxoglutarate /fe-dependent oxygenases  were either up- or down-regulated .
five genes associated directly with the phenylpropanoid pathway  were identified: flavonoid 6-hydroxylase  was upregulated, while isoflavone 2′-hydroxylase , flavonoid 3′-hydroxylase , flavonoid 3′, 5′-hydroxylase , and dihydroflavonol 4-reductase  were downregulated in high root isoflavonoid cultivars.
heatmaps were also generated to display the differential expression of these five phenylpropanoid genes in the high and low isoflavonoid cultivars .table  <dig> differentially expressed genes that are downregulated  in high isoflavonoid cultivars: de gene models were inspected manually using annotations from the soybean database , and the ‘go’ annotations associated with their arabidopsis homologs ; they were compiled into families for consideration based on function, and potential involvement in the phenylpropanoid pathway and as underlying factors for isoflavonoid content
the log <dig> fold changes for the gene models in each of the four pair-wise comparisons between cultivars are detailed in the table
fig.
<dig> heatmaps of significantly differentially expressed phenylpropanoid genes.
lists of genes up- or down-regulated consistently were mined for gene models annotated with putative function in the phenylpropanoid pathway.
gmf6h was upregulated, while gmi2′h, gmf3′h, gmf3′5′h, gmdfr were downregulated in high root isoflavonoid cultivars.
data for all three replicates are shown for each soybean cultivar.
panels a and b show heatmaps of ac colombe or conrad with ac glengarry while panels c and d show heatmaps of ac colombe or conrad with pagoda
to investigate the significance of the represented go annotations in the list of genes up- and down-regulated in high isoflavonoid cultivars, panther over-representation test was performed for go-slim biological processes, molecular function and cellular component  <cit>  .
cell structure and morphogenesis elements were overrepresented in the list of upregulated genes.
genes involved in secondary metabolism, specifically those encoding membrane-bound proteins with hydrolase and oxidoreductase activity were overrepresented in the downregulated list.
therefore, a significant number of genes involved in secondary metabolism were downregulated in high isoflavonoid cultivars.
enrichment of specific metabolic pathways in the list of genes up- and down-regulated was analyzed using phytomine.
the upregulated genes were significantly enriched for phenylalanine metabolism .
downregulated genes were not significantly enriched for any single pathway; however, the classification for flavonoid genes was divided into three categories: flavonoid; flavone and flavonol; stilbenoid, diarylheptanoid and gingerol biosynthesis .
the isoflavone-metabolizing i2′h was incorrectly classified under the lattermost category , while f3′h was not included in flavonoid biosynthesis.
quantitative expression analysis to confirm differential expression candidates
to determine the exact expression levels of differentially expressed candidate genes, and to confirm the predictions of the rna sequencing results, a qrt-pcr was performed.
rna samples utilized in the making of the sequencing libraries were used for reverse-transcription.
primer sequences were designed to amplify unique regions within target genes: gmf6h, gmi2′h, gmf3′h, gmf3′5′h and gmdfr .
the expression of gmf3′h and gmdfr were significantly lower in high isoflavonoid cultivars .
there was no significant difference in the expression of gmf6h, gmi2′h and gmf3′5′h between the two sets of cultivars .fig.
<dig> expression analysis of candidate phenylpropanoid genes in four soybean cultivars.
total rna was extracted from the roots of soybean cultivars conrad, ac colombe, ac glengarry and pagoda, and was used in qpcr analysis.
expression of the following genes was analyzed: gmf6h <dig>  gmi2′h, gmf3′h, gmf3′5′h and gmdfr, using specific primers.
relative expression corresponds to mean gene expression in three biological replicates, with technical triplicates.
error bars indicate sem.
values were normalized against the reference genes gmcons <dig> and gmcons <dig>  asterisks  indicate significant differences between the samples as determined by student’s t -test
discussion
the transcriptomic exploration of underlying differences responsible for isoflavonoid content in two sets of soybean cultivars has yielded a large suite of transcripts encoding metabolic enzymes, transcriptional regulators, metabolite transporters and other potentially significant genes.
using rna sequencing, differential gene expression analysis, and subsequent qpcr analysis of candidate phenylpropanoid metabolic genes, we have identified two de genes, gmf3′h and gmdfr that could play a significant role in isoflavonoid accumulation.
the enzyme, f3′h, competes with ifs for the branch-point substrate naringenin, while dfr is involved in downstream flavonoid production .
the down-regulation of these putative genes in high isoflavonoid cultivars could describe a mechanism of diverting naringenin, from flavonoid towards isoflavonoid biosynthesis.
functional work would be required to substantiate the link between gene expression and isoflavonoid biosynthesis and/or accumulation.
differentially expressed phenylpropanoid genes
attempts at reconstituting the isoflavonoid pathway in heterologous systems have underlined the importance of competition for flavanone substrates , particularly naringenin.
this has been described as a ‘bottleneck’  <cit> ; biosynthesis of isoflavonoids increases in mutant backgrounds with reduced or absent flow of substrates into flavonoid metabolism  <cit> .
another mechanism affecting the level of isoflavone aglycones and their glucoside conjugates is conversion into downstream metabolites such as phytoalexins and signaling molecules  <cit> .
five differentially expressed genes were identified with phenylpropanoid functional annotations: gmf6h <dig>  gmi2′h, gmf3′h, gmf3′5′h, and gmdfr.
these genes and their corresponding positions in the respective pathways are displayed in fig.
<dig>  gmi2′h  is a putative isoflavonoid pathway gene annotated as a cytochrome p <dig>  family  <dig> , involved in the nadph-dependent conversion of isoflavone into 2′-hydroxyisoflavone  <cit> .
substrates for gmi2′h include the isoflavones, daidzein and genistein, and the isoflavone derivative, formononetin.
therefore, gmi2′h encodes a cyp capable of channeling isoflavones into downstream metabolites, and eventually pterocarpan phytoalexins  <cit> .
downregulation of gmi2′h in soybean roots  might lead to decreased flux of isoflavone aglycones into subsequent pathways, thereby accounting for higher levels of daidzein and genistein in conrad and ac colombe .
however, differential expression was not confirmed by qpcr analysis .
further downregulated phenylpropanoid genes identified by de analysis include gmf3′h, gmf3′5′h, and gmdfr .
all three are involved in the processing of naringenin to flavonoid metabolism.
gmf3′h  and gmf3′5′h  encode putative cyps that determine the hydroxylation of the b-ring of flavonoids .
together these cyps, f3′h and f3′5′h, are responsible for the production of dihydroflavonols, which is subsequently converted to leucoanthocyanidin by the action of dfr  <cit> , another flavonoid gene downregulated in high isoflavonoid cultivars .
the combined downregulation of these flavonoid pathway genes could play a role in the flux of flavanone substrates  into the competing isoflavonoid biosynthesis branch and the increased production of genistein.
gmf3′h and gmdfr were found, by qpcr analysis, to be significantly downregulated in high isoflavonoid cultivars; gmf3′5′h was reduced, but not significantly.
the majority of isoflavonoid content in soybean is comprised of daidzein and genistein and their β-glycoside/malonyl derivatives; however, glycitein and its derivatives form the third, and smallest component of this metabolite class.
flavonoid 6-hydroxylase  <dig>   transcript was found at a higher level in high isoflavonoid cultivars  , and encodes for a cyp enzyme that catalyzes the a-ring hydroxylation of liquiritigenin, synthesizing  <dig> ,4′-trihydroxyflavanone  <cit> .
gmf6h <dig> is associated with two other highly similar isoforms gmf6h <dig> and gmf6h <dig> of the cyp71d <dig> family  <cit> .
the intermediate produced in this reaction is further catalyzed by ifs  to produce  <dig> , <dig> ′-trihydroxyflavanone, which undergoes a dehydration , and methylation  event to produce glycitein  <cit> .
it has been reported that gmf6h <dig> is the isoform responsible for glycitein content of soybean seeds, which is almost-exclusively limited to the hypocotyls , as it was the sole member expressed in the tissue of question  <cit> .
expression of gmf6h <dig> was concurrent with isoflavonoid accumulation in the hypocotyl ; it was also absent from the seeds of soybean cultivars that are null-mutants for glycitein accumulation.
strikingly, gmf6h <dig> unlike the other members of the gmf6h family lacks the characteristic n-terminal hydrophobic transmembrane domain that anchors cyp enzymes into the er  <cit> .
gmf6h <dig> gene expression, as analyzed by qpcr , did not correlate with high isoflavonoid accumulation .
the discrepancy between qpcr and differential expression analysis results for gmf6h <dig>  gmi2′h and gmf3′5′h could indicate that the de results were not accurate to the same level as qpcr analysis, as these genes fell in the ‘baseline’ p-value score  .
alternatively, the gene expression and function of a different isoform of the f6h or i2′h gene families could correlate with high isoflavonoid content, as assigning reads to highly similar genetic regions is problematic.
enzymatic candidates for a role in specialized metabolism
table  <dig> and table  <dig> are a compilation of manually selected candidates from the high stringency  de analysis that might be associated with the isoflavonoid content of the four cultivars in question.
in addition to gmdfr that is significantly downregulated, there are several cyps, 2-odds and fad/nadp-dependent oxidoreductases in the list that could be top candidates for involvement in specialized metabolism.
overrepresentation and pathway enrichment statistics  show that genes involved in secondary metabolism, particularly those encoding hydrolase and oxidoreductase activity are significantly enriched in the ‘low isoflavonoid’ gene list.
membrane integral proteins are also overrepresented within this list, which could impact the formation of isoflavonoid or flavonoid multi-enzyme complexes or ‘metabolons’ .
2-odds are a very versatile family of dioxygenases that are superior to cyps in the range of reactions they catalyze, from hydroxylations to ring fragmentation  <cit> .
they are ubiquitously distributed throughout nature, and have very important roles in core and specialized metabolic pathways  <cit> , including histone demethylation  <cit> , and the biosynthesis of gibberellic acid  <cit> , flavonoids  <cit> , benzylisoquinoline alkaloids  <cit> , and glucosinalates  <cit> , to name a few.
it is important to note that flavanone 3β-hydroxylase  is a 2-odd, converting a 2s-flavanone to dihydroflavonols by a c- <dig> hydroxylation  <cit> .
this is an enzyme that competes with isoflavonoid biosynthesis, and downregulation of a corresponding gene could underlie increased isoflavonoid content.
there are several reactions in the flavonoid and isoflavonoid pathways that have yet to be elucidated and table  <dig> and table  <dig> could provide a starting point towards gene discovery.
among them there are four putative 2-odd genes, two that are upregulated: glyma.10g <dig> .p and glyma.17g <dig> .p, and two downregulated members: glyma.08g <dig> .p and glyma.14g <dig> .p.
other metabolic enzymes represented in the list of high stringency de genes include a downregulated homolog of cyp71b <dig> , and several fad/nad-binding oxidoreductases  .
the latter group has been less present within the current literature on phenylpropanoid pathway; however an nadph/fad-dependent enzyme was shown to be responsible in the additional hydroxylation of positions  <dig> and/or  <dig> of the flavonol a-ring  <cit>  this provides a precedence for fad/nad-linked oxidoreductases being involved in the content and composition of phenylpropanoids, and should be further researched.
stress response in high isoflavonoid cultivars
isoflavonoid biosynthesis is inducible by several abiotic and biotic stimuli.
therefore, it would be interesting to further investigate the  <dig> ‘high isoflavonoid’ genes that are annotated as being responsive to stress .
subsets of these genes are annotated as being involved in disease resistance such as, dirigent proteins, which dictate the stereochemistry of other proteins  <cit> , and several nucleotide binding disease resistance proteins.
the latter included nb-lrr proteins that are associated with r gene function  <cit> .
several of these genes were annotated as being root hair specific.
as root hairs are often the site for pathogen entry, the convergence of tissue-specialization and function could be indicative of a role for these genes in inhibiting pathogen infection or colonization.
as mentioned earlier, conrad and ac colombe are also favored for their increased resistance to p. sojae  and other soybean pathogens  <cit> .
therefore, the differential expression of genes involved in stress response, localized to the root hairs could be suggestive of a mechanism underlying the improved resistance to pathogen infection.
transcriptional regulation: transcription factors and chromatin regulators
subcellular localization of the ‘high isoflavonoid’ genes to the nucleus highlights the potential role of the cognate proteins in transcriptional regulation.
the list of  <dig> nuclear-localized genes includes  <dig> transcriptional regulators , including members of the bhlh superfamily of proteins and myb transcription factors.
physical interaction and regulatory synergy has been shown between these family of proteins in regulating gene expression  <cit> .
combinatorial plant gene regulation might be a factor in the coordinated redirection of flux in specialized metabolism, leading to isoflavonoid content variation in soybean cultivars.
a gras family transcription factor  was significantly upregulated  in high isoflavonoid cultivars .
based on homology to arabidopsis transcription factors it is putatively denoted for stress transduction.
as isoflavonoid production is closely related to stress response, this transcription factor could be an important candidate for further study.
another function for nuclear-localized proteins could be the regulation of chromatin structure, which thereby effects broad transcriptomic changes.
two genes upregulated  and annotated for such function were identified: a set7/ <dig> family, histone h3k4-specific methyltransferase protein and a histone deacetylase  <dig>  family protein .
the former, a histone-modifier protein, cannot be immediately associated with activation or repression of transcriptional regulation, as methylation of h3k <dig> would have to be regarded in the larger landscape and context of methylations  <cit> .
gmhda <dig>  on the other hand, could be associated with the downregulation of certain genes, explaining the transcriptomic changes between cultivars.
deacetylation conferred by elements such as gmhda <dig> can increase dna-histone affinity and condense the overall structure, rendering it inaccessible to transcription machinery  <cit> .
hypoacyetylated chromatin is silent, or has reduced gene expression.
therefore, gmhda <dig> might be a histone-modifying element responsible for some of the overall changes in the root transcriptome that coincide between the high isoflavonoid cultivars, conrad and ac colombe.
vacuolar sequestration: a putative regulatory mechanism for isoflavonoid accumulation
isoflavone aglycones produced by an enzyme complex tethered to the er, are conjugated by transferases, allowing their appropriate channeling into the vacuole  <cit> .
from this storage point, isoflavonoids can be sequestered outside of the metabolic flux and released upon demand for stress response or as signaling molecules for symbiosis with nitrogen-fixing rhizobia  <cit> .
interestingly, there were  <dig> up- and  <dig> down-regulated genes annotated for transferase activity;  <dig> up- and  <dig> down-regulated genes with putative transporter activity; the suite of transferases and transporters being simultaneously regulated could denote a preference towards the conjugation and storage of a certain class of metabolites .
however, there is still ambiguity over the specificity of vacuolar transporters involved in the channeling of specialized metabolites, and broad substrate acceptance within classes of compounds has been reported .
interestingly, one of the candidate upregulated genes identified in the de analysis of soybean cultivars, was an uncharacterized member of the mate family , and hereafter referred to as gmmate <dig>  mate transporters are important families of proteins in the glycosylation and the subsequent transport of isoflavones  <cit> .
gmmate <dig> was consistently upregulated in high isoflavonoid cultivars, suggesting a possible role in the sequestration of isoflavonoids.
the predicted localization of this protein was to the secretory pathway, based on the presence of signal peptide.
this could imply localization to the er, golgi bodies or the vacuole, with the last option being highly probable, given its functional annotation as a mate transporter.
the deduced amino acid sequence of this gene had an amino acid sequence similarity of  <dig> % with the tt <dig> homolog in arabidopsis, based on clustal w  <cit> , the latter being a known proanthocyanidin vacuolar transporter  <cit> .
mate transporters in medicago truncatula, mate <dig> and mate <dig>  are involved in transport of phenylpropanoids.
mate <dig> is a tonoplast epicatechin 3′-o-glucoside  transporter  <cit> .
mate <dig> is a flavonoid transporter involved in vacuolar sequestration of anthocyanins and other flavonoids in flowers and leaves.
mate <dig> transporter prefers malonylated flavonoid glucosides and is co-expressed with three genes encoding malonyltransferases  <cit> .
the discovery of an isoflavonoid-specific mate transporter, in turn, could be used as a tool to investigate co-expressed transferases, depicting a more complete image of the conjugation and transport of isoflavonoids.
furthermore, up-regulation of gmmate <dig> could indicate a mechanism for the vacuolar transport and sequestration of isoflavonoids and/or precursors, in competition with parallel pathways for metabolic flux.
CONCLUSIONS
the coordinated regulation of genes encoding enzymes, transporters, transcription factors and other molecular elements could lead to the increase of root isoflavonoids in conrad and ac colombe.
these transcriptomic elements can help overcome bottlenecks in isoflavonoid production for soybean cultivars, legumes, and heterologous systems.
genes annotated for disease and stress response could be important in conferring soybean resistance to p. sojae infection in conrad and ac colombe.
functional genomics should be employed to characterize the role of such transcriptomic elements in isoflavonoid production and disease resistance.
additional files
additional file 1:
table s <dig>  sequence of oligonucleotides used for qpcr of soybean genes in high and low root isoflavonoid cultivars to confirm differential expression analysis from rnaseq.
additional file 2:
figure s <dig>  estimated intra-cultivar dispersion for all genes.
dispersion values plotted as a function of average normalized expression for each gene model.
biological variation displayed for  conrad  ac colombe  ac glengarry and  pagoda in three biological replicates.
additional file 3:
figure s <dig>  m-a plot of log <dig> fold change as a function of average expression of gene.
significantly differentially expressed genes  are represented in red.
plots displayed for pair-wise comparisons  ac glengarry with ac colombe  ac glengarry with conrad  pagoda with ac colombe  pagoda with conrad.
additional file 4:
figure s <dig>  gene frequency as a function of benjamini-hochberg adjusted p-values.
histogram of p-values for all genes assessed for significant differential expression between two cultivars  ac glengarry with ac colombe  ac glengarry with conrad  pagoda with ac colombe  pagoda with conrad.
additional file 5:
figure s <dig>  genes ranked by total read count against –log <dig> of the p-value.
significance scores divided into three levels: , , and .
average gene counts divided into percentiles.
four plots representing significance of differential expression between cultivars  ac glengarry with ac colombe  ac glengarry with conrad  pagoda with ac colombe  pagoda with conrad.
additional file 6:
table s <dig>  list of genes upregulated in high  as compared with low  root isoflavonoid content cultivars.
differentially expressed genes  in the four comparisons between high and low cultivars were analyzed for overlap , generating a set of  <dig> candidates upregulated consistently in high root isoflavonoid cultivars.
additional file 7:
table s <dig>  list of genes downregulated in high  as compared with low  root isoflavonoid content cultivars.
differentially expressed genes  in the four comparisons between high and low cultivars were analyzed for overlap , generating a set of  <dig> candidates downregulated consistently in high root isoflavonoid cultivars.
additional file 8:
table s <dig>  list of genes upregulated in high  as compared with low  root isoflavonoid content cultivars.
highly differentially expressed genes  in the four comparisons between high and low cultivars were analyzed for overlap , generating a set of  <dig> candidates upregulated consistently in high root isoflavonoid cultivars.
additional file 9:
table s <dig>  list of genes downregulated in high  as compared with low  root isoflavonoid content cultivars.
list of genes downregulated in high  as compared with low  root isoflavonoid content cultivars.
highly differentially expressed genes  in the four comparisons between high and low cultivars were analyzed for overlap , generating a set of  <dig> candidates downregulated consistently in high root isoflavonoid cultivars.
these genes were annotated using the soybean database and have been compiled below.
additional file 10:
figure s <dig>  heatmaps of genes upregulated in high root isoflavonoid cultivars.
read counts for gene models were normalized across cultivars.
four heatmaps representing genes upregulated in high isoflavonoid cultivars as compared with low root isoflavonoid cultivars were generated :  ac glengarry with ac colombe  ac glengarry with conrad  pagoda with ac colombe  pagoda with conrad.
additional file 11:
table s <dig>  genes upregulated in ‘high isoflavonoid’ cultivars  were analyzed for overrepresentation of panther go-slim classifications for: biological process, molecular function and cellular component.
additional file 12:
table s <dig>  genes downregulated in ‘high isoflavonoid’ cultivars  were analyzed for overrepresentation of panther go-slim classifications for: biological process, molecular function and cellular component.
additional file 13:
table s <dig>  genes upregulated in ‘high isoflavonoid’ cultivars  were analyzed for pathway enrichment.
additional file 14:
table s <dig>  genes downregulated in ‘high isoflavonoid’ cultivars  were analyzed for pathway enrichment.
additional file 15:
table s <dig>  genes annotated with disease or stress response and upregulated in ‘high isoflavonoid’ cultivars.
additional file 16:
table s <dig>  genes annotated for transcriptional regulation and upregulated in ‘high isoflavonoid’ cultivars.
abbreviations
chichalcone isomerase
chrchalcone reductase
chschalcone synthase
cypcytochrome p450
dedifferentially expressed
gmdfrsoybean dihydroflavonol 4-reductase
gmf3′hsoybean flavonoid 3′-hydroxylse
gogene ontology

 soybean is a paleopolyploid with two duplication events, approximately and million years ago. the majority of phenylpropanoid and specialized metabolism enzymes that have been characterized belong to multi-gene families. functional divergence within these families have been conjectured to be partially responsible for the array of plant specialized metabolites. in the phenylpropanoid pathway chalcone synthase produces the chalcone nucleus that is then utilized in downstream metabolite biosynthesis. there are at least nine chs genes in soybean that share to % sequence identity at the amino acid level. members of the chs gene family are differentially expressed, respond to different stimuli, and have shown functional divergence. gmchs and gmchs are critical for isoflavonoid biosynthesis and accumulation in soybean seeds. these two enzymes have shown differential localization; gmchs was localized to the cytoplasm, while gmchs was localized to the nucleus and cytoplasm. the level of differentiation in the expression of almost identical genes in the chs family, and the putative functional specialization of their cognate proteins underlines the complexity associated with multi-gene families in large genomes such as soybean. genetic and functional variation have also been evidenced in other closely linked enzyme families such as chalcone isomerase, chalcone reductase and isoflavonoid transferases. the characterization of multi-gene families is important in the study of overlapping branches in the phenylpropanoid pathway. specialization can explain mechanisms that mediate competition between shared enzymes and metabolites. further, it can describe the evolutionary path that leads to new legume or species-specific metabolites such as isoflavonoids. the competition between the flavonoid and isoflavonoid branches of the phenylpropanoid pathway has been described in the shared usage of flavanone substrate, naringenin, and enzymes, such as chs and chi. the evolution of enzymatic capacity is apparent in the evolution of the chi fold from protist homologs to the legume-specific catalysis of isoliquiritigenin and naringenin chalcone to liquiritigenin and naringenin, respectively. isoflavonoid biosynthesis is a legume-specific branch of the diverse phenylpropanoid pathway. this class of specialized metabolites is involved in pathogen inhibition and nitrogen-fixing symbiosis. as part of human diet isoflavonoids are linked with a reduction in the risk of cardiovascular disease and hormone-dependent cancers. the isoflavonoid pathway is in direct competition with the concurrent flavonoid pathways for flavanone substrates. efforts to manipulate isoflavonoid biosynthesis or to engineer the pathway in non-legumes have underlined this effective bottleneck to metabolite accumulation.fig. phenylalanine recruited into a diverse network of metabolism including the production of isoflavones and other specialized metabolites. competition for the common substrate naringenin is highlighted, with downstream metabolites in yellow. chalcone synthase produces the chalcone, naringenin chalcone, and in legumes chs coactively with chalcone reductase forms isoliquiritigenin chalcone. chalcone isomerase converts the two chalcones to their flavanone derivatives: liquiritigenin and naringenin; the latter, as stated above, can be utilized by isoflavones synthase or a whole host of other enzymes, to produce a variety of metabolites including the isoflavones. the multiple arrows indicate two or more steps in the pathway. other enzymes labeled: pal, phenylalanine ammonia lyase; c4h, cinnamate 4-hydroxylase; 4cl, 4-coumarate:coa ligase; 2hid, 2-hydroxyisoflavanone dehydratase; iomt, isoflavone o-methyltransferase; f6h, flavonoid 6-hydroxylase the present study applied a transcriptomic approach to investigate the factors that underlie isoflavonoid content variation. the root was chosen as the organ of study; due to the dual importance of isoflavonoids in the root as a) signaling molecules for nod gene induction in nitrogen-fixing symbiotic bacteria, and b) as phytoalexins and anti-microbial agents warding off pathogen infection. four soybean cultivars were chosen for this study: two with increased resistance to stem and root rot disease caused by phytophthora sojae, conrad and ac colombe; two with susceptibility to p. sojae, ac glengarry and pagoda. multiple sources have cited increased resistance in conrad soybeans to a broad range of pathogen infection. here we report that the root isoflavonoid contents of conrad and ac colombe are significantly higher than that of ac glengarry and pagoda. the variation in isoflavonoid levels in these four soybean cultivars was used as a basis for a transcriptomic approach to find the underlying mechanisms of metabolic variation. five differentially expressed genes encoding flavonoid and isoflavonoid metabolic enzymes were found. this list was further refined by quantitative pcr to two genes encoding flavonoid-metabolizing enzymes: flavonoid 3′-hydroxylase and dihydroflavonol 4-reductase. these genes were significantly downregulated in high isoflavonoid cultivars. transcriptomic downregulation of genes from overlapping phenylpropanoid pathways could divert common substrates, such as naringenin, toward isoflavonoid production. further de genes could play a role in transcriptional regulation, metabolism, transport, and plant stress response. the lattermost could be responsible for the favorable resistance traits of conrad and ac colombe. methods plant material soybean seeds from the cultivars, conrad, ac colombe, ac glengarry and pagoda, were planted in pro-mix bx mycorrhizae™ soil in a growth chamber with a 16 h light cycle at 25 °c, and an 8 h dark cycle at 20 °c, with 60-70% relative humidity, with light intensity 100–150 μmol m−2s− rna and isoflavonoid extraction; hplc analysis two week-old soybean roots were harvested, flash-frozen with liquid nitrogen and ground to a fine powder. three biological replicates were prepared for each of the four cultivars. rna was extracted from the twelve samples according to the instructions in the rneasy mini kit. rna integrity and concentration were recorded using the agilent. the same root samples were used for isoflavonoids extraction following the method described previously for seed isoflavonoid extraction and analyzed by high-performance liquid chromatography. quantitative rt-pcr analysis the same rna samples used for rna sequencing were also used in a reverse transcription reaction using the thermoscript™ rt-pcr system. primer sequences for qpcr are listed in additional file 1: table s all reactions were performed in three technical triplicates, and fold expression was normalized to the reference genes cons and cons. illumina sequencing, data quality and mapping root tissue total rna from the four cultivars, three replicates per cultivar, was sequenced on an illumina hiseq at the dna technologies unit of the plant biotechnology institute using 100 bp paired-end runs. contaminating adaptor sequences were removed using a custom perl script, reads subjected to 3′ end trimming. reads from each library were mapped against the glycine max transcriptome v, wm a v using the burrows-wheeler aligner. pcr duplicates and reads aligned with a low mapping quality were dropped, and uniquely mapped reads per transcript counted using samtools. differential gene expression analysis counts for reads mapping to each transcript were imported into the r statistical environment for differential expression analysis with the deseq package. transcripts were normalized by counts across biological replicates and cultivars and the bottom 10% of low expressing transcripts dropped. dispersion estimates were calculated for each gene using deseq functions to estimate transcript expression variation across replicates for a given cultivar. for the purposes of gene discovery two sets of de genes were produced based on ‘baseline’ p <  or ‘high’ p <  significance, representing moderate and high stringency. p-values were adjusted using benjamini-hochberg multiple-testing correction and all diagnostic and analysis plots generated using deseq or r functions. heatmaps were generated using the r programming package deseq, and the heatmap. function in gplots, accessed from the cran library. raw gene expression data were normalized for replicate and library differences in read counts and coverage by variance stabilized transformation, using the deseq package, prior to heatmap generation. functional annotation of differentially expressed transcripts de genes were annotated using the soybean transcriptome annotation. arabidopsis orthologs were used to associate ‘go’ annotations to soybean genes using the tair. over-representation analysis was performed using the ‘panther over-representation test’. pathway enrichment analysis was carried out using the phytomine tool ‘pathway enrichment’,. data from kegg and plantcyc, www.plantcyc.org) resources were used to conduct the analysis, with the glycine max database selected as the reference. the benjamini-hochberg statistical test was used to adjust p values for multiple testing.  stem and root rot resistant soybean cultivars contain high isoflavonoid content in roots four soybean cultivars were assayed for isoflavonoid content: two with reported resistance, conrad and ac colombe, and two with susceptibility to p. sojae, ac glengarry and pagoda. total isoflavonoid content includes daidzein, genistein, glycitein and their corresponding glycosides as measured by hplc analysis. as shown in fig.  conrad and ac colombe showed higher total isoflavonoid content in the roots, as compared to ac glengarry and pagoda.fig. root isoflavonoid content of four soybean cultivars. total isoflavonoid content, including the six isoflavonoids, daidzein, genistein, glycitein, daidzin, genistin and glycitin, as measured by hplc. isoflavonoids were extracted from the roots of four soybean cultivars: conrad, ac colombe, ac glengarry and pagoda. error bars indicate sem of three independent experiments data quality and coverage of soybean transcriptome in four soybean cultivars the root transcriptomes of four soybean cultivars were studied, including high, conrad and ac colombe, and low, ac glengarry and pagoda, root isoflavonoid content. root mrna from three replicates for each of the four cultivars was sequenced, using 100 bp paired-end reads. the resulting data were then mapped against the soybean transcriptome v. mapping and quality information including total reads per rna sequencing library, base quality ≥ , mapped reads, uniquely mapped reads and gene coverage are summarized in table  coverage of the protein-coding loci ranged from to %. biological variation for all the gene models in each cultivar was displayed using dispersion graphs. dispersion plots have, overall ‘tight’ dispersion, which is indicative of good quality rna sequencing libraries for the three biological replicates per cultivar.table rna sequencing quality and coverage of soybean transcriptome coverage of the soybean transcriptome in soybean roots of conrad, ac colombe, ac glengarry and pagoda biological replicates ranged from to %. library replicate: the soybean cultivar replicate sequenced; total reads: the total number of sequence reads obtained; q bases: the percentage of bases with >  quality after 3′ end trimming; mapped reads: the number of reads that mapped to one of the soybean transcript models; uniquely mapped reads: number of mapped reads mapping uniquely to a transcript model with a mapping quality ≥ 20; genes hit: number of gene models with one or more reads mapped to its transcript as a percentage of the total number of gene models differential expression analyses between high and low root isoflavonoid cultivars alignments of rna sequencing libraries to the soybean transcriptome were imported into r and assessed for differential gene expression using deseq. the benjamini-hochberg adjusted p-value distribution for all genes assessed for differential expression between high and low cultivars showed a large number of significantly differentially expressed genes with a uniform dispersion. to categorize and rank the level of expression variation among these genes, the log p-value distribution for all the genes was binned and plotted in terms of significance levels:,, and. an abundance of gene models being significantly, and highly significantly, differentially expressed in each of the comparisons is representative of the genetic differences between the cultivars. four lists were generated from pair-wise differential expression analyses of low against high root isoflavonoid cultivars: ac glengarry with ac colombe; ac glengarry with conrad; pagoda with ac colombe; pagoda with conrad. these four lists were filtered to include genes significantly differentially expressed and are compiled in additional file 6: table s additional file 7: table s additional file 8: table s and additional file 9: table s the lists generated were assessed to find gene models that were consistently up- or down-regulated in all four differential expression analyses. the overlap studies are illustrated in four, four-way venn diagrams for up- and down-regulation at both baseline and high p-values. the numbers in the ovals represent the number of gene models in the four lists, and their corresponding overlap is also enumerated. the core overlaps of and are the number of genes consistently up- and down-regulated, respectively. these lists of genes can be studied for their correlation with isoflavonoid biosynthesis and a putative role in accumulation. the list of up- and down-regulated genes will hereafter be referred to as ‘high isoflavonoid’ and ‘low isoflavonoid’ genes. this does not denote a role in the pathway, but an association with the isoflavonoid content of the cultivars.fig. overlap study of up- and down-regulated genes in high root isoflavonoid cultivars. venn diagrams depicting the number of upregulated and downregulated genes and their overlap in four pair-wise differential expression studies: ac glengarry with ac colombe, ac glengarry with conrad, pagoda with ac colombe, pagoda with conrad. venn diagrams are based on differential expression analysis with ‘baseline’ p <  or ‘high’ p <  significance. therefore, a core of and genes are consistently up- and down-regulated, respectively; the higher p-value reduces this list to and highly up- and down-regulated genes, respectively genes that fall outside of this core overlap would include genes up- or down-regulated in one or more of the comparisons but not all four. these genes might represent cultivar-specific differences in the root transcriptome that are not consistently differentially expressed, within the parameters of the study. heatmaps were generated using raw read counts for the set of ‘high isoflavonoid’ genes, in pair-wise comparisons between low and high root isoflavonoid cultivars. the heatmaps display the higher expression of the ‘high isoflavonoid’ genes in the cultivars, ac colombe and conrad, as compared with ac glengarry and pagoda. functional and structural annotation of differentially expressed genes arabidopsis genes orthologous to the soybean differentially expressed genes were used to find corresponding ‘go’ annotations from the tair. figure  depicts the putative molecular and biological function, and cellular compartmentalization of the protein products for these genes. figures 4a, b, and c are representative of p < , and figs. 4d, e, and f are representative of p <  de genes.fig. gene product annotations of upregulated and downregulated genes in high root isoflavonoid cultivars. a, d cellular compartment; b, e biological processes; c, f molecular function annotation, using the gene ontology database annotations of tair identifiers that are homologous to differentially expressed soybean genes. a-c are upregulated and d-f downregulated genes in high root isoflavonoid cultivars. dark blue and light blue bars denote differentially expressed genes at a p-value of < and <, respectively aside from cytoplasmic proteins, the majority of ‘high isoflavonoid’ gene products were predicted to localize to the nucleus, while the plasma membrane is the prominent cellular compartment with ‘low isoflavonoids’. twelve genes were annotated as transcription factors and three as nucleic acid binding proteins, with the sum of dna-dependent transcriptional elements coming to seventeen ‘high isoflavonoid’ genes. an aprataxin-like / bhlh protein was the only ‘low isoflavonoid’ gene product associated with transcription factor activity, putatively involved in dna metabolism or repair. five ‘high’ and fifteen ‘low isoflavonoid’ gene products were associated with the er and golgi, as part of the secretory pathway. among the ‘high isoflavonoid’ genes, were annotated as being responsive to stress. these genes included member encoding the dirigent family of proteins. several nucleotide binding disease resistance proteins were identified, containing c-terminal leucine-rich repeat domains fused to central nucleotide-binding domain. fifteen ‘high’ and sixteen ‘low isoflavonoid’ gene products upregulated in high isoflavonoid cultivars were annotated for transferase activity. furthermore, eight ‘high’ and nine ‘low isoflavonoid’ gene products with putative transporter activity were identified.table differentially expressed genes that are upregulated in high isoflavonoid cultivars: de gene models were inspected manually using annotations from the soybean database, and the ‘go’ annotations associated with their arabidopsis homologs; they were compiled into families for consideration based on function, and potential involvement in the phenylpropanoid pathway and as underlying factors for isoflavonoid content the log fold changes for the gene models in each of the four pair-wise comparisons between cultivars are detailed in the table the lists of de genes were manually inspected for genes putatively involved in specialized metabolism and in particular the isoflavonoid pathway. several members of oxidoreductase families, including fad-dependent oxidoreductases, cyps, and 2-oxoglutarate /fe-dependent oxygenases were either up- or down-regulated. five genes associated directly with the phenylpropanoid pathway were identified: flavonoid 6-hydroxylase was upregulated, while isoflavone 2′-hydroxylase, flavonoid 3′-hydroxylase, flavonoid 3′, 5′-hydroxylase, and dihydroflavonol 4-reductase were downregulated in high root isoflavonoid cultivars. heatmaps were also generated to display the differential expression of these five phenylpropanoid genes in the high and low isoflavonoid cultivars.table differentially expressed genes that are downregulated in high isoflavonoid cultivars: de gene models were inspected manually using annotations from the soybean database, and the ‘go’ annotations associated with their arabidopsis homologs; they were compiled into families for consideration based on function, and potential involvement in the phenylpropanoid pathway and as underlying factors for isoflavonoid content the log fold changes for the gene models in each of the four pair-wise comparisons between cultivars are detailed in the table fig. heatmaps of significantly differentially expressed phenylpropanoid genes. lists of genes up- or down-regulated consistently were mined for gene models annotated with putative function in the phenylpropanoid pathway. gmf6h was upregulated, while gmi2′h, gmf3′h, gmf3′5′h, gmdfr were downregulated in high root isoflavonoid cultivars. data for all three replicates are shown for each soybean cultivar. panels a and b show heatmaps of ac colombe or conrad with ac glengarry while panels c and d show heatmaps of ac colombe or conrad with pagoda to investigate the significance of the represented go annotations in the list of genes up- and down-regulated in high isoflavonoid cultivars, panther over-representation test was performed for go-slim biological processes, molecular function and cellular component. cell structure and morphogenesis elements were overrepresented in the list of upregulated genes. genes involved in secondary metabolism, specifically those encoding membrane-bound proteins with hydrolase and oxidoreductase activity were overrepresented in the downregulated list. therefore, a significant number of genes involved in secondary metabolism were downregulated in high isoflavonoid cultivars. enrichment of specific metabolic pathways in the list of genes up- and down-regulated was analyzed using phytomine. the upregulated genes were significantly enriched for phenylalanine metabolism. downregulated genes were not significantly enriched for any single pathway; however, the classification for flavonoid genes was divided into three categories: flavonoid; flavone and flavonol; stilbenoid, diarylheptanoid and gingerol biosynthesis. the isoflavone-metabolizing i2′h was incorrectly classified under the lattermost category, while f3′h was not included in flavonoid biosynthesis. quantitative expression analysis to confirm differential expression candidates to determine the exact expression levels of differentially expressed candidate genes, and to confirm the predictions of the rna sequencing results, a qrt-pcr was performed. rna samples utilized in the making of the sequencing libraries were used for reverse-transcription. primer sequences were designed to amplify unique regions within target genes: gmf6h, gmi2′h, gmf3′h, gmf3′5′h and gmdfr. the expression of gmf3′h and gmdfr were significantly lower in high isoflavonoid cultivars. there was no significant difference in the expression of gmf6h, gmi2′h and gmf3′5′h between the two sets of cultivars.fig. expression analysis of candidate phenylpropanoid genes in four soybean cultivars. total rna was extracted from the roots of soybean cultivars conrad, ac colombe, ac glengarry and pagoda, and was used in qpcr analysis. expression of the following genes was analyzed: gmf6h gmi2′h, gmf3′h, gmf3′5′h and gmdfr, using specific primers. relative expression corresponds to mean gene expression in three biological replicates, with technical triplicates. error bars indicate sem. values were normalized against the reference genes gmcons and gmcons asterisks indicate significant differences between the samples as determined by student’s t -test discussion the transcriptomic exploration of underlying differences responsible for isoflavonoid content in two sets of soybean cultivars has yielded a large suite of transcripts encoding metabolic enzymes, transcriptional regulators, metabolite transporters and other potentially significant genes. using rna sequencing, differential gene expression analysis, and subsequent qpcr analysis of candidate phenylpropanoid metabolic genes, we have identified two de genes, gmf3′h and gmdfr that could play a significant role in isoflavonoid accumulation. the enzyme, f3′h, competes with ifs for the branch-point substrate naringenin, while dfr is involved in downstream flavonoid production. the down-regulation of these putative genes in high isoflavonoid cultivars could describe a mechanism of diverting naringenin, from flavonoid towards isoflavonoid biosynthesis. functional work would be required to substantiate the link between gene expression and isoflavonoid biosynthesis and/or accumulation. differentially expressed phenylpropanoid genes attempts at reconstituting the isoflavonoid pathway in heterologous systems have underlined the importance of competition for flavanone substrates, particularly naringenin. this has been described as a ‘bottleneck’; biosynthesis of isoflavonoids increases in mutant backgrounds with reduced or absent flow of substrates into flavonoid metabolism. another mechanism affecting the level of isoflavone aglycones and their glucoside conjugates is conversion into downstream metabolites such as phytoalexins and signaling molecules. five differentially expressed genes were identified with phenylpropanoid functional annotations: gmf6h gmi2′h, gmf3′h, gmf3′5′h, and gmdfr. these genes and their corresponding positions in the respective pathways are displayed in fig.  gmi2′h is a putative isoflavonoid pathway gene annotated as a cytochrome p family, involved in the nadph-dependent conversion of isoflavone into 2′-hydroxyisoflavone. substrates for gmi2′h include the isoflavones, daidzein and genistein, and the isoflavone derivative, formononetin. therefore, gmi2′h encodes a cyp capable of channeling isoflavones into downstream metabolites, and eventually pterocarpan phytoalexins. downregulation of gmi2′h in soybean roots might lead to decreased flux of isoflavone aglycones into subsequent pathways, thereby accounting for higher levels of daidzein and genistein in conrad and ac colombe. however, differential expression was not confirmed by qpcr analysis. further downregulated phenylpropanoid genes identified by de analysis include gmf3′h, gmf3′5′h, and gmdfr. all three are involved in the processing of naringenin to flavonoid metabolism. gmf3′h and gmf3′5′h encode putative cyps that determine the hydroxylation of the b-ring of flavonoids. together these cyps, f3′h and f3′5′h, are responsible for the production of dihydroflavonols, which is subsequently converted to leucoanthocyanidin by the action of dfr, another flavonoid gene downregulated in high isoflavonoid cultivars. the combined downregulation of these flavonoid pathway genes could play a role in the flux of flavanone substrates into the competing isoflavonoid biosynthesis branch and the increased production of genistein. gmf3′h and gmdfr were found, by qpcr analysis, to be significantly downregulated in high isoflavonoid cultivars; gmf3′5′h was reduced, but not significantly. the majority of isoflavonoid content in soybean is comprised of daidzein and genistein and their β-glycoside/malonyl derivatives; however, glycitein and its derivatives form the third, and smallest component of this metabolite class. flavonoid 6-hydroxylase transcript was found at a higher level in high isoflavonoid cultivars, and encodes for a cyp enzyme that catalyzes the a-ring hydroxylation of liquiritigenin, synthesizing,4′-trihydroxyflavanone. gmf6h is associated with two other highly similar isoforms gmf6h and gmf6h of the cyp71d family. the intermediate produced in this reaction is further catalyzed by ifs to produce, ′-trihydroxyflavanone, which undergoes a dehydration, and methylation event to produce glycitein. it has been reported that gmf6h is the isoform responsible for glycitein content of soybean seeds, which is almost-exclusively limited to the hypocotyls, as it was the sole member expressed in the tissue of question. expression of gmf6h was concurrent with isoflavonoid accumulation in the hypocotyl; it was also absent from the seeds of soybean cultivars that are null-mutants for glycitein accumulation. strikingly, gmf6h unlike the other members of the gmf6h family lacks the characteristic n-terminal hydrophobic transmembrane domain that anchors cyp enzymes into the er. gmf6h gene expression, as analyzed by qpcr, did not correlate with high isoflavonoid accumulation. the discrepancy between qpcr and differential expression analysis results for gmf6h gmi2′h and gmf3′5′h could indicate that the de results were not accurate to the same level as qpcr analysis, as these genes fell in the ‘baseline’ p-value score. alternatively, the gene expression and function of a different isoform of the f6h or i2′h gene families could correlate with high isoflavonoid content, as assigning reads to highly similar genetic regions is problematic. enzymatic candidates for a role in specialized metabolism table  and table  are a compilation of manually selected candidates from the high stringency de analysis that might be associated with the isoflavonoid content of the four cultivars in question. in addition to gmdfr that is significantly downregulated, there are several cyps, 2-odds and fad/nadp-dependent oxidoreductases in the list that could be top candidates for involvement in specialized metabolism. overrepresentation and pathway enrichment statistics show that genes involved in secondary metabolism, particularly those encoding hydrolase and oxidoreductase activity are significantly enriched in the ‘low isoflavonoid’ gene list. membrane integral proteins are also overrepresented within this list, which could impact the formation of isoflavonoid or flavonoid multi-enzyme complexes or ‘metabolons’. 2-odds are a very versatile family of dioxygenases that are superior to cyps in the range of reactions they catalyze, from hydroxylations to ring fragmentation. they are ubiquitously distributed throughout nature, and have very important roles in core and specialized metabolic pathways, including histone demethylation, and the biosynthesis of gibberellic acid, flavonoids, benzylisoquinoline alkaloids, and glucosinalates, to name a few. it is important to note that flavanone 3β-hydroxylase is a 2-odd, converting a 2s-flavanone to dihydroflavonols by a c- hydroxylation. this is an enzyme that competes with isoflavonoid biosynthesis, and downregulation of a corresponding gene could underlie increased isoflavonoid content. there are several reactions in the flavonoid and isoflavonoid pathways that have yet to be elucidated and table  and table  could provide a starting point towards gene discovery. among them there are four putative 2-odd genes, two that are upregulated: glyma.10g.p and glyma.17g.p, and two downregulated members: glyma.08g.p and glyma.14g.p. other metabolic enzymes represented in the list of high stringency de genes include a downregulated homolog of cyp71b, and several fad/nad-binding oxidoreductases. the latter group has been less present within the current literature on phenylpropanoid pathway; however an nadph/fad-dependent enzyme was shown to be responsible in the additional hydroxylation of positions and/or of the flavonol a-ring this provides a precedence for fad/nad-linked oxidoreductases being involved in the content and composition of phenylpropanoids, and should be further researched. stress response in high isoflavonoid cultivars isoflavonoid biosynthesis is inducible by several abiotic and biotic stimuli. therefore, it would be interesting to further investigate the ‘high isoflavonoid’ genes that are annotated as being responsive to stress. subsets of these genes are annotated as being involved in disease resistance such as, dirigent proteins, which dictate the stereochemistry of other proteins, and several nucleotide binding disease resistance proteins. the latter included nb-lrr proteins that are associated with r gene function. several of these genes were annotated as being root hair specific. as root hairs are often the site for pathogen entry, the convergence of tissue-specialization and function could be indicative of a role for these genes in inhibiting pathogen infection or colonization. as mentioned earlier, conrad and ac colombe are also favored for their increased resistance to p. sojae and other soybean pathogens. therefore, the differential expression of genes involved in stress response, localized to the root hairs could be suggestive of a mechanism underlying the improved resistance to pathogen infection. transcriptional regulation: transcription factors and chromatin regulators subcellular localization of the ‘high isoflavonoid’ genes to the nucleus highlights the potential role of the cognate proteins in transcriptional regulation. the list of nuclear-localized genes includes transcriptional regulators, including members of the bhlh superfamily of proteins and myb transcription factors. physical interaction and regulatory synergy has been shown between these family of proteins in regulating gene expression. combinatorial plant gene regulation might be a factor in the coordinated redirection of flux in specialized metabolism, leading to isoflavonoid content variation in soybean cultivars. a gras family transcription factor was significantly upregulated in high isoflavonoid cultivars. based on homology to arabidopsis transcription factors it is putatively denoted for stress transduction. as isoflavonoid production is closely related to stress response, this transcription factor could be an important candidate for further study. another function for nuclear-localized proteins could be the regulation of chromatin structure, which thereby effects broad transcriptomic changes. two genes upregulated and annotated for such function were identified: a set7/ family, histone h3k4-specific methyltransferase protein and a histone deacetylase family protein. the former, a histone-modifier protein, cannot be immediately associated with activation or repression of transcriptional regulation, as methylation of h3k would have to be regarded in the larger landscape and context of methylations. gmhda on the other hand, could be associated with the downregulation of certain genes, explaining the transcriptomic changes between cultivars. deacetylation conferred by elements such as gmhda can increase dna-histone affinity and condense the overall structure, rendering it inaccessible to transcription machinery. hypoacyetylated chromatin is silent, or has reduced gene expression. therefore, gmhda might be a histone-modifying element responsible for some of the overall changes in the root transcriptome that coincide between the high isoflavonoid cultivars, conrad and ac colombe. vacuolar sequestration: a putative regulatory mechanism for isoflavonoid accumulation isoflavone aglycones produced by an enzyme complex tethered to the er, are conjugated by transferases, allowing their appropriate channeling into the vacuole. from this storage point, isoflavonoids can be sequestered outside of the metabolic flux and released upon demand for stress response or as signaling molecules for symbiosis with nitrogen-fixing rhizobia. interestingly, there were up- and down-regulated genes annotated for transferase activity; up- and down-regulated genes with putative transporter activity; the suite of transferases and transporters being simultaneously regulated could denote a preference towards the conjugation and storage of a certain class of metabolites. however, there is still ambiguity over the specificity of vacuolar transporters involved in the channeling of specialized metabolites, and broad substrate acceptance within classes of compounds has been reported. interestingly, one of the candidate upregulated genes identified in the de analysis of soybean cultivars, was an uncharacterized member of the mate family, and hereafter referred to as gmmate mate transporters are important families of proteins in the glycosylation and the subsequent transport of isoflavones. gmmate was consistently upregulated in high isoflavonoid cultivars, suggesting a possible role in the sequestration of isoflavonoids. the predicted localization of this protein was to the secretory pathway, based on the presence of signal peptide. this could imply localization to the er, golgi bodies or the vacuole, with the last option being highly probable, given its functional annotation as a mate transporter. the deduced amino acid sequence of this gene had an amino acid sequence similarity of % with the tt homolog in arabidopsis, based on clustal w, the latter being a known proanthocyanidin vacuolar transporter. mate transporters in medicago truncatula, mate and mate are involved in transport of phenylpropanoids. mate is a tonoplast epicatechin 3′-o-glucoside transporter. mate is a flavonoid transporter involved in vacuolar sequestration of anthocyanins and other flavonoids in flowers and leaves. mate transporter prefers malonylated flavonoid glucosides and is co-expressed with three genes encoding malonyltransferases. the discovery of an isoflavonoid-specific mate transporter, in turn, could be used as a tool to investigate co-expressed transferases, depicting a more complete image of the conjugation and transport of isoflavonoids. furthermore, up-regulation of gmmate could indicate a mechanism for the vacuolar transport and sequestration of isoflavonoids and/or precursors, in competition with parallel pathways for metabolic flux.  the coordinated regulation of genes encoding enzymes, transporters, transcription factors and other molecular elements could lead to the increase of root isoflavonoids in conrad and ac colombe. these transcriptomic elements can help overcome bottlenecks in isoflavonoid production for soybean cultivars, legumes, and heterologous systems. genes annotated for disease and stress response could be important in conferring soybean resistance to p. sojae infection in conrad and ac colombe. functional genomics should be employed to characterize the role of such transcriptomic elements in isoflavonoid production and disease resistance. additional files additional file 1: table s sequence of oligonucleotides used for qpcr of soybean genes in high and low root isoflavonoid cultivars to confirm differential expression analysis from rnaseq. additional file 2: figure s estimated intra-cultivar dispersion for all genes. dispersion values plotted as a function of average normalized expression for each gene model. biological variation displayed for conrad ac colombe ac glengarry and pagoda in three biological replicates. additional file 3: figure s m-a plot of log fold change as a function of average expression of gene. significantly differentially expressed genes are represented in red. plots displayed for pair-wise comparisons ac glengarry with ac colombe ac glengarry with conrad pagoda with ac colombe pagoda with conrad. additional file 4: figure s gene frequency as a function of benjamini-hochberg adjusted p-values. histogram of p-values for all genes assessed for significant differential expression between two cultivars ac glengarry with ac colombe ac glengarry with conrad pagoda with ac colombe pagoda with conrad. additional file 5: figure s genes ranked by total read count against –log of the p-value. significance scores divided into three levels:,, and. average gene counts divided into percentiles. four plots representing significance of differential expression between cultivars ac glengarry with ac colombe ac glengarry with conrad pagoda with ac colombe pagoda with conrad. additional file 6: table s list of genes upregulated in high as compared with low root isoflavonoid content cultivars. differentially expressed genes in the four comparisons between high and low cultivars were analyzed for overlap, generating a set of candidates upregulated consistently in high root isoflavonoid cultivars. additional file 7: table s list of genes downregulated in high as compared with low root isoflavonoid content cultivars. differentially expressed genes in the four comparisons between high and low cultivars were analyzed for overlap, generating a set of candidates downregulated consistently in high root isoflavonoid cultivars. additional file 8: table s list of genes upregulated in high as compared with low root isoflavonoid content cultivars. highly differentially expressed genes in the four comparisons between high and low cultivars were analyzed for overlap, generating a set of candidates upregulated consistently in high root isoflavonoid cultivars. additional file 9: table s list of genes downregulated in high as compared with low root isoflavonoid content cultivars. list of genes downregulated in high as compared with low root isoflavonoid content cultivars. highly differentially expressed genes in the four comparisons between high and low cultivars were analyzed for overlap, generating a set of candidates downregulated consistently in high root isoflavonoid cultivars. these genes were annotated using the soybean database and have been compiled below. additional file 10: figure s heatmaps of genes upregulated in high root isoflavonoid cultivars. read counts for gene models were normalized across cultivars. four heatmaps representing genes upregulated in high isoflavonoid cultivars as compared with low root isoflavonoid cultivars were generated: ac glengarry with ac colombe ac glengarry with conrad pagoda with ac colombe pagoda with conrad. additional file 11: table s genes upregulated in ‘high isoflavonoid’ cultivars were analyzed for overrepresentation of panther go-slim classifications for: biological process, molecular function and cellular component. additional file 12: table s genes downregulated in ‘high isoflavonoid’ cultivars were analyzed for overrepresentation of panther go-slim classifications for: biological process, molecular function and cellular component. additional file 13: table s genes upregulated in ‘high isoflavonoid’ cultivars were analyzed for pathway enrichment. additional file 14: table s genes downregulated in ‘high isoflavonoid’ cultivars were analyzed for pathway enrichment. additional file 15: table s genes annotated with disease or stress response and upregulated in ‘high isoflavonoid’ cultivars. additional file 16: table s genes annotated for transcriptional regulation and upregulated in ‘high isoflavonoid’ cultivars. abbreviations chichalcone isomerase chrchalcone reductase chschalcone synthase cypcytochrome p450 dedifferentially expressed gmdfrsoybean dihydroflavonol 4-reductase gmf3′hsoybean flavonoid 3′-hydroxylse gogene ontology 
a transcriptomic approach using soybean cultivars with high and low root isoflavonoid content was used to find elements that underlie this variation. coordinated regulation of genes involved in flavonoid metabolism could redirect flux into the isoflavonoid branch of the phenylpropanoid pathway, by reducing competition for the flavanone substrate. these candidate genes could help identify mechanisms to overcome the endogenous bottleneck to isoflavonoid production, facilitate biosynthesis in heterologous systems, and enhance crop resistance against pathogenic infections. two genes, encoding the flavonoid-metabolizing enzymes, flavonoid 3′-hydroxylase and dihydroflavonol 4-reductase, had lower expression levels in high isoflavonoid cultivars. they play a role in signaling for symbiosis with nitrogen-fixing bacteria and inhibiting pathogen infection. in addition, genes annotated with stress and disease response were upregulated in high isoflavonoid cultivars.
 isoflavonoids are a class of specialized metabolites found predominantly in legumes. they play a role in signaling for symbiosis with nitrogen-fixing bacteria and inhibiting pathogen infection.  a transcriptomic approach using soybean cultivars with high and low root isoflavonoid content was used to find elements that underlie this variation. two genes, encoding the flavonoid-metabolizing enzymes, flavonoid 3′-hydroxylase and dihydroflavonol 4-reductase, had lower expression levels in high isoflavonoid cultivars. these enzymes compete with isoflavonoid biosynthetic enzymes for the important branch-point substrate naringenin and its derivatives. differentially expressed genes, between the two sets of cultivars, encode transcription factors, transporters and enzymatic families of interest, such as oxidoreductases, hydrolases and transferases. in addition, genes annotated with stress and disease response were upregulated in high isoflavonoid cultivars.  coordinated regulation of genes involved in flavonoid metabolism could redirect flux into the isoflavonoid branch of the phenylpropanoid pathway, by reducing competition for the flavanone substrate. these candidate genes could help identify mechanisms to overcome the endogenous bottleneck to isoflavonoid production, facilitate biosynthesis in heterologous systems, and enhance crop resistance against pathogenic infections. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords soybeansecondary metabolismphenylpropanoidisoflavonoidroottranscriptomerna sequencinghttp://dx.doi.org/ /501100000040agriculture and agri-food canada j-000151dhaubhadel sangeeta issue-copyright-statement© the author 2017 