Mining transcriptomic data to study the origins and evolution of a plant allopolyploid complex

Allopolyploidy combines two progenitor genomes in the same nucleus. It is a common speciation process, especially in plants. Deciphering the origins of polyploid species is a complex problem due to, among other things, extinct progenitors, multiple origins, gene flow between different polyploid populations, and loss of parental contributions through gene or chromosome loss. Among the perennial species of Glycine, the plant genus that includes the cultivated soybean (G. max), are eight allopolyploid species, three of which are studied here. Previous crossing studies and molecular systematic results from two nuclear gene sequences led to hypotheses of origin for these species from among extant diploid species. We use several phylogenetic and population genomics approaches to clarify the origins of the genomes of three of these allopolyploid species using single nucleotide polymorphism data and a guided transcriptome assembly. The results support the hypothesis that all three polyploid species are fixed hybrids combining the genomes of the two putative parents hypothesized on the basis of previous work. Based on mapping to the soybean reference genome, there appear to be no large regions for which one homoeologous contribution is missing. Phylogenetic analyses of 27 selected transcripts using a coalescent approach also are consistent with multiple origins for these allopolyploid species, and suggest that origins occurred within the last several hundred thousand years.


INTRODUCTION
Polyploidy (whole genome duplication, WGD) is a key process in plant evolution. All seed plants are fundamentally polyploid, with a second WGD event shared by all flowering plants (Jiao et al., 2011), and additional events found in many lineages (see http://genomevolution.org/wiki/index.php/Plant paleopolyploidy) (Soltis et al., 2009). It has been estimated that 15% of all flowering plant speciation events involve polyploidy (Wood et al., 2009). Systematists generally recognize autopolyploidy and allopolyploidy as distinct types of polyloidy events, based on the level of divergence of the diploid genomes that formed the polyploid. The terms are best thought of as describing elements of a continuum that ranges from the doubling of a single genome (autopolyploidy), to the incorporation of differentiated genomes in a single nucleus by hybridization of different species (allopolyploidy). From a genetic perspective, allopolyploids are characterized by diploid-like meiotic behavior and limited interaction between the two homoeologous genomes. The duplicated chromosomes of an autopolyploid (and, to a lesser extent, a newly formed allopolyploid; Ramsey & Schemske, 2002) initially can associate randomly, leading to polysomic segregation, but it is generally assumed that this is a transient state; diploidization leads to the eventual presence of homoeologous genomes. It is difficult, if not impossible to determine from the genomes of older polyploids (paleopolyploids, mesopolyploids) how differentiated their progenitor genomes were in large part due to the frequent absence of extant diploid progenitors for comparative purposes.
The initial "fixed hybrid" condition of an allopolyploid erodes over time as homoeologous loci are lost (Lynch & Conery, 2000;Maere et al., 2005); this process of "fractionation" is thought to occur preferentially from one subgenome, but the precise mechanisms remain unknown (Schnable & Freeling, 2011;Freeling et al., 2012). In addition to the loss of genes, the process of concerted evolution can result in the replacement of a gene from one genome by its homoeologue, notably through gene conversion (e.g., Wang et al., 2007). The earliest stages of polyploid evolution may contribute disproportionately to gene loss and genomic rearrangement through genomic shock (McClintock, 1984). For example, some individuals of the ca. 100 year-old allopolyploid, Tragopogon miscellus, have lost entire chromosomes of one parent (Chester et al., 2012). Diversity in polyploids can be due to mutational divergence from parental diploids, but also due to multiple origins produced by different polyploidization events between different genotypes of the same diploid species (Symonds, Soltis & Soltis, 2010). Questions concerning how polyploids originate (e.g., single vs. multiple origins), how they partition their variation (e.g., as a single lineage united by gene flow vs. as separate lineages formed from different genotypes of the same progenitor species), and how much of the initial parental contributions they retain are among the major questions in polyploid evolutionary research .
High-throughput sequencing produces massive amounts of genome-wide data, and thus has great potential for systematic and evolutionary studies in general (Gilad, Pritchard & Thornton, 2009). The ready availability of genomic and transcriptomic data has opened new opportunities for studying the evolution of polyploids (Bombarely et al., 2012;Grover, Salmon & Wendel, 2012;Ilut et al., 2012;Dufresne et al., 2014) at the scale of whole genomes. However, it is not trivial to extract relevant information from short read sequencing data, particularly for allopolyploids, where the interest is often in deconvoluting the complex genome into its two homoeologous subgenomes (Grover, Salmon & Wendel, 2012;Ilut et al., 2012). Moreover, the field of systematics has what has been called a new paradigm for studying species relationships, involving genealogical approaches (Edwards, 2009). Genealogical methods have lately begun to be applied to both autopolyploids (Arnold, Bomblies & Wakeley, 2012;Hollister et al., 2012) and allopolyploids (e.g., Slotte et al., 2011;Jones, Sagitov & Oxelman, 2013;Slotte et al., 2013). The confluence of these two developments promises to accelerate the study of polyploid evolution.
The genus Glycine includes the cultivated soybean (G. max) and its wild progenitor (G. soja), both annual species native to northeastern Asia, as well as approximately 30 perennial species native to Australia classified as subgenus Glycine (Ratnaparkhe, Singh & Doyle, 2011). Like many plant species, Glycine has a complex history of polyploidy: in addition to events shared with all angiosperms (Jiao et al., 2011) and eudicots (Jiao et al., 2012), the soybean genome retains evidence from a WGD around 50 million years ago (MYA) shared with a large subset of legumes (Blanc & Wolfe, 2004;Schlueter et al., 2004;Cannon et al., 2010), and particularly from a more recent polyploidy event that increased the chromosome number of the ancestor of all extant Glycine species from 2n = 20 to 2n = 40 (Shoemaker, Schlueter & Doyle, 2006;Doyle & Egan, 2010;Schmutz et al., 2010;Doyle, 2012). This Glycine-specific WGD occurred between the estimated time of homoeologous gene divergence in the soybean genome (10-13 MYA; e.g., Egan & Doyle, 2010;Schmutz et al., 2010), and around 5 MYA, when the annual and perennial species diverged from an already-polyploid common ancestor (Doyle & Egan, 2010).
In addition to these older events, eight perennial Glycine species are allopolyploids with 2n = 78 or 80, hypothesized to have arisen by hybridization involving various combinations of eight extant diploid species, several of them multiple times and involving both progenitors as chloroplast genome donors . Various lines of evidence culminated in these hypotheses of reticulate relationships, which are shown in Fig. 1 for the six species that are part of the G. tomentella sub-complex . Chromosome number polymorphism (2n = 38,40,78,80) was observed in what was initially considered a single taxon, Glycine tomentella (Newell & Hymowitz, 1978). Patterns of sterility and partial chromosome pairing in artificial crosses among G. tomentella plants were consistent with the presence of shared homoeologous diploid genomes among polyploids (Grant, Brown & Grace, 1984;Doyle et al., 1986;Singh, Kollipara & Hymowitz, 1998). Isozyme studies of diploid and allopolyploid G. tomentella led to the characterization of numerous "races" designated either "D" for diploid, or "T" for tetraploid (Doyle & Brown, 1985;Singh, Kollipara & Hymowitz, 1998). Morphological complexity, presumably due to the reticulate nature of the complex, has slowed nomenclatural recognition of what are clearly species in the biological sense. More recently, molecular phylogenetic studies assumed a dominant role in refining hypotheses of relationships (Hsing et al., 2001;Singh, Kollipara & Hymowitz, 1998;Brown et al., 2002;Doyle et al., 2002;, and corroborated earlier hypotheses concerning the origins of polyploids from among the diploid (2n = 38,40) "genome groups" that were also initially defined by artificial hybridization studies and later by molecular studies (see Ratnaparkhe, Singh & Doyle, 2011). However, these DNA-level studies were based on only two molecular markers: the internal transcribed spacers of the 18S-5.8S-26S nuclear ribosomal gene cistron (nrDNA ITS) and the low copy nuclear gene, histone H3D. Relationships of chloroplast genomes are broadly consistent with these results (Hsing et al., 2001), but are complicated by incongruence with nuclear markers, likely due to a combination of incomplete lineage sorting and introgression . Thus, a genome-wide perspective on the  (Ratnaparkhe, Singh & Doyle, 2011) are given for diploids. Species used in this study (G. tomentella D1, G. tomentella D3, G. syndetika D4, G. canescens, G. clandestina, G. dolichocarpa T2, G. tomentella T1 and G. tomentella T5) are shown in green. origin and evolution of the G. tomentella complex, including estimates of dates of origin, has been lacking.
A better understanding of the origin and evolution of the Glycine allopolyploid complex will complement its exploitation in studying the impact of allopolyploidy on a range of morphological and physiological characters (Coate & Doyle, 2010;Coate et al., 2012;Ilut et al., 2012;Coate & Doyle, 2013;Hegarty et al., 2013). Here we apply phylogenetic and coalescent methods to a transcriptomic dataset from three of these allopolyploid species and their diploid progenitors that was originally generated to study the effects of polyploidy on their ability to cope with stress from excess light (Coate & Doyle, 2013).
For each light treatment, all tissue was collected in a single morning and immediately frozen in liquid nitrogen. Total RNA was isolated from pooled leaf tissue using the Plant RNeasy Kit with on-column DNase treatment (Qiagen, Valencia, CA, USA). Single-end RNA-Seq libraries were constructed following the Illumina mRNA-seq Sample Preparation Kit protocol (Illumina, San Diego, CA, USA), with the following modifications: (1) two rounds of polyA selection were performed using the Dynabeads mRNA DIRECT Kit (Life Technologies, Carlsbad, CA, USA); (2) RNA was fragmented for 2 min at 70 • C using the RNA fragmentation reagents kit (Life Technologies); and (3) Illumina PE adapters were replaced with custom-made adapters containing 3nt barcodes in order to facilitate multiplexing of samples (see Coate & Doyle, 2013 for adapters and Table S1 for the barcode sequences). Sequencing was performed on either the GAIIx or HiSeq 2000 platform (Illumina), generating 88 nt or 100 nt reads, respectively. Equimolar amounts of three (GAIIx) or four (HiSeq 2000) barcoded libraries were combined and sequenced per channel.

Read processing and single nucleotide polymorphism (SNP) calling
Reads were processed with Fastq-mcf (Aronesty, 2013) to trim low quality extremes (min. quality 30) and remove short reads (min. read length 50 bp). They were aligned to the soybean genome (version 1.0, downloaded from www.phytozome.net/soybean) using Bowtie2 (Langmead & Salzberg, 2012) with the default parameters. Mapping files from the same accession were merged. Reads without preferential mapping (same score for two or more mapping hits) and with a mapping score below 20 were removed. SNP calling was performed using Samtools (Li et al., 2009). SNPs supported with read coverage below 5 were removed. VCF files were combined and formatted to Structure and Hapmap formats using the Perl script MultiVcfTool (https://github.com/aubombarely/GenoToolBox/blob/ master/SeqTools/MultiVcfTool).

Homoeologue read identification and transcript-guided assembly
For homoeologous SNP identification, a consensus diploid transcriptome was rebuilt for each of the species groups (A, with G. clandestina and G. canescens accessions; D1, with G. tomentella D1 accessions; D3, with G. tomentella D3 accessions; and D4, with G. syndetika accessions) using Samtools (Li et al., 2009) and Gffread from the Cufflinks software package (Trapnell et al., 2010). A progenitor reference set was created for each of the polyploid species joining the diploid transcriptome sets (T1 = D1 + D3, T2 = D3 + D4 and T5 = A + D1). Reads from the polyploid species were mapped with these references using Bowtie2. Sam mapping files were processed to identify reads according the preferential mapping with each of the progenitors using the Perl script, SeparateHomeolog2Sam (https://github.com/aubombarely/GenoToolBox/blob/master/ SeqTools/SeparateHomeolog2Sam). Reads with mapping score AS and XS = 0 (No SNPs) were kept and used to rebuild the polyploid transcriptomes using Samtools (Li et al., 2009) and Gffread (from the Cufflinks package Trapnell et al., 2012). Once the reads were separated according its preferential mapping, they were mapped back to the soybean genome. SNPs were called as described above.

Population structure analysis
The programs Structure (Pritchard, Stephens & Donnelly, 2000) and fineStructure (Lawson et al., 2012) were used to analyze population structure of the two SNP datasets, with and without polyploid SNPs separated by homoeologue, described above. For Structure, each of the datasets was divided into three subsets of 20,000 SNPs selected with a random function incorporated in the MultiVcfTool. 5 replicates were run for each of the subsets with a burn-in of 10,000 and a number of MCMC repetitions of 10,000, from K = 1 to K = 15 using the default parameters (λ = 1, assuming uniform distribution of allele frequencies, Pritchard, Stephens & Donnelly, 2000). Admixture was selected. The optimal number of clusters was identified based on the rate of change in the log probability of data between successive K values (Evanno, Regnaut & Goudet, 2005). Results at K = 6 were verified with a re-analysis using a burn-in of 100,000 generations. Results were visualized using R (barplot function).
The two SNP datasets were divided into 20 different subsets each mapping to one soybean reference chromosome for FineStructure analysis. Analyses were performed following the instructions from the fineStructure web for the unlinked model (http:// www.maths.bris.ac.uk/∼madjl/finestructure/data example.html). Results were presented as a heatmap of distances between each of the accessions. A principal component analysis (PCA) was performed over the same distance matrix using fineStructure software. The PCA figure was created using R.

Reconstruction of phylogenies using concatenated SNPs
SNPs from the dataset in which SNPs from allopolyploids were partitioned into their two homoeologues ("homoeologue data set") and were concatenated to create a supermatrix with 36 operational taxonomic units (OTUs). The two homoeologous gene copies from each allopolyploid were treated as individual OTUs; for example the D1 and D3 homoeologues of T1 individuals were treated as D1T1 and D3T1, respectively. G. max, accession William82 was used as outgroup. The alignment files were produced changing the SNPs Hapmap format to fasta using a Perl script. The resulting matrix was used in two analyses. First, maximum likelihood (ML) was used, implemented in PhyML (Guindon & Gascuel, 2003) with GTR as the substitution model; 100 bootstrap replicates were conducted. Second, in order to visualize reticulations in the dataset, a network method, NeighborNet, was implemented in the SplitsTree package (Huson & Bryant, 2006) with the default parameters. Trees were visualized and drawn using FigTree (Rambaut, 2012).

Gene-based analyses
A subset of transcripts was selected for phylogenetic and network analyses based on the following criteria: No more than 10% of Ns for the guided assembly consensus sequence in any of the accessions after the homoeologue read identification; alignments with at least 1000 bp; and genes with their corresponding G. max homologue identified as an existing pair retained from the most recent (ca. 5-10 million years; Doyle & Egan, 2010) Glycine WGD event, as compiled by Du et al. (2012). Sequence alignments were based on the transcriptome-guided assembly. Sequence for each of the genes was collected with a Perl script (FastaSeqExtract, GenoToolBox script package), concatenated and changed to the required sequence alignment format using a BioPerl script (bp sreformat.pl). The 95 alignments selected were used in an exploratory phylogenetic analysis using the Bayesian MCMC method, BEAST (Drummond et al., 2012) (HKY substitution model, 10,000,000 MCMC). Alignments that produced trees in which G. max was not sister to perennial Glycine species in the consensus tree were removed. Generally the removed alignments showed tree topologies with two large clades with long branches, indicating the possibility of inclusion of paralogous genes from the older whole genome duplication (ca. 50 MY, common to the Leguminosae; reviewed in Doyle, 2012) instead the orthologue.
27 genes selected after this filtering were analyzed using three different methods: (1) Phylogenies were reconstructed using ML using PhyML (Guindon & Gascuel, 2003) with 1,000 bootstraps. jModelTest2 was used to choose the best substitution model (Darriba et al., 2012). According to the Bayesian Information Criterion (BIC) HKY was the preferred model (40% of the genes), followed by K80 (26% of the genes; Table S2).
(2) Networks were constructed using NeighborNet in SplitsTree4 with the default parameters (Huson & Bryant, 2006). (3) Bayesian analysis was performed using BEAST v2.0 (Drummond et al., 2012). The two homoeologous gene copies from each allopolyploid were treated as individual OTUs as in the concatenated analysis, and G. max, accession William82 was again used as outgroup. Based on the jModelTest2 results, HKY was used as the substitution model. The MCMC chain was set to 100,000,000 MCMC generations, taking samples every 1000 generations. Divergence ages were estimated by scaling the tree root (divergence between G. max and perennials) to 5 Myr (Egan & Doyle, 2010). All trees were drawn using FigTree (Rambaut, 2012).

Species tree reconstruction
Species tree reconstruction under the coalescent was performed using the 27 selected genes in *BEAST (Drummond et al., 2012). The 24 accessions, including two homoeologues for each allopolyploid accession, were grouped in 11 operational taxonomic units (OTUs) for this analysis: G. max was used as outgroup. Based on jModelTest2 results, HKY was used as substitution model. The MCMC chain was set to 100,000,000 MCMC generations, taking samples every 1000 generations. Divergence dates were estimated as described above. All the trees were drawn using FigTree (Rambaut, 2012).

Phylogenomics dataset generation
Between 7 and 60 million reads from leaf transcriptomes of 24 accessions representing 8 Glycine perennial species were mapped to the Glycine max genome (v1.0) (Schmutz et al., 2010). Reads mapped to 22,500-25,000 genes (∼40% of soybean gene models; Table 1); this represents between 4.5 and 11.6% of the genome. 200,000-965,000 single nucleotide polymorphisms (SNPs) were identified relative to G. max; 6.3-12.6% of SNP positions were polymorphic in diploid species (G. clandestina, G. canescens, G. tomentella D1 (referred as D1 hereafter), G. tomentella D3 (referred as D3) and G. syndetika (referred as D4)), and 18.4-28.8% in polyploid species (G. tomentella T1 (referred as T1), G. tomentella T5 (referred as T5) and G. dolichocarpa T2 (referred as T2); Table 2). The interpretation of these positions as standard heterozygosity is complicated by the recent (5-10 MYA: Doyle & Egan, 2010) WGD in the ancestral Glycine genome. In a gene for which soybean has lost one of the homoeologous copies from this event, but the perennial species for which it is serving as reference has retained both copies, polymorphic SNPs may be due to reads from two different homoeologous loci in the perennial, rather than two alleles at a single locus. Low levels of conventional heterozygosity are expected in Glycine species, because of their strongly selfing reproductive biology, with much reproduction occurring through cleistogamous (closed, selfing) flowers. The much higher percentage of polymorphic positions in polyploid individuals (T1, T2, T5) likely is also due to the mapping of reads from two homoeologous copies to a single target, in this case due to recent polyploidy: for example, mapping reads from tetraploid (2n = 80) T2 to a single locus in the diploid (2n = 40) G. max reference genome will result in reads from both its D3 and D4 homoeologous subgenomes mapping to the same target, increasing the chance of observing a polymorphism at a given site. Separating reads from T1, T2, and T5 polyploid individuals was possible where the read has at least one SNP that could be related to one homoeologous genome contributor (e.g., D3 and D4 differed by a SNP) and this difference was retained in the D3 and D4 homoeologous genomes of T2; diploid-distinguishing polymorphism (DDP; see Ilut et al., 2012). Between 11.4 and 20.8% of reads were assigned to one of the progenitors (Table 3).
Between 124,984 and 399,884 SNPs were produced for each accession. The filtering of the missing data produced 237,243 and 75,958 polymorphic positions for all the accessions before and after the homoeologous read assignment, respectively. SNPs per chromosome ranged from 7,455 (chromosome 14) to 16,494 (chromosome 8) and from 2,288 (chromosome 14) to 5,300 (chromosome 8) before and after the homoeologous read assignment, respectively. SNPs per species group ranged from 21,830 (D1 species) to 26,438 (A species, G. canescens and G. clandestina) ( Table 4).
Transcriptome-guided assemblies produced between ∼1,800 and ∼6,600 full-length sequences (as mapped to the G. max gene models) for each diploid accession. For polyploid subtranscriptomes this number was much lower because only reads that mapped preferentially to one of the diploid consensus species and reads that mapped equally but with no polymorphism (perfect match) were used during the transcriptome-guided assembly. Any read that mapped equally to two or more positions with one or more polymorphisms was discarded because it was impossible to assign it to any of the diploid progenitors, reducing the mapping coverage of the reference gene models. Between ∼350 and ∼1,350 full length sequences were assembled for the T1, T2, and T5 polyploid homoeologous subtranscriptomes of which between 4 and 19% were duplicated genes from the 5 to 10 MYA WGD event in the common ancestor of Glycine species (Schmutz et al., 2010). For phylogenetic analysis, full length sequences are not needed so a phylogenetic analysis dataset was created with 27 genes (see Material and Methods for the criteria used to generate this dataset; Table 5).

Genome-wide distribution of homoeologous SNPs
For each allopolyploid accession, the ca. 120,000-400,000 SNPs (Table 3) that could be identified to the homoeologous subgenome were mapped to the soybean reference genome (Schmutz et al., 2010). This produced a map that is analogous to chromosome painting (genomic in situ hybridization, GISH) experiments using the reads from which the SNPs were derived, which we term "electronic chromosome painting" (e-chromosome painting). Similar patterns were seen for all accessions, with high densities of SNPs at the ends of each soybean chromosome and far lower densities in pericentromeric regions (Fig. 2). This pattern is expected using reads from transcriptome data, because of the sparse distribution of genes in pericentromeric regions of the soybean genome (Schmutz et al., 2010). Notably, in all allopolyploid accessions, SNPs from both homoeologues were distributed across the entire genome, and no regions were identified in which SNPs from only one homoeologue were mapped ( Fig. 2; Figs. S1-S10).

Population structure analyses
Structure (Pritchard, Stephens & Donnelly, 2000) was first run using all available SNPs, without separating SNPs from polyploids into homoeologous groups. Structure was run from K = 1-15; K = 6 was identified as one of the optimal preferred values of K using the delta K method of Evanno, Regnaut & Goudet (2005; Fig. S11). Five of these six groups corresponded to diploid taxa: D1, D3, D4, G. canescens, and G. clandestina (Fig. 3A). The sixth group was represented only as a minor component in D4 accession 2073. Diploid accessions showed little or no evidence of admixture, with the exception of D4 accession 2073 (Fig. 3). In contrast, all polyploid accessions were admixed, each with  (Joly et al., 2004). A second Structure analysis was conducted with each polyploid accession treated as two separate OTUs, using the homoeologue dataset ( Table 2). As with the previous analysis, the analysis was run for K = 1-15. The Evanno method (Evanno, Regnaut & Goudet, 2005) identified K = 6 and 9 as the preferred values (Fig. S11). In the case of K = 9 the group representation shows the same structure than the K = 6 (Fig. S12). Results for diploids were similar to those obtained in the previous analysis (Fig. 3B). Subgenomes from natural allopolyploids and the synthetic T5 allopolyploid (A58) were shown to belong exclusively to diploid groups, with little or no evidence of admixture, indicating that the SNP filtering into homoeologous contributions was successful. Complementary to the second Structure analysis, the data were analyzed using ChromoPainter and FineStructure (Lawson et al., 2012). ChromoPainter produces a co-ancestry matrix (as a measure of the ancestry sharing between individuals) based on the haplotype information provided by shared chunks (regions) of biallelic markers between individuals (Lawson et al., 2012). The two SNP datasets were filtered by selecting only the biallelic markers, producing a subset with 220,952 and 71,610 SNPs (before and after homoeologous read assignment, respectively) distributed along all 20 soybean chromosomes. Regions identified by ChromoPainter for each accession ranged from 516 (D4 2321) to 567 (G. clandestina 1253) and from 202 (D4 1300 and 2321) to 221 (D4 2073) (before and after homoeologous read assignment respectively). Principal component analysis (PCA) and population relationship analysis using a Bayesian approach were performed over the co-ancestry matrix using FineStructure (Lawson et al., 2012). PCA before homoeologous read assignment (Fig. 4A) shows seven well-differentiated groups, one per species with the exception that G. canescens and G. clandestina clustered together. Diploid species formed the vertices of a trapezoid. A-genome species (G. canescens, G. clandestina and D4) formed a more dispersed group than either D1 or D3. Each polyploid species fell between its putative diploid progenitors, consistent with each being an admixture (fixed hybrid). After the homoeologous read assignment (Fig. 4B), each of the polyploid subgenomes clustered with its diploid progenitors, producing three clear clusters: D1, D3, and A-genome (comprising G. canescens, G. clandestina and D4, as expected). Heatmaps were used to visualize the population relationships produced by FineStructure, complementing the information shown by the PCA figures. The heatmap before homoeologous read assignment (Fig. 4C), showed four intense regions (red, magenta and blue colors) corresponding to the four species groups of the PCA (Fig. 4A). Each polyploid showed the expected similarity to its progenitors; similarly, as expected the two G. clandestina accessions were more similar to one another than either was to G. canescens. Also, T5 A58, the artificial polyploid produced from a cross between G. canescens 1232 and D1 1316, showed the expected relationships with these accessions. Other T5 polyploids also showed a stronger signal from D1 1316 than from other D1 accessions. T2 accessions did not show any stronger signal with any particular D3 accession than with others, but they did with the D4 accessions 1300 and 2321, relative to 2073. T1 accessions 1288 and 1763 also showed a stronger signal with particular D1 and D3 accessions, whereas T1 accession 1361 showed a weaker signal with the D1 and D3 accessions included here. After the homoeologous read assignment (Fig. 4D), some of these signals were intensified, such as the relationship between T5 and D1 subgenomes and particular D1 accessions, but other relationships that were suggested when all SNPs were considered were not observed (for example there is not a stronger signal of D1 1316 with the T5 accessions). These differences may be due to the methodology used for the homoeologous read assignment.

Phylogeny and network analysis of concatenated SNPs
Phylogenetic and network analyses were conducted using the homoeologue dataset, with SNPs concatenated to create a single supermatrix. The maximum likelihood (ML) tree, rooted with G. max, identified four subclades comprising two major clades: (1) the A-genome, with subclades of D4 vs. G. clandestina and G. canescens; and (2) the D-genome (D3) and E-genome (D1) (Fig. 5A). Each of the subclades showed a different pattern with respect to diploid and tetraploid subgenome relationships. In the canescens/clandestina clade, the A-subgenome of the synthetic allopolyploid (A58) was sister to the accession from which it was created (G. canescens 1232), as expected, though with deeper coalescence than expected from an artificial hybrid; the two natural T5 allopolyploids were sister to G. clandestina 1126, as expected from other data (e.g., Doyle et al., 2002). In the D4 clade, diploid accession 2073 was sister to all remaining accessions, a unique placement consistent with its apparently admixed nature (Fig. 3A). The polyploid subgenomes formed a paraphyletic group, with the two diploid accessions sister to the D4 subgenome of one T2 accession (1134). A similar pattern was seen in the D3 subclade, where T2 accessions formed a paraphyletic group, and all four diploid accessions formed a clade sister to T2 accession 1134. Also embedded within the T2 accessions was a clade consisting solely of T1 accessions. T1 accessions also formed a monophyletic group within the D1 clade, where natural T5 accessions and D1 accessions also formed monophyletic groups. Surprisingly, there was not a sister relationship between the D1-subgenome of synthetic allopolyploid A58 and the D1 accession from which it was formed (1316). Similar topologies were produced by neighbor-joining analysis (data not shown).
NeighborNet was used to analyze the full homoeologue dataset to identify minority patterns of relationships in the data. When rooted with G. max, the topology (Fig. 5B) was very similar to the ML tree (Fig. 5A), even having such features as the sister relationship of D4 2073 to other D4 accessions, and the monophyly of T1 homoeologues in both the D1 and D3 clades. There was clear evidence of character support for alternative relationships, but those relationships were minor in comparison with the major phylogenetic signal.

Gene-based phylogenetic and network analyses
Gene trees were constructed for the 27 genes (described in the Material and Methods) using several different phylogenetic and network methods. Similar topologies for trees from individual genes were obtained with BEAST and PhyML. All 27 trees showed the split between the A-genome clade and the D1/D3 clade seen in the ML tree reconstructed from concatenated SNPs (Fig. 5A). However, many individual gene trees showed unexpected groupings of one or more accessions, particularly within the A-genome clade, where several trees grouped accessions from G. canescens with G. syndetika-D4 instead of with G. clandestina (for example ML and BEAST trees for the gene Glyma04g39670, Figs. S17 and S45). Relationships within the major subclades varied among the 27 gene trees. For example, nine of the 27 trees showed separate clades for G. canescens (plus the A58 sequence) and G. clandestina (e.g., Figs. 6A and 6C), but in only three of them did diploid species form monophyletic groups (Figs. S13-S67). Overall, there were far more departures from expectations in the A-genome clade than in the D1/D3 clade.
There were numerous cases where alleles from diploid accessions formed monophyletic groups (e.g., 12 of 27 BEAST topologies had alleles from all four D3 accessions in a clade, often with high posterior probability). At some loci, alleles from one or more polyploids formed monophyletic clades; for example, at Glyma06g18640 (Fig. S50), all taxa, including both homoeologous subgenomes of each polyploid, formed separate clades, with the exception of G. clandestina. However, this was unusual, and paraphyletic groupings of alleles were common, particularly in polyploids. For example, at 26 of 27 loci, T2-D3 alleles were not monophyletic, at least some having closer relationships to D3 or T1-D3 alleles, and in gene Glyma01g35620, T5-D1 1969 was most closely related to D1 1156 whereas T5-D1 1487 was most closely related to D1 1157 and D1 1316 (Fig. S41). On the assumption that alleles in tetraploids all originated from diploid progenitor species, such paraphyletic relationships suggest the input of alleles from different genotypes of diploid progenitors, due either to multiple origins or, alternatively, to continued gene flow from diploids after polyploid formation, perhaps involving unreduced gametes.
The BEAST trees, calibrated with the 5 MYA divergence of G. max and the perennial subgenus (Innes et al., 2008), allowed dates of allele divergence to be estimated. Among comparisons of interest are the minimum divergences between alleles from a tetraploid and alleles from its diploid progenitor (e.g., T2-D3 vs. D3) or alleles from the same progenitor in a second tetraploid (e.g., T2-D3 vs. T1-D3); the latter represent "diploid" alleles as well, under the assumption that there has been no gene flow between the two tetraploids, something that is reasonable for G. tomentella tetraploids (e.g., Doyle et al., 1986). Minimum distances between polyploid and diploid alleles (over)estimate the time of entry of that allele into the polyploid, which is typically assumed to be an origin of the polyploid (Doyle & Egan, 2010). Minimum dates (Table S3) were 0.31 MY for T1 (measured at the D1 locus), 0.29 MY for T5 (measured at the D1 locus), and 0.38 MY for T2 (measured at the D3 locus). Error bars on these estimates, however, were substantial.
NeighborNet (implemented in SplitsTree 4; Huson & Bryant, 2006) was used to construct networks for each of the 27 genes. Several networks showed patterns consistent with intragenic recombination; the Pairwise Homoplasy Index (PHI) of Bruen, Philippe & Bryant (2006), also implemented in SplitsTree, was significant for 11 of the 27 genes (data not shown). The dominant patterns in NeighborNet topologies were similar to the overall pattern shown in phylogenetic analyses of the 27 genes, and thus to results for the full homoeologous SNP dataset. As with other methods, NeighborNet networks suggested multiple inputs of alleles from diploid progenitors into polyploids (e.g., gene Glyma02g11580, Fig. 6C).

Species tree reconstruction under the coalescent
Species trees were reconstructed using the coalescent approach implemented in *BEAST (Heled & Drummond, 2010), which used information contained in the individual gene trees from the 27 genes described above. The overall *BEAST tree (Fig. 7A) topology was similar to that of trees from concatenated SNPs. By definition, each of the allopolyploid homoeologous genomes was a single OTU despite the possibility of independent origins; each of these was grouped with its putative progenitor species. Within the D1 genome clade, the T1 and T5 polyploids were sisters to one another; similarly, T1 and T2 were sisters in the D3 clade. The DensiTree output (Fig. S40) indicated considerable uncertainty only within the D3 clade, where both other possible topologies (T2 sister to D3, T1 sister to D3) appeared in a substantial number of trees. As expected, divergence dates of polyploids from their diploid progenitors estimated by *BEAST were higher than minimum estimates from the 27 individual loci, all being greater than 300,000 years (Fig. 7A).
To enhance the utility of this model group, it is important to move to a genome-wide understanding of their biology. As noted above, origins of the Glycine allopolyploids were hypothesized initially from crossing data and more recently from gene phylogenies, but inferences have been made from only two nuclear genes. Both of these markers supported the hypotheses of fixed hybridity of Glycine allopolyploid species. However, it is not known to what extent the entire genomes of these plants retain contributions from both parental diploid species in the face of potential loss due to initial genomic shock (McClintock, 1984), or other processes such as "genome downsizing" (Leitch & Bennett, 2004), fractionation (Schnable & Freeling, 2011;Freeling et al., 2012), or concerted evolution (e.g., Wang et al., 2007).

Glycine allopolyploids are fixed hybrids throughout their genomes
Analyses using all SNPs identified from the full dataset showed that all three of these allopolyploids are indeed fixed hybrids, combining diploid genomes as depicted in Fig. 1.
Structure results indicated an essentially equal contribution from both parental diploids in all three cases (Fig. 3A); PCA analysis also was consistent with this hypothesis, placing each polyploid approximately midway between its putative progenitors, as expected for an F1 hybrid (Fig. 4A).
In order to determine whether or not the polyploids have contributions from their parents across their entire genomes, reads were partitioned by homoeologous genome and mapped to the soybean reference genome (Schmutz et al., 2010). As portrayed by e-chromosome painting (Fig. 2), it is clear that no individual sampled from any of the three allopolyploid species has any major regions represented by only one homoeologue. Coverage is sparse in pericentromeric and centromeric regions, as expected due to the low density of genes in these regions of the soybean genome (Schmutz et al., 2010). The degree of shared synteny between soybean and these perennial Glycine species is as yet unknown, but regardless of the order of chromosomal segments, it is clear that there has not been significant loss of homoeologous genes. We mapped reads to over 22,000 of the approximately 46,000 genes of the soybean genome (Schmutz et al., 2010). These numbers include both homoeologous copies from the 5-10 MYA polyploidy event that shaped the modern "diploid" (2n = 38,40) Glycine genome. We were able to deconvolute between 4 and 19% of these 22,000 genes into their homoeologous contributions in each of the three recent allopolyploids (e.g., T1D1 and T1D3). Using genomic in situ hybridization (GISH), Chester et al. (2012) showed examples of allopolyploid T. miscellus plants that had all four chromosomes or chromosome segments of one diploid parent (4:0), but also examples of plants with 3:1 ratios of homoeologous chromosomes or chromosomal segments. Our e-chromosome painting method cannot distinguish the 3:1 condition from an equal contribution from both parents segments, so it is possible that such plants exist in our sample.
Structure analysis using the partitioned homoeologous SNPs corroborated results with the full, unpartitioned dataset, in placing each polyploid homoeologous genome with its putative progenitor (Fig. 3B). The FineStructure PCA supported three major groupings, each of which included diploids and the expected polyploid homoeologous subgenomes derived from them (Fig. 4B). The grouping of D4 accessions and two A-genome species (G. canescens and G. clandestina), along with polyploid genomes derived from them, into a single cluster is not surprising, because G. syndetika (D4) is also a member of the A-genome (Ratnaparkhe, Singh & Doyle, 2011). As noted above, genome groups were originally defined on the basis of reproductive compatibility in artificial crosses (Ratnaparkhe, Singh & Doyle, 2011), and indeed G. syndetika (D4) 2073 shows evidence of admixture with G. canescens and G. clandestina (Fig. 3). In contrast, D1 and D3, though both classified as "G. tomentella", belong to two different genome groups (E and D, respectively;Ratnaparkhe, Singh & Doyle, 2011). This greater genetic similarity of the three A-genome species is not reflected in relative divergence dates; for example, the *BEAST analysis dates the divergence between G. syndetika and the two other A-genome species at slightly earlier than the divergence between D1 and D3 (Fig. 7A). Thus, reproductive barriers likely arose earlier in the D1/D3 lineage than within the A-genome. Allopolyploid evolution in Glycine fits "Darlington's Rule" (Darlington, 1937)-that allopolyploids should form between species that are reproductively isolated, often due to chromosomal differences, whereas reproductively compatible diploids instead tend to form homoploid hybrids. No allopolyploids are known to have formed among A-genome species, and only one of the eight known Glycine allopolyploids involves hybridization within a genome group (tetraploid G. tabacina is the product of the most divergent species cross possible within the B-genome; . D1 and D3, which as noted belong to different genome groups, have different chromosome numbers (2n = 38 vs. 40, respectively), which may contribute to their inability to form fertile diploid hybrids. D1 has also formed allopolyploids with D5A, another 2n = 40 "G. tomentella"; however, reproductive incompatibility also occurs between 2n = 40 G. tomentella taxa (Doyle et al., 1986), and other allopolyploids in the complex combine genomes of two 2n = 40 taxa (Fig. 1).

Gene histories, allele divergence times, and sources of genetic diversity in polyploids
Gene trees from the 27 loci were selected that met criteria designed to provide orthologues. These genes are highly transcribed with sufficient characters for phylogeny reconstruction, and inferences of polyploid origins mostly conformed to expectations based on previous work using the low copy nuclear locus, histone H3D (Brown et al., 2002;Doyle et al., 2002;González-Orozco et al., 2012), the nrDNA ITS (Singh, Kim & Hymowitz, 2001;, and chloroplast noncoding sequences (Hsing et al., 2001). The use of BEAST and *BEAST (Heled & Drummond, 2010) allowed us to estimate divergence times of alleles and species for the first time for some of these taxa. Dating polyploid origins is complicated by numerous factors (Doyle & Egan, 2010). For one thing, if the polyploid has arisen recurrently, then there is no single date that marks "the" origin. Given that polyploids are often invasive (e.g., Pandit, Pocock & Kunin, 2011), and the Glycine tomentella allopolyploids appear to be recently formed based on sharing identical histone H3D and nrDNA ITS alleles with their putative progenitors (Doyle et al., 2002;, we have speculated that they could have originated as a response to ecological disturbance due to human colonization of Australia, around 40,000 years ago (Hudjashov et al., 2007;Pugach et al., 2013). The relevant date for testing this anthropogenic disturbance hypothesis would be the oldest origin of each polyploid. However, because it is unlikely that a polyploid allele and any of a set of diploid progenitor alleles will coalesce at exactly the time of polyploid origin, distances for any given polyploid event will be overestimates of the actual time of origin. Further complicating matters, the error bars on our BEAST divergence estimates were large relative to the estimates themselves. Nevertheless, because even the minimum estimates of allele divergences between diploids and tetraploids are around 0.3 MY, it appears likely that these G. tomentella allopolyploids are hundreds of thousands rather than tens of thousands of years old. *BEAST estimates should be averages of all origins of a polyploid taxon, and these, too are several hundred thousand years for each allopolyploid. Thus, it appears likely that these polyploid species were present in Australia before humans arrived there. The fact that these three species, and possibly other allopolyploid members of the complex, may have evolved at roughly the same time is intriguing. In the ca. 5 MY since the perennial members of Glycine diverged from the annual lineage (Egan & Doyle, 2010), there is no evidence of polyploidy until these species were formed, apparently well within the last 1 MY. Perhaps the onset of severe aridity in Australia around 3 MYA, heralding the change to the present extreme wet-dry glacial cycles (Crisp, Cook & Steane, 2004) could have provided ecological opportunities for polyploids. It will be interesting to refine our estimates through increased sampling of these three triads, and to obtain estimates for the other five allopolyploid species.