A Phylogenomic Assessment of Ancient Polyploidy and Genome Evolution across the Poales.

Comparisons of flowering plant genomes reveal multiple rounds of ancient polyploidy characterized by large intragenomic syntenic blocks. Three such whole-genome duplication (WGD) events, designated as rho (ρ), sigma (σ), and tau (τ), have been identified in the genomes of cereal grasses. Precise dating of these WGD events is necessary to investigate how they have influenced diversification rates, evolutionary innovations, and genomic characteristics such as the GC profile of protein-coding sequences. The timing of these events has remained uncertain due to the paucity of monocot genome sequence data outside the grass family (Poaceae). Phylogenomic analysis of protein-coding genes from sequenced genomes and transcriptome assemblies from 35 species, including representatives of all families within the Poales, has resolved the timing of rho and sigma relative to speciation events and placed tau prior to divergence of Asparagales and the commelinids but after divergence with eudicots. Examination of gene family phylogenies indicates that rho occurred just prior to the diversification of Poaceae and sigma occurred before early diversification of Poales lineages but after the Poales-commelinid split. Additional lineage-specific WGD events were identified on the basis of the transcriptome data. Gene families exhibiting high GC content are underrepresented among those with duplicate genes that persisted following these genome duplications. However, genome duplications had little overall influence on lineage-specific changes in the GC content of coding genes. Improved resolution of the timing of WGD events in monocot history provides evidence for the influence of polyploidization on functional evolution and species diversification.

Introduction

Paleopolyploidy or ancient whole-genome duplication (WGD) events have occurred across the eukaryotic tree of life and have been hypothesized to have had major impacts on life history innovations and organismal diversification (Levin 1983; Taylor and Raes 2004; Soltis et al. 2009; De Smet and Van de Peer 2012). WGDs have been especially widespread throughout flowering plant history (Cui et al. 2006; Vanneste et al. 2014), and ancient polyploidy has been associated with the origin of seed plants (spermatophyta) and angiosperms (Jiao et al. 2011). In fact, multiple rounds of WGD have been inferred for many plant lineages (Tang et al. 2010; Jiao et al. 2011, 2012; McKain et al. 2012; Yang et al. 2015) including those leading to the model plant species Arabidopsis thaliana (Blanc and Wolfe 2004), the tree species Populus trichocarpa (Tuskan et al. 2006), legumes (Fabaceae; [Cannon et al. 2015]), and grasses (Poaceae; [Paterson et al. 2004; Tang et al. 2010]). Though many recent polyploid lineages have relatively lower net diversification rates than diploid lineages, genome duplications have undoubtedly contributed to molecular, ecological, and phylogenetic diversification during angiosperm evolution (Mayrose et al. 2011; Soltis et al. 2014). However, a mechanistic understanding of how polyploidy has contributed to evolutionary innovations and diversification in angiosperm history requires precise phylogenetic placement of the nature and timing of these ancient WGDs (Soltis et al. 2009).

Polyploid can occur through either somatic genome doubling within meristematic tissue, zygotes, or young embryos, or with the fusion of unreduced gametes (Ramsey and Schemske 1998). Both autopolyploidy (multiple copies of the same genome) and allopolyploidy (multiple copies of different genomes, usually of different species) have prevalence in angiosperm history (Barker et al. 2016). The formation of stable, sexually reproducing polyploid populations is thought to involve two steps: 1) mating between individuals with reduced and unreduced gametes, sometimes producing individuals with odd polyploidy levels, the so-called triploid bridge (Mallet 2007), and 2) subsequent fusion of unreduced gametes producing a stable tetraploid, or higher even-level polyploid, such as hexaploid wheat (Levin 2002). Morphological, physiological, and ecological diversification among polyploid populations can follow as a consequence of differential gene loss or subfunctionalization and neofunctionalization of retained duplicates (i.e., homeologs).

Polyploidy is hypothesized to have played a major role in the evolution of key innovations through the proliferation of novel genes and gene interactions (Ohno 1970; Levin 1983; Van de Peer et al. 2009). Polyploidy was recently shown to be a driver of innovation and novelty in secondary metabolites in the mustard family (Edger et al. 2015), demonstrating a connection between genome duplication, adaptation, and diversification. Additionally, polyploidy has been implicated as spurring epigenetic changes (Yoo et al. 2014), immediate and sustained variation in gene expression (Yoo et al. 2013), and reorganization of gene networks (De Smet and Van de Peer 2012). The process of duplicate gene loss (fractionation) can exhibit biases toward some functional classes (Maere et al. 2005; Doyle et al. 2008; Freeling 2009; Tang et al. 2012), and potentially contribute to phenotypic variation due to differential duplicate retention. In some grass family clades (such as the tribe Andropogoneae), polyploidy has been shown to be a recurring evolutionary process with multiple independent polyploid events occurring over a relatively short period of time (Estep et al. 2014).

Early analyses of the rice genome revealed a WGD event designated as rho (Paterson et al. 2004; Yu et al. 2005) inferred to have predated the divergence of the BOP and PACMAD clades within the grass phylogeny (Grass Phylogeny Working Group 2001; Blanc and Wolfe 2004; Paterson et al. 2004; Schlueter et al. 2004; Yu et al. 2005; Wu et al. 2008; Grass Phylogeny Working Group II 2012; Soreng et al. 2015). The rho event has been implicated as contributing to the success of Poaceae through the role of duplicated MADS-box genes in the development of the spikelet (Preston and Kellogg 2006; Preston et al. 2009) and the role of duplicated genes involved in the development of starch-rich seeds (Wu et al. 2008; Comparot-Moss and Denyer 2009). Despite the rho event’s impact on functional and phylogenetic diversification of the grasses, the absence of genomic-scale analyses for early diverging grass lineages (pre-divergence of BOP and PACMAD clades) and nongrass graminids has left the precise placement of rho uncertain (Paterson et al. 2004; Soltis et al. 2009).

Later comparative genomic analyses performed with refined computational tools elucidated an even older genome duplication in grass genomes designated as sigma that was estimated to have occurred sometime prior to the diversification of Poaceae (Paterson et al. 2009; Tang et al. 2010). Analyses of the Musa genome identified three lineage-specific WGDs within Zingiberales and determined that the rho and sigma events occurred within the Poales lineage but prior to the divergence of the PACMAD and BOP clades (D’Hont et al. 2012).

Based on synteny analysis of the rice and sorghum genomes, Tang et al. (2010) hypothesized a third paleopolyploid event prior to sigma but such an event was not inferred in the Musa genome analyses (D’Hont et al. 2012). Comparative synteny analyses of oil palm (Elaeis guineensis) and eudicot genomes have recently placed this event, designated tau, as early in the evolutionary history of monocots (Jiao et al. 2014), but after the divergence of the Alismatales (Ming et al. 2015). Clearly, polyploidy has contributed to molecular variation and evolution within grass and other monocot lineages, but a full understanding of the impact of the rho, sigma, and tau events requires more precise estimation of the timing of these genome duplications relative to branching points in the monocot phylogeny.

