Domestication of Soybean

Soybean [Glycine max (L.) Merr.] is a remarkable plant, producing both high quality oil and protein and is one of the primary row crops in the United States. Although soybean is relatively new to western agriculture, it has been under cultivation for > 3000 years [1, 2]. The transition from wild Glycine soja to cultivated Glycine max was the result of ancient plant breeders/farmers selecting for a large number of domestication-specific traits (photoperiod insensitivity, lack of shattering, lack of lodging, seed size increases, seed set increases, etc.). Dramatic changes in seed oil/protein content and fatty acid composition have apparently also been selected for during domestication, either directly or indirectly [3, 4].

Genetics of soybean seed coloration

The compounds responsible for soybean seed coat and hilum color in soybean are derivatives of phenylpropanoid pathway [9–11] (Figure 1). The wild type condition of black seed coats is primarily due to two anthocyanidin glycosides (anthocyanins): cyanidin-3-monoglucoside and delphinidin-3-monoglucoside [10, 11]. In lines which feature brown seed coats, only cyanidin is apparently present at maturity [10]. Aside from the cosmetic and aesthetic aspect of coloration, anthocyanins are thought to have diverse human health promoting capabilities [12].

The action of UDP-glucose:flavonoid 3-O-glucosyltransferase enzymes is a critical step in anthocyanin accumulation

Two anthocyanin glycosides form the predominant colored compounds in black seed coats: cyanidin-3-monoglucoside and delphinidin-3-monoglucoside [10]. These are formed through the action of UDP-glucose:flavonoid 3-O-glucosyltransferase (UF3GT) enzymes, which specifically transfer a glucose moiety from UTP to the 3' position of cyanidin and delphinidin (recently reviewed in [13], Figure 1). This glycosylation is thought to increase the stability and solubility of the cyanidin molecule [14]. In lines with brown seed coats (r), cyanidin accumulates, though high levels of proanthocyanidins are also present [10]. Recently, two highly similar co-expressed UF3GT genes (Glyma07g30180 and Glyma08g07130) were determined to be expressed in seed coats of black seeded soybean lines, and these genes have been demonstrated to specifically transfer a glucose moiety to the cyanidin molecule at the 3'-hydroxyl group, resulting in the formation of cyanidin-3-glucoside [15].

The Inhibitorlocus

Seed coat color is primarily under control of the Inhibitor locus, which has at least four classically defined genetic alleles [8], listed here from the most dominant to the least: I (largely colorless seeds) > i ⁱ (color restricted to hilum) > i ^k ("saddle;" color in hilum and spreading slightly beyond the hilum) > i (seeds completely black). Inhibitor acts in a dominant, gain-of-function manner with maternal-effect inheritance, and results in seed coats appearing pale yellow due to the absence of anthocyanins [10]. Both the dominant Inhibitor allele (I) and the i ⁱ alleles have been shown to be due to naturally occurring, gene-silencing effects derived from linked but independent Chalcone Synthase (CHS) gene clusters (chromosome 8, LG A2) that generate siRNA which target CHS gene transcripts specifically within the seed coat for degradation [16–22].

The genetics of soybean hilum coloration

Lines which have the dominant I allele can still exhibit some traces of color within the hilum, with the specific hilum coloration due to the allelic status at three other genetic loci: Tawny, R, and W1 [6, 8] (Table 1). Hilum tissue is not maternally-derived, in contrast to the seed coat [23]. In lines with the recessive (i) allele, seed coat color is brown, imperfect black, buff or black, dependent on the allelic status of the Tawny, R and W1 loci (Table 1).

The Tawny locus has two pleiotropic effects: homozygosity for the gray (t) allele results in gray pubescence at maturity and, in lines carrying the combination of the i ⁱ allele of the Inhibitor locus, a functional R gene, and purple flowers (W1), seed which feature "imperfect black" hila (Table 1). Alternatively, gray pubescent (t) lines carrying the i ⁱ allele of the Inhibitor locus, a functional or nonfunctional R, and white flowers (w1) produce seed which feature buff hila [8] (Table 1). The phenotypic effects of the recessive allele of Tawny have been discerned to be due to loss of function mutations affecting a flavonoid 3' hydroxylase gene (Glyma06g21920) [24]. At the chemical level, this is the result of a reduction in the accumulation of anthocyanins within the hilum, and the presence of pelargonidin (Figure 1), which does not accumulate in lines carrying the wild type version of the Tawny locus [10].

The recessive allele of the Rlocus is responsible for brown hilum/seed coats

Another locus, classically termed R, also interacts epistatically with the Tawny and Inhibitor loci (as well as the W1 locus) to control hilum and seed coat colors [8] (Table 1). Lines with a functional Tawny gene and homozygous for the recessive allele of the R locus possess either brown seed coats or brown hilum, dependent on the allelic status of the Inhibitor locus (i or i ⁱ respectively). Although the genetics behind this trait were well resolved shortly after the rediscovery of Mendel's work in the 1920s [6], the molecular genetic basis has not been ascertained. Despite this, the ease of phenotyping has resulted in the inclusion of this locus in the development of genetic maps for soybean [25–27].

Epistasis for genes involved in soybean coloration

Epistatic and pleiotropic interactions are the norm for genes involved in soybean coloration (Table 1). For example, loss of function mutations affecting a flavonoid 3'5'-hydroxylase gene (w1, F3'5'H, Figure 1) have been demonstrated to result in two phenotypes: white flowers and loss of purple pigment in hypocotyls [28]. The allelic status of the W1 locus, when combined with the recessive gray allele of the Tawny locus, determines if seed coats or hila are colored "imperfect black" or "buff" (Table 1) [8].

Approaches to identify the r locus, which results in brown hilum/seed coats

Loss of function mutations affecting a gene involved in the terminal end of the anthocyanin biosynthetic pathway have been suggested as the cause of the recessive brown seed coat/hilum phenotype (Figure 1). Possible candidates have included UF3GT, Anthocyanidin Synthase (ANS) and/or Dihydroxyflavone Reductase (DFR) genes. However, no correlation has been found between the genomic locations of any UF3GT, DFR or ANS gene and the location of the R gene [29]. Alternately, a transcription factor or other regulatory element could be responsible for the brown hilum/seed coat phenomenon. The objective of this work was to identify the specific gene and causative basis behind the phenomenon of brown hilum/seed coat coloration, historically defined as the R locus, in soybean.

RIL population development