Among the genome-wide molecular changes that may have been influenced by polyploidy, the GC profile of protein-coding genes has been hypothesized to have influenced the synonymous substitution rate of retained duplicates (Wang et al. 2005). A bimodal distribution of GC composition among coding genes, and 5′ to 3′ gradients in GC content of coding genes have long been identified as divergent characteristics of grass genomes relative to eudicot genomes (Carels et al. 1998). There have been conflicting reports on whether these are features of an ancestral monocot genome (Clément et al. 2015) or have been independently derived in the grasses (Carels and Bernardi 2000; Kuhl et al. 2004, 2005) and other plant lineages (Escobar et al. 2011; Serres-Giardi et al. 2012). The functional implications of variation in GC content among monocot lineages have been discussed in light of associations between whole-genome GC content, genome size and architecture, and environmental conditions (Smarda et al. 2014). However, the relative importance of mutational processes, including GC-biased gene conversion (Eyre-Walker and Hurst 2001; Duret and Galtier 2009), and selective processes is controversial (Shi et al. 2006; Tatarinova et al. 2010, 2013; Serres-Giardi et al. 2012; Glémin et al. 2014; Clément et al. 2015). The debate has been fueled by the observed negative correlations between GC-content and gene length, exon number (Carels and Bernardi 2000; Wang et al. 2004; Shi et al. 2006; Tatarinova et al. 2010), and gene body methylation, together with positive correlations with recombination rates and variance in gene expression (Tatarinova et al. 2010; Serres-Giardi et al. 2012). Here, we ask whether ancient polyploidy may have contributed to the evolution of variation in GC composition among coding genes in monocot genomes.

Following the classic work of Bowers et al. (2003), we employ a phylogenomic approach to resolve the timing of rho and sigma relative to speciation events within the Poales and tau earlier in monocot history. Transcriptome data for species sampled within each family in Poales were generated and assembled transcripts included in gene trees to expand taxon sampling beyond species with available genome sequences. Species were chosen for our phylogenomic analyses to increase sampling density and uniformity across the Poales phylogeny. Broader taxon sampling is considered to increase the accuracy of phylogeny (Zwickl and Hillis 2002; Leebens-Mack et al. 2005). Gene trees were estimated for all gene families with paralog pairs mapping to syntenic blocks (i.e., syntelogs) within the rice and sorghum genomes, and the timing of duplication relative to speciation events was inferred through interrogation of the gene trees. The combination of synteny-based approach and the phylogenomic approach is the best method to assay the origins of WGD duplicates. At the same time, putative single copy gene families were used to rigorously estimate phylogenetic relationships within the Poales and related monocot lineages, improving precision over previous reconstructions with plastid sequences. Gene trees were reconciled with the inferred species phylogeny to estimate the timing of rho, sigma, and tau and the influence of these events on the evolution of key innovations. Finally, variation in GC composition across the Poales was analyzed and tested to determine if such variation relates to the relative placement of these polyploid events.

Materials and Methods

Taxon Sampling

Relationships within the order Poales have recently been investigated through analyses of gene sequences extracted from plastid genome sequences (Givnish et al. 2010), and we used phylogenies from that study to guide sampling of all major clades and families within the order. RNAs were isolated from mixed vegetative and reproductive tissues for 25 species, and transcriptome assemblies were combined with available genome data sets to carry out the analyses described below. Sample identities and data descriptions can be found in supplementary table S1, Supplementary Material online along with doi numbers for RNA Seq reads, assemblies, multiples sequence alignments, and gene trees deposited in NCBI's Short Read Archive (SRA) (PRJNA313089) database and DRYAD (XXXXX).

RNA Isolation and Sequencing

RNA was isolated from fresh young leaf or apical meristematic tissue using the RNeasy Plant Mini Kit (Qiagen, Valencia, CA). Samples were kept on liquid nitrogen prior to isolation. RNA was eluted into a final volume of 100 µl of RNase-free water. RNA total mass and quality were estimated using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA). Samples were deemed acceptable if RIN (RNA Integrity Number) scores were greater than 8.0.

The TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA) was used to construct paired-end libraries with an average fragment length of 300 base pairs (bp). Sample libraries were sequenced on a GAIIx Genome Analzyer (Illumina), at University of Missouri-Columbia’s DNA Core Facility, and on a HiSeq 2000 (Illumina) at the BGI Americas sequencing facility in Davis, CA.

Transcriptome Assembly