The generation of the F₆ RIL mapping population, derived from a cross between Jake X PI 283327, was previously described [30]. Jake (PI 643912) has tawny pubescence, purple flowers, and shiny yellow seed with black hila (i ⁱ T R W1)[31]. The brown hila line, PI 283327 has tawny pubescence, purple flowers, and yellow seed with brown hila (i ⁱ , T, r, W1) (USDA GRIN germplasm collection, accessed 06/22/2011 (http://www.ars-grin.gov/npgs/). The reference cultivar Williams 82, for which the genome sequence was determined [32], has tawny pubescence, white flowers and yellow seed with black hila (i ⁱ , T, R, w1) [33].

Bulk segregant analysis of selected RIL lines

A total of 100 F₆ RIL lines were selected from a Jake X PI 283327 cross in which segregation for hilum color had occurred (50 possessed black hila, 50 had brown hila) and seed from each were pooled to form two bulks. Only RILs that were definitively black or brown were used in the bulks, with ambiguous or mixed RILs not included. The seeds (1 per RIL) were ground utilizing a coffee grinder to generate a fine powder. The grinder was cleaned thoroughly between grindings. DNA was isolated using a DNeasy Plant Maxi Kit (Qiagen, Inc., Valencia, CA) according to manufacturer's recommendations. Bulk DNA was concentrated using standard ethanol precipitation procedure to yield a final concentration of 3.52 micrograms mL^-1 (black bulk) or 2.40 micrograms mL^-1 (brown bulk). Bulk DNA was used with Universal Soybean Linkage Panel (USLP) as previously described [34].

Simple Sequence Repeat (SSR) markers

All SSR primer pairs from within the newly delimited R locus region, drawn from a bioinformatically defined list, were also examined for potential utility in fine-mapping [35]. Fine mapping PCR was performed in 20 microliter reactions as previously described [36] and PCR products were separated on 2% agarose gels. Genotypic classes were assigned by visual comparison to PCR reactions using DNA from parental lines. Only those SSR primer pairs which showed obvious, easily scored size polymorphism between the two parents (PI 283327 and Jake) were used in subsequent analysis. SSR primers pairs which displayed polymorphism within the newly defined R region, and which could theoretically be used to select for this trait, are listed in Additional File 1.

DNA isolation, PCR and sequencing of candidate genes from pureline seed

DNA was isolated using a DNeasy plant mini kit (Qiagen), and 5-50 ng of DNA were used in PCR with Ex taq (Takara) with gene specific primers (Additional File 1) under the following conditions: 95°C for 5 minutes, followed by 40 cycles of 95°C for 30 seconds, 59°C for 30 seconds, and at 72°C for 1 minute per 1 kbp of predicted product size. Following PCR, products were examined on a 1% agarose gel by electrophoresis and sent for sequencing at the University of Missouri DNA core facility. Sequence traces were downloaded, imported into Contig Express model of the VectorNTI Advance 11 software (Invitrogen, Carlsbad, CA, USA), assembled and manually evaluated for polymorphisms. Putative polymorphisms were verified by a second, independent PCR and sequencing reaction.

Selection of diverse lines from the germplasm repository

136 lines were selected for sequencing of the putative R gene, drawn either from a previously established list of diverse germplasm [37] or were individually selected from the USDA GRIN germplasm collection (http://www.ars-grin.gov/npgs/) to ensure a broad geographic distribution with a range of hilum and seed coat colors. Certain color classes were only minimally investigated, due to epistatic interactions which precluded novel information (e.g. yellow seed coat with buff hila, see Table 1). A full listing of the 136 lines examined for the allelic status of the R gene/Glyma09g36990 is listed in Additional File 2. For a subset of ten lines, all three exons were examined by sequencing (including the 5' UTR, 3'UTR, the 1^st intron and the majority of the 2^nd intron, although portions of the 2^nd intron are highly repetitive AT-rich and recalcitrant to PCR and sequencing). These lines were: PI 84970 (Hokkaido Black, black seed coats), PI 518671 (Williams 82, yellow seed coats, black hila), PI 643912 (Jake, yellow seed coats, black hila), PI 548461 (Improved Pelican, yellow seed coats, brown hila), PI 548389 (Minsoy, yellow seed coats, brown hila), PI 438477 (Fiskeby 840-7-3, yellow seed coats, brown hila), PI 180501 (Strain #18, yellow seed coats, brown hila), PI 283327 (Pingtung Pearl, yellow seed coats, brown hila), PI 240664 (Bilomia No. 3, yellow seed coats, brown hila), PI 567115 B (MARIF 2782, black seed coats). Because all mutations identified were found to affect the 1^st or 2^nd exons, we elected to only sequence the first and second exons (as well as 5' UTR, the 1^st intron, and a portion of the 2^nd intron) in the remaining 126 lines.

qRT-PCR

Expression analysis on seed coat, cotyledon or leaf total RNA (DNAse-treated using Turbo DNase (Ambion, Austin, TX, USA)) was performed as described [38] with minor modifications. The RT-PCR mix was supplemented with 0.2X Titanium Taq polymerase (BD Biosciences, Palo Alto, CA) to improve primer efficiency. Following the reverse transcriptase reaction, amplification was 95°C for 15 min, then 35 cycles of 95°C for 20 seconds, 60°C for 20 seconds, and 72°C for 20 seconds. Primers used in this work are listed in Additional File 1. The reference gene used to normalize data was CONS6 [39] and raw Ct values were first applied to efficiency curves developed for each primer set utilizing Williams 82 genomic DNA, then normalized to the expression of the reference gene and expressed as a percent of CONS6.

Numerous researchers have reported reliable data from qRT-PCR utilizing RNA from mature yellow seed coat tissue. However RT-PCR using RNA derived from brown seed coat tissue was challenging, likely owing to the known effect of interference due to proanthocyanins [10]. The use of a simple PCR Inhibitor removal column (Zymo, Irvine, CA, USA) remedied this difficulty, resulting in acceptable qRT-PCR data derived from mRNA isolated from maturing brown seed coat tissue.

We also investigated CHS7/8 using a primer pair previously described [18]; however the results were highly variable in both cotyledon and seed coat tissues with no significant expression level differences detected between the brown and black seed coat samples (data not shown).

Bulk Segregant Analysis

In order to identify the gene responsible for the r locus effect (brown hilum or seed coat color), we initially utilized the bulk segregant analysis (BSA) [40] method on RILs from a population derived from the cross of soybean cultivar Jake with the PI 283327 which had segregated for the R gene alleles with the USLP array [34]. Although this technique confirmed the previously identified location of the R locus [25, 26], the extremely broad window identified (data not shown, ~4.2 Mbp, based on the Williams 82 sequence) failed to further delimit the boundaries of the R locus.

Identification of lines featuring recombination events within the delimited Rregion

Three primer pairs from the novel SSR set (BARCSOYSSR 09_1475, 09_1501 and 09_1566 were examined for all 100 RIL lines. For the majority, the hilum color phenotype was correlated with the expected parental polymorphic band. We also observed seven individual RILs which possessed recombination events within the region identified on chromosome Gm09/LG K (Figure 2A). We examined these seven RILs using all novel polymorphic SSRs markers within this region, and compared the marker genotype to the RIL phenotype (Table 2 Figure 2A). Our methodology allowed us to fine-map the location of the R gene to a predicted region of less than 200 kbp with only 100 RIL lines. This region in Williams 82 contains 23 predicted open reading frames, with another 3 genes annotated as pseudogenes (Figure 2B).

Identification of four R2R3 MYB genes as candidates for the Rlocus

BLAST searches using the 26 candidate genes were performed against NCBI (http://www.ncbi.nlm.nih.gov/) and TAIR (http://www.arabidopsis.org) databases to search for candidate genes. BLAST searches revealed four tandem genes which featured homology to the R2R3 MYB transcription factor gene family: Glyma09g36970, Glyma09g36980, 09g36990 and Glyma09g37010. R2R3 MYB genes have been shown to control flux through the phenylpropanoid pathway, and mutants in multiple species are associated with changes in fruit, flower and/or seed color (recently reviewed in [41]). These four tandem R2R3 MYB genes are highly similar (~80-90% nucleotide identity, excluding presumed intronic sequence) and may have arisen due to a tandem gene amplification event(s). Strikingly, none of these genes appears to have been identified in recent seed focused studies using RNAseq methods [42, 43].

Expression analysis of R2R3 gene candidates

Because soybean hilum tissue is extremely small and difficult to accurately dissect from seeds in non-pigmented stages, we utilized a large seeded soybean line with brown seed coats (PI 567115 B) and a large seeded line with black seed coats (PI 84970) to examine mRNA expression. In order to assess whether a subset of these four tandem genes were pseudogenes and/or expressed in seed coat tissue (Glyma09g36970 is annotated as a pseudogene in the current whole genome shotgun sequence build), we utilized qRT-PCR. Only one of these candidate R2R3 MYB genes, Glyma09g36990, was expressed in any of the tissues examined (leaf, seed cotyledons, and seed coats). Gene transcripts from Glyma09g36990 were present in the seed coats of both a brown seeded and a black seeded cultivar. However, this gene was not expressed in either cotyledon tissue (Figure 3A) or in leaves (data not shown). It is not clear if the other three R2R3 MYB genes in the cluster are expressed in other tissues. Nor is the role these genes play in soybean physiology known, if any.

Curiously, the Williams 82 Glyma09g36990 gene model was predicted to possess four exons, in contrast to the canonical 3 exons identified for authentic R2R3 MYB transcription factor genes [44, 45]. To characterize the authentic expressed sequence, RT-PCR was used to analyze full length cDNA for comparison to the reference Williams 82 gene model. The authentic gene is slightly larger than that the predicted Glyma09g36690 gene model and possesses three exons (Additional File 3), in concordance with that reported for other R2R3 MYB genes [44, 45].

Analysis of Glyma09g36990 for potential causative polymorphisms

PCR and Sanger sequencing of exons (and partial intronic sequence) was used to evaluate the Glyma09g36990 gene for polymorphisms in a selection of lines: Jake (black hilum), PI 283327 (brown hilum), Williams 82 (black hilum), PI 84970 (black seed coats) and PI 567115 B (brown seed coats). We discovered a single-base deletion within exon 2 in PI 283327 and PI 567115 B that results in a frameshift mutation (C377-, relative to the start codon) (Figure 4, details in Additional File 3). The open reading frame for Glyma09g36990 was allelic between Williams 82, Jake and PI 84970.

We then elected to examine a broad geographic distribution of lines (136 in total, Additional File 2) from the available soybean germplasm corresponding to all of the known seed coat and hilum color classes. From this pool, we identified three additional presumed loss of function mutations: G343-, resulting in frameshift; G95C TGG > TCG (W32S) missense in conserved residue; AGgt > AGtt (g404t) disrupts conserved mRNA splice recognition site (Figure 4, further details in Additional File 3).

In all cases where we observed an intact open reading frame, we noted the phenotype of imperfect black hilum (i ⁱ R t W1), buff hilum (i ⁱ R t w1), black hilum (i ⁱ R T) or black seed coat (i R T), dependent on the allelic status of the Inhibitor and Tawny loci (Additional File 2). Any of these four loss of function alleles resulted in either brown hilum (i ⁱ r T), brown seed coat (i r T) or buff hila (i ⁱ r t). In all cases, we observed a perfect association between the presence of one of the four loss of function alleles and brown hilum or brown seed coats, barring the presence of the epistatic dominant I allele or gray pubescence, both of which can mask the effect of the r allele, resulting in yellow or buff hila (Additional File 2). These epistatic interactions (and masking in the case of Inhibitor) are due to the placement of the step affected by the R2R3 MYB gene at the terminal end of the anthocyanin biosynthesis pathway (Figure 1). Any one of the loss of function mutations affecting the R gene are necessary and sufficient for brown seed coat and/or hilum coloration. However, the phenotypic effect can be masked or modulated by the presence of certain alleles of the Inhibitor and Tawny loci (Table 1 Additional File 2).

Time-course of mRNA expression for Glyma09g36990 and two phenylpropanoid biosynthetic enzymes

If the candidate R gene is controlling expression of a gene which forms a rate limited step in anthocyanin production, we hypothesized that a correlation would exist between 1) R gene expression levels, 2) the appearance of color compounds, and 3) the expression of ANS and/or UF3GT genes in developing seed coats. We examined a time course of seed coat and seed cotyledons by qRT-PCR (Figure 3) for expression of three genes: the R gene candidate, ANS, and UF3GT. Seed coats from the large seeded line with brown seed coats (PI 567115 B) and one with black seed coats (PI 84970) were investigated for quantitation of steady state transcripts. We selected four time-points corresponding to the development of pigmentation during seed growth and maturation for PI 84970 (black seed coats) and PI 567115 B (brown seed coats) (Additional File 4).

Although there are apparently two UF3GT genes expressed in seed coats in soybean (Glyma07g30180 and Glyma08g07130), only one of these genes (Glyma08g07130) is not expressed in cotyledon tissue [15]. We elected to focus on this gene for qRT-PCR, as we noted a virtual absence of ANS or R gene expression in cotyledons (Figure 3A and 3B).

We observed a near-perfect coefficient of correlation (R² = 0.96) between the level of expression (relative to an internal control CONS6) of the putative R gene and a UF3GT gene (Glyma08g07130) (Figure 3A and 3C). In contrast, we observed a weak correlation between expression of the R gene and ANS gene expression (R² = 0.66) in the black seed coat line (Figure 3A and 3B). In the brown seeded line PI 567115 B, no significant correlation was found between R gene expression levels and either ANS or UF3GT expression levels (Figure 3A-C). During early and mid-development stages R gene expression is similar in both black and brown seed coat lines, though R gene expression declined during the last stages of development of the brown seeded line, in contrast to the high expression noted for the black seed coat lines (Figure 3A). In striking contrast to the increase in expression of ANS and UF3GT during seed coat maturation of the black seed coat line, only negligible ANS and UF3GT expression was observed in the brown seed coat line as seeds approached maturity (Figure 3B and 3C).

These findings confirmed our hypothesis that loss of function mutations within Glyma09g36690, an R2R3 MYB gene, are correlated with reduced expression of a UF3GT gene and ANS genes and with the brown hilum/seed coat phenotype. It remains to future work to determine the specific DNA sequence targeted by the soybean R2R3 MYB R gene product and its specific interactions in complexes with basic-helix-loop-helix (bHLH) transcription factors and WD40 proteins. It is unclear if the R gene product acts to promote transcriptional activation of both ANS and UF3GT genes, or if activation of ANS gene expression is due to an indirect effect.