Illumina data generated in this study and downloaded from the NCBI SRA were assembled in the same fashion. FastQC v0.10.1 (Andrews 2010) was used to identify adapter contamination and Cutadapt v1.1 (Martin 2011) was used to clean contaminated sequences. Two mismatches per 10 bp were allowed and a minimum overlap of at least 10 bp was required for contaminant identification and trimming. Cleaned reads were further trimmed using a custom perl script (https://github.com/mrmckain/poales-polyploidy, last accessed March 22, 2016) that trimmed reads from the ends until there were three consecutive bases with a quality score of 20 or more. Finally, reads with a median quality score less than 22, more than three uncalled bases, or trimmed lengths less than 40 bp were removed. We retained the orphaned as well as paired reads for downstream assembly and analyses.

Cleaned and filtered data sets were assembled using Trinity (Release 2012-06-08; Grabherr et al. 2011) with default parameters. Reads were aligned to the Trinity assembly using bowtie v0.12.8 (Langmead et al. 2009) through the alignReads.pl script available in the Trinity distribution. Output from this script was then piped into the run_RSEM.pl script (also packaged with Trinity), which utilizes RSEM v1.2.0 (Li and Dewey 2011) to quantify transcript abundance. Per kilobase of exon per million fragments mapped (FPKM) was estimated for each component (gene), and the percentage of mapped fragments corresponding to each isoform within a component was estimated using the summarize_RSEM_fpkm.pl script packaged within Trinity and an average fragment length setting of 300 bp. Isoforms that had 1% or less of all fragments mapped to a component were removed to create a data set that represented well-supported transcripts. Roche 454 transcriptome reads for E. guineensis (GenBank accession numbers: SRX059258–SRX059263) were assembled with MIRA (Chevreux et al. 2004) using default parameters.

Transcriptome Translation and Gene Family Circumscription

Translation of the transcriptome assemblies was conducted using a series of custom perl scripts (https://github.com/mrmckain/RefSeq, last accessed March 22, 2016). TBLASTX matches to a database of gene models from 22 genomes used in the gene family circumscription developed by the were identified for each assembled transcript, and matches with an e-value less than 1e-10 were used to guide translations with GeneWise (part of the Wise2 v2.2.0 package; Birney et al. 2004). The longest GeneWise translations for each transcript were used. In rare cases where internal stop codons were identified, they were spliced out of the transcript assemblies. BLASTP searches against the 22-genome amino acid database were then performed with the inferred amino acid sequences in order to assign them into the 22-genome orthogroups (i.e., gene families circumscribed by OrthoMCL (Li et al. 2003) as in the Amborella genome paper ().

Species Tree Estimation

Species trees were estimated using concatenated and coalescence-based approaches. A set of 970 single/low copy gene families were identified in our orthogroups (Duarte et al. 2010) and used as a starting set for identifying strictly defined single copy genes in our taxon set. Peptide sequences within these orthogroups were aligned using MUSCLE v3.8.31 (Edgar 2004), and coding sequences (CDSs) were aligned onto the amino acid alignments using PAL2NAL v13 (Suyama et al. 2006). Alignments were filtered using two criteria. First, columns in alignments were removed if gaps were observed in more than 90% of sequences (rows in the alignment matrix). Second, transcript translations (rows) were removed if they covered less than 30% of the total alignment length for the gene family. Maximum likelihood (ML) gene trees were estimated using RAxML v7.3.0 (Stamatakis 2006) with the general time reversible + Γ substitution model and 500 bootstrap replicates. Gene trees were rooted to outgroup taxa (Amborella trichopoda and Vitis vinifera) found in each gene family. Orthogroup trees were analyzed for gene copy count for each species using a modified version of the clone-reducer script from Estep et al. (2014) (https://github.com/mrmckain/clone_reducer, last accessed March 22, 2016) as used in Cannon et al. (2015).

The primary goal of the species tree estimation was to identify relationships of major lineages (i.e., branches coming off of the backbone of the phylogeny) rather than more shallow relationships. Putative isoforms for a species were collapsed into a consensus sequence if they formed a species-specific clade in the gene tree with a BSV (bootstrap support value) of 50 or greater. Collapsing close isoforms (possible splice variants, alleles, or recent duplicates) allowed us to maximize the number of single copy orthogroups. Further, for the species tree analysis we reduced certain lineages to a single species representative in order to avoid coalescence issues that were potentially lowering the total number of single copy orthogroups. This included reducing the genera with multiple samples (Cyperus, Typha, Juncus, and Mayaca) and the family Centrolepidaceae to either single species or combining them into a single sample. These groups were identified as the hotspots where coalescence of putative isoforms was occurring within the reduced lineage but not always within the species. We repeated the clone reducer methodology after these simplification steps and identified 234 high occupancy orthogroups that maintained genes in single copy across all taxa.

To estimate the concatenation-based species tree, each aligned gene family was concatenated for a species. If the species did not have a representative in the gene family, the taxon was represented by Ns in the alignments for the missing gene. All 234 orthogroups were represented in the concatenated analysis. The topology was reconstructed using RAxML with the GTR + Γ model over 500 bootstrap replicates and rooted relative to A. trichopoda.

The coalescence-based analysis was conducted using ASTRAL v.4.7.6 (Mirarab et al. 2014). Trees that were used in the ASTRAL analysis were generated in RAxML as above. The bootstrapping option of ASTRAL was used for 100 replicates.

A consensus topology was created to represent the shared relationships of both methodologies. Nodes were collapsed where the topology disagreed. This did not affect the estimation of the timing of WGD events.

Ks Frequency Plot Estimation and Estimation of Ks-Derived Putative Paralogs

K frequency plots were estimated using the FASTKs pipeline (https://github.com/mrmckain/FASTKs, last accessed March 22, 2016). Translated transcriptome data sets were blasted against themselves using an e-value cutoff of 1e-40. Putative pairs were filtered if they: 1) exhibited 100% identity in their alignment, 2) had less than 300 base pairs of an alignment length, and 3) had less than 40% identity. Amino acid sequences for putative paralog pairs were then aligned using MUSCLE, and back translated to CDS using PAL2NAL. K, K, and K/K were estimated for the aligned pairs using codeml in PAML v.4.8 (Yang 2007) using the same parameters used by McKain et al. (2012). Normal mixture models were estimated for K values using the mclust v.5.0.2 (Fraley and Raftery 2002; Fraley et al. 2012) as implement in R. K plots were used as a secondary confirmation of rho and sigma placement in the species tree (supplementary fig. 3, Supplementary Material online).

Estimation of Syntelogs and Gene Tree Estimation for Multicopy Gene Families

Synteny blocks and retained duplicate genes derived from the rho, sigma, and tau WGDs (i.e., syntelogs) were identified in the Oryza sativa and Sorghum bicolor genomes as described by Tang et al. (2010). The rho synteny blocks were identified using a chaining distance of 40 genes. Following this step, we reconstructed the “pre-rho” putative ancestral regions (PARs; Tang et al. 2010) by interleaving the genes between the rho duplicate pairs. A second synteny search was then performed on the pre-rho PARs to identify the sigma + tau synteny blocks using a relaxed chaining distance of 60 genes.

Jiao et al. (2014) inferred tau synteny blocks through comparisons of PARs estimated for oil palm (E. guineensis) and the eudicot, sacred lotus (Nelumbo nucifera). Each of these genomes was known to have included at least one round of paleopolyploidy since divergence and Jiao et al. (2014) estimated PARs for these two genomes using the methods described above. Tau synteny blocks were identified as pairs of oil palm PARs lining up with a single sacred lotus PAR. Syntelog sets identified in the oil palm tau synteny blocks were then compared with rice and sorghum genes to infer a total of 2,248 duplicate gene pairs in the sorghum and rice genomes derived from the tau WGD (Jiao et al. 2014). We performed phylogenomic analyses on orthogroups with theses tau duplicates and compared them with the sigma + tau syntelogs we inferred in analyses of the rice and sorghum pre-rho PARs (see above).

In some instances, syntelog pairs were sorted to different gene families as circumscribed by OrthoMCL (Li et al. 2003; ). Because we were interested in identifying the last common ancestral (LCA) node in gene trees for each syntelog set, these orthogroups were combined for alignment and phylogenetic analysis. Gene family alignment and subsequent tree estimation were conducted as described above for species tree estimation. If a gene tree did not contain a nonmonocot outgroup (Amborella or Vitis) homolog, it was dropped from further analysis.

Queries of the resulting gene tree analyses were conducted using PUG (https://github.com/mrmckain/PUG, last accessed March 22, 2016) with the consensus species tree as the guide tree. PUG uses an algorithm that queries putative paralogs in a gene tree identifying species congruence relative to the species tree. Taxa that were removed prior to species tree estimation were inserted back into the tree for estimation of WGD placement. The timing of duplication events relative to the speciation events shown in the consensus species tree (fig. 1) was estimated by querying the rooted gene trees for the species represented in clades defined by an LCA node for each syntelog set. Syntelog sets with poorly supported LCA nodes demonstrating BSVs less than 50 were not used to infer the timing of WGDs. For the remaining syntelog sets, placement of duplication events was based on the species composition within the clade defined by the LCA node in the gene trees and the species composition within its sister lineage. This second criterion for placing duplication events was implemented to minimize the impact of incomplete gene sampling in the transcriptome assemblies. For an LCA node placement to be accepted, the following criteria had to be met: 1) a minimum of two taxa had to be present above the duplication node (i.e., at least one more taxon in addition to rice or sorghum syntelog pairs), 2) no taxa found above the hypothesized WGD position in the species tree could be found in the lineage sister to the putative LCA defined clade in the gene tree, and 3) the lineage sister to the LCA defined clade had to contain at least one taxon from the sister lineage to the hypothesized WGD node in the species tree. Gene tree topologies were also inspected manually to verify the results of the automated analysis. Counts of LCAs (Last Common Ancestor for hypothesized WGD clade) with well supported placement were made for each node in the species tree and split into two categories based on the BSV of the LCA: 1) 50 ≤ BSV < 80 and 2) BSV ≥ 80.

F

Consensus species tree from concatenated and coalescence-based analysis of 234 single copy orthogroups with results from gene tree querying of putative paralogs. (A) Mapping results of synteny-derived paralogs from the rice and sorghum genomes displayed as total number of unique last common ancestor (LCA) nodes with BSV ≥80. Results show placement of rho, sigma, and tau WGD events. (B) Mapping results of K plot-derived paralogs with 22 (number for sigma event) or more total unique LCAs and boostrap values ≥80 for Poales species only. (C) Mapping results of gene tree-derived paralogs with 235 (number for sigma event) or more total unique LCAs and boostrap values ≥80 for monocot species only. Previously published WGD events are identified and placed on the tree, including a shared Zingiberaceae event (gamma), a palm event, and an Agavoideae event, represented as gold diamonds. If previously named, the Greek character representing the event is also displayed. Higher support for rho, sigma, and tau is identified relative to the synteny-derived paralogs and other potential WGD events in Juncus, Cyperus, and Restionaceae are also identified.

For each orthogroup, all possible pair combinations of transcripts/annotated genes for each species, except Amborella and Vitis, were estimated. These pairs were used in a separate analysis with PUG to identify putative WGD events by accumulated signal at all nodes.

GC Content Analysis

GC content was analyzed for all orthogroups with identified putative paralogs pairs and for other orthogroups with at least ten taxa present. A total of 13,798 orthogroups were queried. The PAL2NAL alignments were used for GC estimation. For each orthogroup, sequences were removed if they were not within 1 standard deviation (SD) of the average length of the orthogroup alignment in order to better assess 5′ to 3′ gradients in GC content along complete or nearly complete transcript assemblies.

Total GC percentage was calculated for each transcript along with 5′ to 3′ gradients in GC composition. For each orthogroup, the alignment length was estimated as number of codons and GC composition were estimated for nonoverlapping windows corresponding to 1% of the alignment. For each species, the GC content of each nonoverlapping 1% interval was calculated as the total GC percentage for that window across all orthogroups.

Mean and SD of GC composition (both total and GC3) were estimated for each taxon. A student’s t-test was used to determine if the mean total GC percentage was different from the mean GC3 in each taxon. Hartigan’s dip test (HDT) for unimodality, as implemented in diptest v.0.75-7 in R, was used to test if the distribution of total GC content and GC3 for genes across all taxa exhibited unimodality. For taxa where the unimodality null hypothesis was rejected, an inspection of the distribution of %GC across all orthogroups suggested the data were bimodal. A kmeans clustering analysis of these seven taxa using two clusters identified cluster centers at 46.7% and 63.7% GC for total GC and 45.8% and 82.5% for GC3. Bimodal distributions were particularly evident for GC3, with 18 taxa exhibiting bimodality with kmeans centers at 45.7% and 60.6% GC for total GC and 44.4% and 76.1% for GC3.

For each of the seven taxa found to have a bimodal distribution for total GC composition, kmeans clusters (n = 2) were calculated. Clustered transcripts were used to identify abundance of low GC and high GC transcripts in orthogroups. Orthogroups were then classified into three categories: 1) high GC—75% or more of transcripts for a given taxon in an orthogroup are clustered as high GC, 2) low GC—75% or more of transcripts for a given taxon in an orthogroups are clustered as low GC, and 3) mixed GC—remaining orthogroups that do not fall into other classes. Results from each taxon were compared and sets of true high GC and true low GC were identified across all bimodally distributed taxa. GC composition distributions for all transcripts in both high and low GC classes were compared. We then tested for significant differences in high vs. low GC classification distributions for each taxon and all taxa using t-tests. A contingency test was used to test for significant associations between orthogroup GC content classifications and the retention of tau, sigma, or rho duplicates across orthogroups.

Results

Transcriptome Assemblies, Translations, and Gene Family Sorting

In total, over 1.04 million transcripts from 35 species (an average of 36,758 transcripts per taxon) were assembled after adaptor trimming and quality filtering. Transcript assemblies with an average length of 1,005 bp (supplementary table S2, Supplementary Material online) were translated, and sorted into gene families as circumscribed through OrthoMCL clustering of gene models from 22 sequenced and annotated plant genomes including rice and sorghum (Li et al. 2003; Wall et al. 2008; ). Resulting orthogroups (i.e., OrthoMCL clusters) including sequences from up to all 35 species included in the study (supplementary table S1, Supplementary Material online) were used for species tree estimation and resolution of timing of ancient WGDs.

Single Copy Gene Family Phylogeny Analysis

Multiple sequence alignments for all 234 orthogroups found to have no more than one gene copy in any of the sampled genome or filtered transcriptomes were concatenated and an ML tree was estimated from the resulting super-matrix (Materials and Methods; supplementary fig. S1, Supplementary Material online). A coalescence-based species tree estimate was generated from the 234 orthogroup trees using ASTRAL (Mirarab et al. 2014; supplementary fig. S2, Supplementary Material online). Minimal discordance was found between the super-matrix and ASTRAL trees, and all relationships with bootstrap values (BSV) greater than 90% were considered reliable and recovered in both analyses (supplementary figs. S1 and S2, Supplementary Material online). These robustly estimated relationships (fig. 1) are largely consistent with previous ML super-matrix analyses of 81 plastid genes (Givnish et al. 2010; Barrett et al. 2016) with a few exceptions. The basal lineage of Poales differs between the nuclear and plastid phylogenies. Although the plastome analyses placed Bromeliaceae sister to the rest of Poales, our nuclear gene analyses recovered strong support for Typhaceae as sister to the rest of Poales in both the concatenated and coalescence-based analysis. Further, Givnish et al. (2010) found strong support for an Eriocaulaceae + Mayacaceae clade and Xyridaceae separate, whereas our analyses suggest Lachnocaulon (Eriocaulaceae) and Xyris (Xyridacae) form in a clade that does not include Mayaca (Mayacaceae). However, the support for a Xyridaceae + Eriocaulaceae clade is very low in the coalescence tree (BSV 44), suggesting rapid diversification at this point within the Poales phylogeny. Finally, we place Joinvillea (Joinvilleaceae) and Ecdeiocolea (Ecdeiocoleaceae) in a clade, sister to the Poaceae (fig. 1) in contrast to the plastome-based resolution of Joinvilleaceae and Ecdeiocoleaceae as successive sister lineages in a grade leading to the Poaceae (Givnish et al. 2010; Barrett et al. 2016). The plastid genome is inherited as a single-locus, so differences between the trees inferred from multiple nuclear genes and plastid genome may be due to incomplete sorting of ancestral variation in the plastid genome haplotypes between speciation events (Maddison 1997). In any event, the discrepancies between our inferred species phylogeny and plastome-based phylogenetic inferences involve single rearrangements across short internodes on the backbone of our species tree estimate. Therefore, discordance among gene histories due to incomplete lineage sorting is expected. We used the well-supported relationships found to be concordant in our nuclear gene super-matrix and coalescence-based analyses (fig. 1) to guide gene tree queries and estimate the timing of the rho, sigma, and tau WGD events, after the reinsertion of removed taxa (Typha angustifolia, Cyperus alternifolius, Aphelia sp.) and the splitting of combined genera (Mayaca—Mayaca fluviatilis and Mayaca sp., Juncus—Juncus effusus and Juncus inflexus).

Delineation of Syntelog Pairs from Sorghum and Oryza Genomes

Following the methodology of Tang et al. (2010), syntenic regions were identified within the genomes of S. bicolor (sorghum) and O. sativa (rice) using a hierarchical approach (see Materials and Methods). A total of 56 syntenic block pairs identified in the rice genome were assigned to the rho event, and 39 rho block pairs were identified in the sorghum genome. These blocks included 4,296 rice syntelog pairs and 3,971 sorghum syntelog pairs. Synteny analysis of the pre-rho blocks resulted in circumscription of 58 sigma + tau blocks in rice and 63 blocks in sorghum. These inferred syntenic blocks associated with sigma or tau include 1,782 and 1,898 syntelog pairs in rice and sorghum genomes, respectively. In addition, a total of 2,248 duplicated rice and sorghum gene pairs were associated with the tau event by Jiao et al. (2014), to add to our exhaustive list of gene pairs for dating purpose. All inferred rho, sigma, and tau syntelog pairs were included in gene tree analyses in order to estimate the timing of duplication relative to speciation events.

Querying of Gene Trees for Placement of WGD Events Using Paralogs from Genomic Data

Orthogroup membership was determined for all 8,267 putative rho and 3,680 sigma + tau syntelog pairs derived from our synteny analysis. In 1,680 instances syntelog pairs were placed in separate orthogroups. In these cases, the homologous orthogroups were combined before sequence alignment and tree estimation. In total, alignments and ML trees were estimated for 2,795 OrthoMCL circumscribed gene families () including one or more syntelog pairs, 2,089 of which included at least one V. vinifera or A. trichopoda gene that could be used as outgroup sequences to root gene duplications that occurred at any point in monocot history. Collectively, these 2,089 gene trees included a total of 6,484 syntelog pairs mapping to 3,203 unique ancestral nodes representing the LCA homeologous gene pairs.

The 3,203 LCA nodes were filtered based on BSV and grouped into two sets, those with 50≤ BSV<80 and those with BSV≥80 (designated as BSV 50 and BSV 80, respectively). Additionally, trees were queried in relationship to the estimated species tree, and only nodes in the gene tree congruent in topology to the species tree were accepted (see Materials and Methods). After filtering based on support values, 836 unique LCA nodes for 1,348 rho block syntelog pairs and 124 LCA nodes for 318 sigma + tau block pairs for were retained.

The estimated timing of duplication events relative to the speciation events depicted in figure 1 were inferred by determining the collection of species represented in clades defined by the syntelog LCA nodes described above. This analysis strongly suggested that the rho WGD occurred on the branch leading to the last common ancestor of all extant grass species including the basal grass Streptochaeta. Trees with rho syntelog pairs included 411 rho LCA nodes with BSV 80 defining clades with Streptochaeta genes but no genes from taxa outside of the Poaceae (fig. 1). Another 107 rho syntelog LCA nodes with BSV of 50 or better showed this same pattern. In contrast, just 123 rho LCA nodes (BSV 80; 143 BSV 50) suggested that rho occurred after Streptochaeta diverged from the rest of the grasses but prior to the diversification of the BOP + PACMAD clades (fig. 1).

As expected, sigma + tau block LCA nodes mapped to two points on the species tree corresponding to the two separate WGD events. Presumed sigma related duplication LCAs were concentrated (26 BSV 80; 11 BSV 50) on the branch leading to the last common ancestor of all extant species within the Poales, whereas the concentration of LCAs mapping to the Asparagales + commelinid clade (50 BSV 80; 14 BSV 50) is interpreted as representing tau (fig. 1).

This timing of the tau WGD was also inferred by querying orthogroup trees including the homeologs identified by Jiao et al. (2014). Of the 2,248 putative tau pairs from that study, we identified 546 pairs (416 BSV 80; 130 BSV 50) representing 56 unique LCAs. Of these 56 LCAs derived from the Jiao et al. study, 18 overlap with LCAs we identified in the sigma + tau synteny blocks, and all map to the lineage leading to the last common ancestor of Asparagales and the commelinids. This finding is in agreement with results recently reported in the Phalaenopsis equestris (Cai et al. 2014) and pineapple genome (Ming et al. 2015) papers.

In summary, phylogenomic analyses of rho, sigma, and tau syntelog pairs improve precision in the placement of each of these WGD events inferred from analyses of Poaceae genomes and provide a “gold standard” secondary source of WGD estimation.

Querying of Gene Trees for Placement of WGD Events Using Paralogs from Transcriptomic Data

In addition to using the paralogs generated from synteny analyses of rice and sorghum, putative paralogs were estimated from K analyses of all of the transcriptomes used in this study (supplementary fig. S3, Supplementary Material online). A total of 20,900 K-derived putative paralogs were queried against orthogroups identified as having rho or sigma synteny-derived putative paralogs. Many of the K-derived putative paralogs were estimated as isoforms by Trinity, and are more than likely true isoforms or alleles especially for the low Ks pairs. Of all putative paralog pairs, 667 (550 BSV 80; 117 BSV 50) representing 343 unique duplicate LCAs (274 BSV 80; 69 BSV 50) were placed within the species tree after filtering. Concentrations of gene duplications placed at the base of the Restionaceae (unique LCAs: 102 BSV 80; 27 BSV 50), at the base of the restiids (unique LCAs: 15 BSV 80; 2 BSV 50), at the base of Poaceae (unique LCAs: 78 BSV 80; 16 BSV 50), and after the divergence of Juncus from cyperids (unique LCAs: 29 BSV 80; 2 BSV 50; fig. 1, supplementary figs. S3 and S4, Supplementary Material online). In sum, the K plots and phylogenomic analyses provide consistent evidence for three additional WGD events placed on lineages leading to the restiids (including Centrolepidaceae, and Restionaceae), the Restionaceae and Juncus (supplementary fig. S4, Supplementary Material online). As discussed by Cui et al. (2006) and others, K plots may not detect ancient genome duplications so these should be interpreted as a minimal set of polyploidy events among the investigated Poales lineages.

A total of 1,870,214 unique combinations of transcripts or annotated genes were identified across all monocot taxa in sampled orthogroups. This number is inflated to what may be considered paralogs due to existence of isoforms and alleles in transcriptome data. Of these, 36,567 (BSV 80; 51,727 BSV 50) pairs were identified to correspond to nodes in the species tree after filtering. These pairs represent 5,455 unique LCAs (BSV 80; 7,107 BSV 50) spanning almost every node of the species tree. We used the lowest represented known event (sigma, 235 unique LCAs represented in this analysis) to set a minimum threshold for identification of WGD events from the gene tree-derived pairs. This threshold allows for identification of a number of events known from either synteny analysis or other publications including: tau (410 BSV 80), sigma (235 BSV 80), rho (610 BSV 80), Zingerberales gamma event (377 BSV 80)(D’Hont et al. 2012), palm WGD event (345 BSV 80) (D’Hont et al. 2012), and Agavoideae WGD event (615 BSV 80) (McKain et al. 2012; fig. 1). In addition to these published events, we also have strong evidence for a Juncus event (423 BSV 80), a Cyperus event (463 BSV 80), and a Restionaceae event (499 BSV 80) (fig. 1). There are two other possible WGD events that fell short of the imposed threshold: a Centrolepidaceae event (200 BSV 80) and a restiid event (184 BSV 80). It is quite possible that our threshold value is simply too conservative and these are true WGD events. Relaxing of the threshold even further would implicate possible events within the cyperids (162 BSV 80), within Bromeliaceae (128 BSV 80), and at the base of the PACMAD clade (121 BSV 80).

Analysis of Changes in GC Content Relative to WGDs

Total GC content (fig. 2) and GC composition at the third codon position (GC3) calculated for each taxon were highly correlated (table 1, supplementary fig. S5, Supplementary Material online). The highest average %GC found was in Brachypodium distachyon (56.3% GC, SD: 8.9%) and the lowest was in J. effusus (44.8% GC, SD: 6.0%) (table 1). Many of the GC composition distributions of Poales taxa appeared to be bimodal in nature (fig. 2), but a stringent HDT rejected the unimodal null hypothesis for only seven species (P ≤ 0.05): Lachnocaulon anceps, Aphelia sp., O. sativa, B. distachyon, Dendrocalamus latiflorus, Aristida stricta, and S. bicolor (table 1). Even among these seven species, there is variation in total GC composition and the degree to which the GC-content distributions are bimodal (fig. 2). When HDTs were conducted on GC3 data, 18 taxa exhibited a bimodal distribution including multiple cyperid taxa, which overall exhibit a much lower GC composition (table 1 and fig. 3). Interestingly, unimodality in total GC and GC3 could not be rejected for Streptochaeta, suggesting a possible loss of this characteristic in that taxon (fig. 2) or an independent trend toward bimodality in the BOP + PACMAD clade.

F

Heatmaps depicting general trends in gene GC composition across all sampled taxa. (A) Distribution of total GC content for all genes across all taxa. The varied GC composition of monocots is highlighted by the differences between low GC species (i.e., Juncus) and high GC species (i.e., Aphelia or grasses). Histogram depicts the GC distribution of Sorghum bicolor demonstrating the heatmap information in a vertical format. (B) Distribution of GC content across genes for all orthogroups and species sampled. A 5′ bias for increased GC percentage is seen.

Table 1

Statistical Tests for Total GC and GC3 Composition across 13,798 Orthogroups for All Taxa Sampled

Species	GC.Mean	GC.SD	GC3.Mean	GC3.SD	Hart.GC	Hart.GC3	T-Test. GCvGC3
Mayaca sp.	0.505	0.063	0.551	0.132	0.712	0.163	0.000
Amborella trichopoda	0.463	0.040	0.461	0.084	0.997	0.822	0.171
Aphelia sp.	0.529	0.088	0.599	0.190	0.000	0.000	0.000
Aristida stricta	0.508	0.095	0.547	0.201	0.003	0.000	0.000
Brachypodium distachyon	0.563	0.089	0.665	0.185	0.000	0.000	0.000
Brocchinia reducta	0.482	0.071	0.509	0.158	0.559	0.014	0.000
Centrolepis monogyna	0.492	0.061	0.530	0.147	0.437	0.002	0.000
Chondropetalum tectorum	0.480	0.061	0.494	0.129	0.529	0.034	0.000
Cyperus alternifolius	0.459	0.049	0.455	0.096	0.918	0.087	0.000
Cyperus papyrus	0.458	0.048	0.450	0.092	0.949	0.008	0.000
Dendrocalamus latiflorus	0.500	0.085	0.534	0.189	0.009	0.000	0.000
Ecdeiocolea monostachya	0.521	0.086	0.582	0.184	0.957	0.263	0.000
Elaies guineensis	0.487	0.071	0.503	0.155	0.762	0.118	0.000
Elegia fenestrata	0.476	0.060	0.489	0.130	0.412	0.008	0.000
Flagellaria indica	0.491	0.071	0.513	0.155	0.977	0.151	0.000
Hosta venusta	0.472	0.060	0.481	0.142	0.760	0.012	0.000
Joinvillea ascendens	0.502	0.080	0.539	0.170	0.955	0.676	0.000
Juncus effusus	0.448	0.060	0.459	0.115	0.776	0.106	0.000
Juncus inflexus	0.455	0.065	0.475	0.131	0.936	0.110	0.000
Lachnocaulon anceps	0.518	0.083	0.572	0.179	0.045	0.000	0.000
Lepidosperma gibsonii	0.471	0.050	0.484	0.096	0.667	0.011	0.000
Mapania palustris	0.471	0.049	0.486	0.099	0.232	0.000	0.000
Mayaca fluviatilis	0.498	0.069	0.538	0.146	0.300	0.018	0.000
Musa acuminata	0.517	0.079	0.577	0.171	0.861	0.658	0.000
Neoregalia carolinae	0.511	0.085	0.562	0.177	0.795	0.183	0.000
Oryza sativa	0.542	0.096	0.618	0.200	0.000	0.000	0.000
Phoenix dactylifera	0.490	0.062	0.515	0.139	1.000	0.994	0.000
Sorghum bicolor	0.553	0.093	0.642	0.194	0.002	0.000	0.000
Stegolepis ferruginea	0.497	0.081	0.533	0.164	0.193	0.004	0.000
Streptochaeta angustifolia	0.506	0.083	0.552	0.183	0.830	0.113	0.000
Typha angustifolia	0.469	0.064	0.467	0.143	0.744	0.268	0.039
Typha latifolia	0.479	0.068	0.494	0.147	0.908	0.179	0.000
Vitis vinifera	0.462	0.038	0.456	0.085	0.985	0.891	0.000
Xyris jupicai	0.520	0.080	0.584	0.167	0.948	0.521	0.000
Yucca filamentosa	0.482	0.067	0.502	0.150	0.452	0.008	0.000
Zingiber officinale	0.468	0.072	0.470	0.156	0.426	0.056	0.229

Note.—Bold values represent a significant p-value of < 0.05.

F

Distributions for percent GC across genes identified to “high GC” (red) or “low GC” (blue) orthogroups for (A) all taxa, (B) Juncus effusus, and (C) Brachypodium distachyon. (A) t-test for distributions of high and low GC orthogroups across all sampled taxa suggests the two sets are distinct (P value = 0.00). (B) The distributions for Juncus effusus demonstrate high overlap of high and low GC orthogroups. t-test of these data suggests that they are distinct sets (P value = 0.00). Juncus effusus transcripts exhibit the lowest overall GC composition across all sampled taxa and transcripts assigned to otherwise high GC composition orthogroups is strikingly lower than the overall distribution of these orthogroups. (C) The distributions for Brachypodium distrachyon exhibit almost nonoverlapping GC values for high and low GC orthogroups. This difference is supported by a t-test (P value = 0.00). Brachypodium distachyon represents the highest GC percentage of all taxa sampled.

Kmeans clusters (n = 2) estimated for the seven taxa identified as having bimodal total GC composition distributions show very similar values for the low and high GC cluster averages (table 2). A total of 13,027 out of 13,798 total orthogroups were found to be either “high GC” or “low GC” across these seven taxa. We identified 6,770 low GC orthogroups and 3,662 high GC orthogroups that were consistent across the seven bimodal taxa. The remaining 2,595 orthogroups varied in GC composition assignments (high GC or low GC) among the seven taxa. Interestingly, contingency tests indicated that orthogroups with high GC content were significantly underrepresented among those retaining syntenic rho, sigma, and tau duplicates in the rice and sorghum genomes (P values 0.027–0.0001; table 3).

Table 2

Kmeans Clustering of Taxa Exhibiting Bimodal Total GC Composition Distribution

Species	Low Mean	High Mean
Aphelia sp.	0.472	0.633
Aristida stricta	0.460	0.651
Brachypodium distachyon	0.500	0.659
Dendrocalamus latiflorus	0.460	0.632
Lachnocaulon anceps	0.466	0.616
Oryza sativa	0.483	0.661
Sorghum bicolor	0.491	0.659

Table 3

Counts of Retained and Lost Paralogs in Gc Classed Orthogroups for Sorghum and Rice with Chi-Squared Test Results

Paralog Source	Event	High GC	Low GC	Mixed	Chi Sq.	P value
Synteny	Rho retained duplicate	35	279	86	80.6502	<0.00001
Synteny	Rho duplicate lost	3,627	6,491	2,509
Synteny	Sigma retained duplicate	0	24	4	14.4838	0.000716
Synteny	Sigma duplicate Lost	3,662	6,746	2,591
Synteny	Tau retained duplicate	5	53	20	18.2902	0.000107
Synteny	Tau duplicate lost	3,657	6,717	2,575
Gene trees	Rho retained duplicate	51	385	133	109.3626	<0.00001
Gene trees	Rho duplicate lost	3,611	6,385	2,462
Gene trees	Sigma retained duplicate	15	163	48	55.9048	<0.00001
Gene trees	Sigma duplicate lost	3,647	6,607	2,547
Gene trees	Tau retained duplicate	75	235	95	19.256	0.000066
Gene trees	Tau duplicate lost	3,587	6,535	2,500

GC distributions of all transcripts for high and low GC orthogroups were found to exhibit normal distributions (Shapiro test; P value 0.00) that were found to be statistically different from each other (fig. 3 t-test, P value=0.00). Because orthogroups were classified as high, low, and mixed GC composition based on only the seven taxa with significant bimodal GC composition as indicated by the HDT, we further tested if the variation in GC between the high and low classes was found across all taxa sampled. For each taxon, the high and low GC distributions were found to be significantly distinct (t-test, P value=0.00; fig. 3). Intra-taxon variation between these classes varied immensely. For example, the difference in GC percentage between the two classes in B. distachyon (highest GC percentage measured) was much greater than difference estimated for J. effusus (lowest GC percentage measured) (fig. 3B and C). The distributions of the same GC classes between these two species were found to be statistically different (t-test, P value=0.00).

As has been described in previous work (Wong et al. 2002; Kuhl et al. 2004, 2005; Tatarinova et al. 2010; Clément et al. 2015), the bimodal GC content distribution observed in grasses and some other monocot genomes is associated with a 5′ to 3′ gradient in the GC-composition of coding sequences. The GC composition gradient is absent or very weak for taxa that were not identified as having multimodal GC3 distributions (fig. 2).

Discussion

Understanding the effect of polyploidy on the evolution of angiosperms requires both investigation of recent polyploids and characterization of genome evolution following ancient WGDs. Detection and phylogenetic mapping of paleopolyploidy events is a necessary first step in the characterization of genome evolution following ancient WGDs. Methods for detecting WGD events include synteny analysis (Bowers et al. 2003; Tang et al. 2007; Paterson et al. 2012), assessment of frequency distributions for synonymous substitutions (K) between duplicate genes (Lynch and Conery 2003; Blanc and Wolfe 2004; Cui et al. 2006; Vanneste et al. 2012), and phylogenomic analyses reconciling duplications in gene trees with known species relationships (Bowers et al. 2003; Jiao et al. 2011; Cannon et al. 2015; Yang et al. 2015). Each of these methods has its strengths and limitations (Bowers et al. 2003; Vanneste et al. 2012). Synteny may decay over time; substitutions at synonymous sites can become saturated, thus reducing the resolution of Ks; and gene tree estimation can be confounded by long-branch attraction artifacts and deep coalescence. These methods can be applied in concert to provide multiple lines of evidence for hypothesized WGD events (Bowers et al. 2003; Paterson et al. 2004, 2010; D’Hont et al. 2012; McKain et al. 2012; ; Jiao et al. 2014; Cannon et al. 2015). As first proposed by Bowers et al. (2003), the timing of WGDs identified through synteny analyses can be resolved through phylogenomic analyses including transcript sequences from taxa without reference genome sequences. Such a phylogenomic approach is becoming more reliable as the increasing availability of transcriptome data from taxa well distributed across the tree of life can be used to ameliorate estimation artifacts caused by model misspecification and long-branch attraction (e.g., Leebens-Mack et al. 2005). Further, as we show in this study, increasing the density of taxon sampling within an organismal phylogeny can yield improved precision in estimates of the timing of WGD events.

Assessing the Impact of WGD Events in Monocot History

The Evolution of GC-Profiles in Monocot Genes

The bimodal distribution of GC composition among protein coding genes in cereal grasses (Poaceae) has been studied for nearly two decades (Carels et al. 1998), and Wong et al. (2002) hypothesized that gene and genome duplications may contribute to this pattern. There has been much discussion about the processes responsible for bimodal GC composition distributions (Carels and Bernardi 2000; Wang et al. 2004; Shi et al. 2006; Tatarinova et al. 2010; Serres-Giardi et al. 2012; Tatarinova et al. 2013; Glémin et al. 2014; Clément et al. 2015), but experimental tests of mechanistic hypotheses remain intractable. Inferences have been drawn mainly from correlations between GC-composition and gene length, exon number, gene expression level, synonymous substitution rates, gene body methylation, and local recombination rates (Carels and Bernardi 2000; Wang et al. 2004; Shi et al. 2006; Tatarinova et al. 2010, 2013; Muyle et al. 2011; Serres-Giardi et al. 2012; Glémin et al. 2014; Clément et al. 2015). Our findings support and extend these empirical findings in an evolutionary context. With respect to polyploidy, we find that genes and gene families exhibiting high GC content tend not to retain homeologs following WGD events. Therefore, due to a bias toward low GC in retained genes, whole gene duplications may actually reduce frequency of high GC genes relative to low GC genes. At the same time, some of the grass species exhibit the most bimodal GC content distribution, suggesting that within these lineages any reduction has been offset by some other process. Our data suggest that this bimodality is driven by an increase in the average GC composition of a distinct group of gene families that are consistently classified as high GC across monocot lineages exhibiting bimodality in GC content distribution. These gene families appear to be ancestrally relatively high in GC content, even in species with low total GC composition (fig. 3).

Greatly improved sampling within Poaceae and Poales reveals multiple reductions and increases in the biomodality of GC composition from what appears to be an ancestrally weak bimodal GC distribution (fig. 2). Within Poaceae, GC content distributions for Aristida, and Dendrocalamus exhibit a reduction of high GC content genes, and the basal grass lineage, Streptochaeta, has a unimodal GC content distribution (fig. 2). Sampled species within the cyperid clade (including Cyperaceae and Juncaceae) exhibit weak bimodal or unimodal GC content distribution (fig. 2 and table 1). In agreement with recently published work on strict ortholog sets (Clement et al. 2015), our analyses of homologs within gene families indicate that high GC composition is not randomly distributed across coding genes, but rather it is a presumably an inherited feature of specific gene families. As Clement et al. (2015) inferred for banana, palms, and yam, the declines in GC content in the lineages described above were due to reduction of GC composition in ancestrally high GC gene families. These observations should motivate more focused comparisons of the high GC gene families we have identified in monocot lineages with contrasting GC composition frequencies. Given observed associations between high GC content and higher recombination rates, lower synonymous substitution rates, and absence of gene body methylation in species that do exhibit bimodal GC content distributions, analysis of these features in Aristida, Dendrocalamus (and perhaps all bamboos), Streptochaeta, and the cyperids could be informative and provide insight into processes governing increases and decreases in GC composition of these specific gene families.

The Evolution of Biosynthetic and Developmental Gene Networks

A primary objective of this study was to determine whether rho occurred prior to the diversification of Poaceae or within the family prior to the diversification of the BOP + PACMAD clade (Soltis et al. 2009). An indicator of the impact of rho on the evolution of grasses would be the relationship between duplicated genes, their origin, and their involvement in novel phenotypes. The starch biosynthesis pathway in the endosperm of grasses is unique in that ADP-glucose, a major component of the pathway, is synthesized not only in the plastid like other angiosperms, but in the cytosol as well (Comparot-Moss and Denyer 2009). This alternative pathway dominates over the ancestral plastid pathway for ADP-glucose production (Beckles et al. 2001), suggesting that it is evolutionarily advantageous. The derived, cytosolic pathway is controlled by genes that were duplicated in concert in Brachypodium, Oryza, Setaria, Sorghum, and Zea, placing the origin prior to an ancient speciation event in the BOP + PACMAD clade (Wu et al. 2008; Li et al. 2012). The starch biosynthesis pathway has not been characterized in basal Poaceae species (including Streptochaeta) or the other graminid families—Ecdeiocoleaceae, Joinvilleaceae, and Flagellariaceae—but placement of the rho WGD in the last common ancestor of all extant Poaceae lineages leads to the prediction that basal lineages in the family (e.g., Streptochaeta) potentially harbor the genes necessary for cytosolic production of ADP-glucose (Wu et al. 2008), while members of the other graminid families would not.

Placement of the rho WGD also elucidates the molecular basis for the evolution of morphological characters. The spikelet serves as a key characteristic of the grasses, though the evolution of the structure has long been a subject of research due to the similarity to closely related graminids (Rudall and Stuppy 2005; Sajo and Rudall 2012) and difficulty in determining homology in inflorescence structure between early diverging Poaceae species and the core Poaceae (Sajo et al. 2008; Preston et al. 2009; Sajo and Rudall 2012). The influence of the rho WGD event on the development of the spikelet is suggested through analyses of the MADS-box transcription factor gene family AP1/FUL, which demonstrated that the FUL gene was duplicated prior to the diversification of Poaceae (including Streptochaeta) but after the divergence of Joinvilleaceae (Preston and Kellogg 2006). Interestingly, the paralogs were not maintained in duplicate in either Streptochaeta or Pharus, both of which are early diverging lineages with distinct inflorescence structures compared with core Poaceae (Preston and Kellogg 2006). These findings, in concordance with the placement of rho in the lineage leading to the Poaceae clade, lead us to hypothesize that the retention of duplicate genes led to the development of the spikelet in grasses, though more validation at the functional level looking at paralog retention and expression in the inflorescence across the Poaceae is required to elucidate this complex question.

Synteny analyses (Tang et al. 2010) and a previous phylogenomic analyses including palm and banana genes (D’Hont et al. 2012) suggested that sigma occurred sometime prior to the diversification of Poaceae but after the divergence of Poales from the other commelinid orders. Divergence time analyses suggested that this event occurred approximately 130 Ma (Tang et al. 2010), which would be prior to the diversification of the Poales crown group over 110 Ma (Magallón and Castillo 2009; Magallón et al. 2015) and the divergence of the Poales lineage from other commelinid orders (∼110 Ma; Magallón et al. 2015).

The impact of the sigma event on the evolution of Poales may be seen in the diversity of the group. The order contains approximately 21,000 species, representative of approximately 33% of all monocot species, and is ecological dominant in a number of habitats (Linder and Rudall 2005; Givnish et al. 2010). Of the approximately 21,000 Poales species, three Poales families—Poaceae (∼11,300), Cyperaceae (∼5,700), and Bromeliaceae (∼3,100)—contribute approximately 20,000 species (Givnish et al. 2010). Though the number of species in Poaceae may be partly attributable to the rho event due to diversification in the BOP + PACMAD clade, the other diversifications could be linked to innovations spurred by gene duplications from the sigma event, such as the evolution of the epiphytic habit in Bromeliaceae (Givnish et al. 2011). Investigation into the evolutionary effects of the sigma event will require further identification of duplicated genes and gene networks derived from the event and deep genomic sampling across Poales.

Though this study was primarily aimed at the phylogenomics of Poales and the placement of the rho and sigma WGD events, we were able to detect the tau WGD event characterized by Jiao et al. (2014) and verify its placement prior to the divergence of Asparagales and the commelinid lineage. Analyses within the pineapple genome paper, are consistent with ours in placing the tau WGD before divergence of the commelinid and Asparagales lineages and suggest that it occurred after the divergence of Alismatales from other monocot lineages (Ming et al. 2015). This earlier WGD in monocot history was previously hypothesized by Tang et al. (2010) and supported in our gene tree analyses by 407 unique gene families. A similarly timed duplication of MADS-box genes has been reported to be shared by the Commelinales and Poales, though sampling and support was not high enough to further estimate when exactly the duplication occurred (Litt and Irish 2003). Further sampling of the monocots is needed to more precisely place tau relative to the rapid diversification of monocot lineages early in the history of the group (e.g., Givnish et al. 2010, Magallón et al. 2015) and to understand possible implications of this event on the history of monocots.

Phylogenomic analyses are becoming much more common as high-throughput sequencing becomes the standard. These types of analyses provide insight into the evolution of large portions of the genome of many taxa not previously sampled in genomic studies. Combining phylogenomic analyses with other methods for identifying WGD events, such as synteny analysis, allows for the unambiguous detection of paleopolyploid events and their phylogenetic placement with improved precision and confidence. Understanding when these events occurred in the history of angiosperms will help elucidate the long-term evolutionary patterns associated with polyploidy and further the understanding of how this group of organisms has become so widespread and diverse.

Supplementary Material

Supplementary tables S1 and S2 and figures S1–S5 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).