Phylotranscriptomics Resolves the Phylogeny of Pooideae and Uncovers Factors for Their Adaptive Evolution

Abstract Adaptation to cool climates has occurred several times in different angiosperm groups. Among them, Pooideae, the largest grass subfamily with ∼3,900 species including wheat and barley, have successfully occupied many temperate regions and play a prominent role in temperate ecosystems. To investigate possible factors contributing to Pooideae adaptive evolution to cooling climates, we performed phylogenetic reconstruction using five gene sets (with 1,234 nuclear genes and their subsets) from 157 transcriptomes/genomes representing all 15 tribes and 24 of 26 subtribes. Our phylogeny supports the monophyly of all tribes (except Diarrheneae) and all subtribes with at least two species, with strongly supported resolution of their relationships. Molecular dating suggests that Pooideae originated in the late Cretaceous, with subsequent divergences under cooling conditions first among many tribes from the early middle to late Eocene and again among genera in the middle Miocene and later periods. We identified a cluster of gene duplications (CGD5) shared by the core Pooideae (with 80% Pooideae species) near the Eocene–Oligocene transition, coinciding with the transition from closed to open habitat and an upshift of diversification rate. Molecular evolutionary analyses homologs of CBF for cold resistance uncovered tandem duplications during the core Pooideae history, dramatically increasing their copy number and possibly promoting adaptation to cold habitats. Moreover, duplication of AP1/FUL-like genes before the Pooideae origin might have facilitated the regulation of the vernalization pathway under cold environments. These and other results provide new insights into factors that likely have contributed to the successful adaptation of Pooideae members to temperate regions.


Introduction
The success of plants depends on their adaptation to changing environments, such as the global cooling climate from the mid-Eocene ($46 Million years ago, [Ma]) to late Oligocene ($27 Ma), which included the contraction of tropics and expansion of temperate regions (Zachos et al. 2001;Eldrett et al. 2009). In response to such cooling climates, some angiosperm families (such as Saxifragaceae and Brassicaceae) have evolved the ability to adapt to seasonal climates with cold periods and are now regarded as temperate lineages (Deng et al. 2015;Huang, Sun, et al. 2016). Although such adaptation to cooling climates can be found in a number of lineages, further diversification of the corresponding taxon groups is a relatively rare accomplishment (Donoghue 2008;Watcharamongkol et al. 2018), as supported by the dramatically lower biodiversity in the temperate zone than those in the tropics (Mittelbach et al. 2007) and large-scale extinctions associated with cooling periods (Ivany et al. 2000). Furthermore, seasonal low temperature remains a major challenge nowadays even for temperate biomes. Hence, how plants have evolved the ability to thrive in temperate regions is of great interest and significance.
The estimation of the dates of origin and diversification of major clades (sublineages) of Pooideae can provide clues about evolutionary adaptation to temperate regions (Pimentel et al. 2017); however, previous molecular inferences yielded a wide range of estimated ages for the origin of Pooideae from the Cretaceous-Palaeocene boundary to the mid-Eocene ($66-46 Ma) (Prasad et al. 2011;Christin et al. 2014;Schubert et al. 2019). Nevertheless, regardless of the estimated time of the Pooideae origin, Pooideae lineages have likely experienced global cooling periods since the Eocene and diversified to a great extent in particularly large groups such as Poeae, and further become dominant in temperate region nowadays Pimentel et al. 2017;Schubert et al. 2019). Furthermore, the distribution of species richness in Pooideae varies dramatically at the tribe level , with a single tribe Poeae having three-fifths of genera (121 of 202) and species ($2,562 of 3,900 spp.). Poeae and three other tribes Bromeae ($165 spp.), Littledaleeae ($4 spp.), and Triticeae ($501 spp.) together are referred to as the "core Pooideae" (Soreng and Davis 1998), including >80% of Pooideae species diversity. Among the remaining 11 tribes, three have >100 species: Stipeae ($527 spp.), Bromeae (165 spp.), and Meliceae (158 spp.), whereas others have <30 species (Duthieeae, Diarrheneae, and Brachyelytreae), or only one species (Ampelodesmeae, Brylkinieae, Phaenospermateae, Nardeae, and Lygeeae). Hence, greater sampling is needed to gain further understanding of divergence times and diversification in Pooideae.
Temperate grasses have also evolved diverse morphologies under a variety of environments, including different inflorescences (panicle, raceme, or spike) and variation in flowers (florets) number within a relatively small and compact "branch," called spikelet, of the inflorescence (Kellogg 2015). Each floret includes bracts (usually two: lemma and palea), lodicules (2-3), stamens, and the ovary. It is possible that some of the diverse morphologies might be associated with adaptive evolution; for example, the Triticeae species (such as wheat and barley) with spike inflorescences are often distributed in regions with arid and cold environments, whereas Meliceae species (such as Glyceria) with inflorescence of the panicle type are found in wetlands with relative warm climates. Hence, analysis of evolutionary histories of these traits can help understand adaptive evolution of Pooideae; however, these have not been extensively investigated across Pooideae. Moreover, molecular studies of adaptation to cold temperatures in Pooideae cereals (e.g., wheat and barley) and relatives (e.g., Brachypodium) showed that CBF/DREB1 and COR related genes can mediate reduction of cellular chilling/freezing damage. In addition, cold-dependent flowering (vernalization) is important for adaptation to cold climates and involves functions of VRN1, VRN2, and VRN3 and other genes (Yan et al. 2003(Yan et al. , 2004(Yan et al. , 2006Sandve and Fjellheim 2010;Zhao et al. 2015;Xu and Chong 2018). These findings suggest that genetic innovation for responses of cold might have contributed to the adaptation of Pooideae lineages to temperate regions (Zhong et al. 2018), but the evolutionary history of traits and genes related to cold adaptation in Pooideae at the tribe/subtribe level has not been examined.
Gene duplication (GD) is a major process supporting functional innovations and evolutionary adaptation to environmental opportunities (Crow et al. 2006;Innan and Kondrashov 2010;Panchy et al. 2016) and can be due to whole-genome duplication (WGD), tandem duplication, and transposon-mediated duplication (Panchy et al. 2016). WGD produces massive raw materials available for subsequent genetic variation (Van de Peer et al. 2017), is very common during angiosperm evolution (Jiao et al. 2011;Ren et al. 2018;Leebens-Mack et al. 2019;Guo et al. 2020;Zhao et al. 2021), and has been associated with increased adaptive potential and organismal diversification (Arnold 1992;Fawcett et al. 2009;Tank et al. 2015;Landis et al. 2018;Wu et al. 2020). Tandem duplication can result from unequal crossing-over between similar gene copies that are positioned closely on a chromosome (Zhang 2003). Such tandemly duplicate copies are widespread in plants, contribute to gene family expansion (Cannon et al. 2004), and can exhibit variation in expression and function among copies ). In addition, transposon-mediated duplication and retroduplication have been found in grasses including rice (Jiang et al. 2004;Wang et al. 2006) and maize (Du et al. 2009). In a previous study, Sandve and Fjellheim (2010) proposed that duplications of cold related genes have played a role in Pooideae adaptation to cool climates during the global cooling period ($33 Ma). However, the number and positions of such GDs in the Pooideae phylogeny remain largely unclear.
Macroevolutionary analyses of adaptation, divergence time, morphological transitions, and gene families in Pooideae depend on a robust Pooideae phylogeny with well-resolved relationships as mentioned previously (Tank et al. 2015;Smith, Brown, Yang, et al. 2018;Wang et al. 2019). The molecular phylogeny of grasses has greatly improved in recent decades (GPWG 2001;Schneider et al. 2011;GPWGII 2012;Kellogg 2015;Soreng et al. 2017;Saarela et al. 2018). For example, an integrated analysis of molecular and morphological data of 62 taxa supported classifying Pooideae into 13 tribes (GPWG 2001) and a subsequent study with greater sampling supported the same 13 tribes and divided the largest tribe, Poeae, into Poeae 1 and Poeae 2, each with multiple subtribes (GPWGII 2012) (supplementary fig. Zhang et al. . doi:10.1093/molbev/msac026 MBE S2, Supplementary Material online). More recent reports recognized ten Pooideae tribes and 15 Poeae subtribes (Kellogg 2015), 14 tribes and 21 subtribes , or 15 tribes and 26 subtribes in Poeae ) (supplementary figs. S1 and S2, Supplementary Material online); the increases in numbers of tribes and subtribes in the latter two studies reflect differences in relationships for some genera and have relied largely on phylogenetic analyses using two plastid genes (matK and ndhF) from 448 grasses and other phylogenetic and morphological information . The studies have divided Poeae into the two clades as reported previously (Soreng et al. 1990, namely Poeae chloroplast group 1 (PCG1) and Poeae chloroplast group 2 (PCG2) (supplementary fig. S2, Supplementary Material online).
However, many of the relationships among Pooideae tribes and subtribes remain unclear or inconsistent (supplementary figs. S1 and S2, Supplementary Material online), including the relationships among the five early divergent tribes (Duthieeae, Brylkinieae, Meliceae, Phaenospermateae, and Stipeae) (supplementary fig. S1, Supplementary Material online) (Schneider et al. 2009(Schneider et al. , 2011GPWGII 2012;Soreng et al. 2015Soreng et al. , 2017. In addition to differences in the number of Poeae subtribes (supplementary fig. S2, Supplementary Material online) (Soreng et al. 2007Kellogg 2015), the PCG1 and PCG2 groups of Poeae subtribes from plastid sequence analyses are not supported by nuclear phylogenetic analyses, with some members of PCG2 nested within PCG1 as well as the paraphyly of other PCG2 subtribes (Quintanar et al. 2007;Schneider et al. 2009;Saarela et al. 2017). To clarify the macroevolutionary history of Pooideae at the tribe and subtribe levels, both taxon sampling representing tribes and subtribes and sufficient phylogenetically informative markers are needed; the latter can be achieved by using a large number of nuclear genes from nextgeneration sequencing (Li et al. 2003;Ebersberger et al. 2009;Zimmer and Wen 2015). Several nuclear phylogenomic studies have shown their unprecedented potential to uncover well-resolved species relationships (Duarte et al. 2010;Zeng et al. 2014;Huang, Sun, et al. 2016;Xiang et al. 2017;Qi et al. 2018;Cai et al. 2019;Leebens-Mack et al. 2019).
Here, we used 157 Pooideae data sets representing all tribes and 92% of Poeae subtribes to identify putative orthologous genes for the reconstruction of Pooideae phylogeny. We further conducted molecular clock analysis using the nuclear Pooideae phylogeny and historical temperature changes as references. Moreover, ancestral character analyses and GD detection provide clues to the possible factors of adaptive evolution of temperate grasses.

Taxon Sampling and Transcriptome Sequencing
We sampled a total of 157 species of Pooideae to represent all 15 tribes, and 24 of 26 subtribes of Poeae (supplementary table S1, Supplementary Material online), following the Pooideae classification (such as supertribes, tribes, and subtribes) of Soreng et al. (2017) unless otherwise specified. For the three largest tribes, the sampled taxa cover $92% of the Poeae subtribes (the largest tribe with 121 genera; two unsampled subtribes together have only six species) and nearly half of genera of each of the tribes Stipeae (13 of 28 genera) and Triticeae (13 of Saarela et al. 2015). All tribes and subtribes with at least two species (eight tribes and 18 subtribes) are also monophyletic with 100% BS except Diarrheneae, which is split into two nonsister lineages (Diarrheneae I and Diarrheneae II; fig. 1 and supplementary figs. S5-S9, Supplementary Material online). Diarrheneae I (Neomolinia, distributed in eastern Asia) is sister to Brachypodieae plus the core Pooideae, whereas Diarrheneae II (Diarrhena, native in North America) is sister to Stipodae (Stipeae and Ampelodesmeae). The relationships among all tribes are consistent and received 100% BS in all five trees, and those among subtribes of Poeae are also strongly supported, with !80% BS in 83% nodes in at least four trees (except Cinninae and Avenula pubescens; see below and supplementary note, Supplementary Material online, for details).

Relationships at the Tribe Level
In our Pooideae phylogeny, the first two divergent lineages are, respectively, Brachyelytreae and the supertribe Nardodae with the monotypic tribes Lygeae and Nardeae. These placements receive 100% BS in all five trees and agree with previous studies ( fig. 1

The Divergence Time Estimation for Temperate Grasses
To trace the Pooideae evolutionary history, we performed molecular dating using a penalized-likelihood method implemented in treePL (Smith and O'Meara 2012) with the concatenated data set of 373 OGs. Because of the uncertain placement of the fossil N13 (supplementary table S4, Supplementary Material online) (Prasad et al. 2011), we used three calibration strategies for the estimation with fossil N13 at two different calibration placements, or without the fossil. The estimated times of same node from three calibration strategies show relatively small differences (<3 My, see supplementary figs. S13 and S14 and table S5, Supplementary Material online, for details). For convenience, we describe the result of calibration 1 hereafter unless otherwise stated. As shown in figure 3 and supplementary figure S13, Supplementary Material online, the Pooideae stem lineage was dated to $68.8 Ma with a 95% confidence interval (CI) from 68.3 to 69.1; it took $20 My for the early Pooideae ancestors to diverge into four extant early divergent clades (namely Brachyelytreae, Nardodae, Duthieeae, and the remaining Pooideae). The molecular dating suggested three periods with concentrated divergences, respectively, of tribes, subtribes, and genera. From the early middle to late Eocene ($48-39 Ma), several groups diverged within $8 My ($8. The Pooideae phylogeny is shown with ages of most clades to the left of the corresponding nodes with mean age and 95% CI (in the square brackets). Names of genera (italic), subtribes (close to genera on the right, only in Poeae), and tribes (blue; on the far right) are shown to the right of branch tips. The number of species in each genus is shown in the parentheses. Squares on the left of genera indicate their climate distribution, temperate (blue), subtropical/tropical (orange), and cold (hollow), stars on the right of genera indicate nonmonophyletic genera. Geological timescale is shown below the time tree. The blue dotted line indicates the rapid reduction of global temperature at the Eocene-Oligocene boundary. The transition of character state is shown with a circle, and the number or letter within it represents the character corresponding to the illustration on the upper left. The blue vertical bands represent the concentrated periods of expansion of tribes and genera, respectively, which correspond to periods with relatively obvious temperature decline. The yellow vertical band represents the concentrated period of expansion of subtribes, with a period of temperate fluctuation after the first cooling period. The global climate curve over the last 65 My (modified from Zachos et al. [2008]) is shown at the top. Time periods of major paleoclimatic events are highlighted: PETM, the Paleocene-Eocene Thermal Maximum; EECO, the Early Eocene Climatic Optimum; MECO, the Mid-Eocene Climatic Optimum; MMCO, the Mid-Miocene Climatic Optimum.  mentary table S5 and figs. S13 and S14, Supplementary Material online). The greater CI ranges from BEAST analysis also suggest that the divergence times for some tribes (such as Diarrheneae I and Diarrheneae II; supplementary fig. S15, Supplementary Material online) could be older and outside the period of cooling temperature ($48 to $33 Ma, from the early middle to late Eocene). The greater uncertainties in the dating results as suggested by larger CI ranges from BEAST, especially those for early divergences, allow for the possibility of older origins of the Pooideae lineages and could weaken their proposed links to paleoclimates.

Transitions in Pooideae of Characters Related to the Adaptation to Cool Climates
Using our highly resolved nuclear phylogeny with a timescale (calibration 1) as a reference, we reconstructed the evolutionary histories of six characters with coding information at the genus level ( . Water limitation has been common in temperate regions and the adaptation to different water availability is crucial for the success of plant species. We found that the Pooideae ancestor was likely a mesophyte requiring moderate amounts of water, and that this state has been maintained in most lineages (supplementary fig. S18, Supplementary Material online). Furthermore, there were shifts from mesophyte to xerophyte (requiring little water), including one at the MRCA of Lygeae þ Nardeae and the other one at MRCA of most Triticeae genera (at least 11 genera here) including the Triticum-Aegilops complex, Secale and Hordeum (supplementary fig. S18, Supplementary Material online).
For the remaining three traits related to spike structure, our results showed that the Pooideae ancestor was likely with a panicle, an apical/sinus awn, and a glabrous ovary ( These changes in trait states might have contributed to the adaptation of Pooideae to cooling climates and contributed to the formation of large groups. For example, a dorsal awn could have promoted adherence of the spikelet to passing animals, which might have been important for seed dispersal in the relatively open habitats of temperate grasslands; in addition, a fur-like hairy ovary indumentum might have protected the ovary from the damage of low temperatures. Compared with previous studies, some results of our ancestral character reconstruction are consistent, but others are different: the origin of the Pooideae ancestor in temperate region with moderate requirement of water (mesophyte) is generally consistent with the results of Edwards and Smith (2010) and Watcharamongkol et al. (2018). The former inferred that the mean annual temperature was 13-16 C, and the mean annual precipitation was 1,000-1,500 mm. The latter further inferred a freezing condition in temperate regions. On the other hand, the proposed closed habitat of the Pooideae ancestor is in line with the result of Bouchenak-Khelladi et al.

Diversification Rate Analyses
To examine the history of diversity changes during Pooideae evolution, we used BAMM (Rabosky 2014;Rabosky et al. 2014) and MEDUSA (Alfaro et al. 2009) to infer the shifts of diversification rate. Here, we focused on the shift above the tribe and subtribe levels, rather than at the genus level, because the taxon sample here includes all tribes and >92% (24/26) of subtribes but just over 40% of genera (81/202 Phylotranscriptomics and Phylogeny of Pooideae . doi:10.1093/molbev/msac026 MBE rate were detected in Pooideae from the two methods. The BAMM analyses detected five upshifts, including one shift close to the MRCA of the core Pooideae (with $3/5 of genera and species of Pooideae), and four others at the MRCA of Meliceae, Stipeae, Triticeae plus Bromeae, and Poeae, respectively. Two of the tribes (Poeae and Triticeae) are the largest tribes among the core Pooideae, whereas two other tribes (Meliceae and Stipeae) are the largest outside the core Pooideae. MEDUSA revealed two upshifts: one near the MRCA of the core Pooideae, the same as one of BAMM results, and the other one at the MRCA of Phaenospermateae and Poeae with all Pooideae (including all groups affected by the five upshifts detected by the BAMM analyses), except the four small early divergent tribes. In addition, the increase in diversification rate seems to be linked to the species richness of large tribes or other large groups; this has also been noted in previous studies Pimentel et al. 2017). For more detailed comparison between this and previous studies, see the supplementary note, Supplementary Material online.

Gene Duplications and Evidence for Potential Genome Duplication Events in Pooideae
To unravel the evolutionary history of GD across Pooideae and to identify phylogenetically clustered GDs, in part as signals for potential WGD, we applied gene-tree and species-tree (supplementary fig. S23, Supplementary Material online) reconciliation as described before ) (see details in Materials and Methods and supplementary note, Supplementary Material online). GDs can result from different processes, including WGD, which would lead to a cluster of GDs mapped to a specific node on the species phylogeny. To be conservative in using GD burst as support for candidate WGD, we employed the following two criteria: 1) GD numbers !600, GD ratio !6% of relevant gene trees; 2) GDs of the ABAB type (i.e., each subclade from the GD contains the species in both two descendent lineages [A and B] from the node) !30% (see Materials and Methods). The analyses resulted in 13 GD clusters (supplementary fig. S24 and table S9, Supplementary Material online), including one (CGD1) at the MRCA of Pooideae and those (CGD2-CGD4) at the MRCA of each of three large subclades, three (CGD5-CGD7) at the MRCA of the core Pooideae or associated with the large tribes Poeae and Triticeae, and six others shared by smaller clades in Stipeae, Triciceae, and Poeae.
To test further whether the GD clusters (CGDs) detected here represent signals for WGD, we analyzed genomic synteny (spatial analysis of gene order) (Vision et al. 2000) using relevant genomes and distribution synonymous substitution rate (Ks) for paralogs (Blanc and Wolfe 2004), as previously performed for other WGDs (Vanneste et al. 2013;Jiao et al. 2014;Yuan et al. 2018). Seven of the 13 CGDs (CGD1-CGD7) did not exhibit syntenic signals among affected genomes analyses (supplementary table S9, Supplementary Material online), suggesting that these CGDs were not from WGD. Furthermore, examination of paralogs in sequenced genomes supports the idea that some of these GDs might have resulted from tandem duplication or TE-mediated duplication (supplementary tables S9 and S10, Supplementary Material online). It is worth noting that one CGD (CGD8) was mapped to the MRCA of Triticum and Aegilops; this is consistent with the hybridization events proposed in previous studies (Marcussen et al. 2014;Li et al. 2015;Gl emin et al. 2019) and discussed in the phylogeny part presented in the supplementary note, Supplementary Material online. The remaining five CGDs are further supported by Ks analyses with a peak between 0.13 and 0.14 (supplementary fig. S25, Supplementary Material online), although the corresponding taxa lack sequenced genomes for synteny analyses. These five CGDs support five candidate WGDs ( , consistent with the previous reports that genes involved in signaling pathway are more likely to be retained after duplication (Freeling 2009;McGrath et al. 2014).
The CGD5 node ( fig. 4) is also associated with an upshift in diversification rate and represents >80% species richness of Pooideae (the core Pooideae); thus, we examined the putative functions of duplicates from CGD5 in more detail, with top 24 gene families with most copies listed in supplementary table S11, Supplementary Material online. The annotated functions of these duplicated genes include responses to abiotic stress and biotic stress, and spikelet development. For example, several CGD duplicates encode homologs of the rice OsWAK11 receptor-like protein kinase, which regulates response to toxic heavy metals and inhibits copper uptake (Xia et al. 2018 (Idanheimo et al. 2014), of DCL2 that participates in antiviral silencing response to turnip crinkle virus (Xie et al. 2004), of AtZAR1 and UGT76B1 that prevent damage by pathogenic bacteria (Baudin et al. 2017;Holmes et al. 2021), of MPL1 that acts against green peach aphid (Louis et al. 2010), of OsMADS1 and ENL1 involved in regulating spikelet morphogenesis (Hara et al. 2015;Zhang et al. 2018), and of UBC26 and AtCBF2 related to abscisic acid (ABA) pathways (Yu et al. 2016;Lin et al. 2019;Fernandez et al. 2020). The increased copy number of these genes due to CGD5 suggests that enhanced stress responses might have contributed to adaptive evolution and expansion of Poeae and Triticodae.
Statistical Tests for the Correlation between GDs, Character Transitions, and Diversification Rates To investigate the association between different types of changes from this study, we performed two statistical tests, which have been used in several previous studies (Allio et al. 2021;Lai et al. 2021;Orr et al. 2021). We used the MuSSE test (Multiple State Speciation and Extinction) (Maddison et al. 2007) to investigate whether duplications have had an impact on diversification rates. The GDs were analyzed in three different combinations: five WGDs, eight CGDs, or five WGDs plus eight CGDs, with the selection of the best-fitting model (supplementary fig. S28 and table S12, Supplementary Material online). Overall, the more a lineage has experienced duplications, the higher the net diversification rates (¼speciation rate minus extinction rate) (supplementary fig. S28c, f, and i, Supplementary Material online). When considering the five WGDs, the higher net rate has probably resulted from a lower extinction rate (supplementary fig. S28b, Supplementary Material online), whereas the addition of CGDs seem to evoke higher speciation rates (supplementary fig. S28d and g, Supplementary Material online). We also performed a statistical test for the correlation between two of the three types of features (i.e., diversification rate shifts vs. trait transition; trait transition vs. GDs; GDs vs. diversification rate shifts) following Tank et al. (2015) (also see Landis et al. [2018] and

Molecular Evolution of AP1/FUL-Like Subfamily of MADS-Box Genes
One of the characteristics of the temperate climate is seasonal change, especially the prolonged low temperature through the winter. For temperate cereals and other temperate plants, the treatment with long exposure to low (winter-like) temperature is called vernalization, which leads to the flowering transition from vegetative growth to reproduction. Vernalization is required for flowering of temperate grasses and others similarly adapted to regions with cold winters, but not for flowering in (sub)tropical grasses such as rice; hence vernalization-dependent flowering was thought to be a newly evolved process in temperate lineages (Li et al. 2016;McKeown et al. 2016). Members of the APETALA1/ FRUITFULL-like (AP1/FUL-like) subfamily of MADS-box genes play roles in response to vernalization in cereals (such as VRN1 in wheat) (Yan et al. 2003(Yan et al. , 2004Trevaskis 2010;Kippes et al. 2015) and in inflorescence and floral meristem identities in rice (FUL1/OsMADS14, FUL2/OsMADS15, FUL3/ OsMADS18, and FUL4/OsMADS20) (Wu et al. 2017). Duplications of AP1/FUL-like genes in Poaceae were previously reported (Preston and Kellogg 2006;Wu et al. 2017), but the phylogenetic positions of duplications are largely unclear because of limited taxon sampling.
Here, to examine the evolution of AP1/FUL-like subfamily in grasses, especially in Pooideae, and to position the GD event(s), we selected 26 Poaceae, 15 other Poales, and seven outgroups (supplementary table S13, Supplementary Material online) to reconstruct the gene family tree. The results here ( fig. 5a and supplementary fig. S30, Supplementary Material online) indicate four grass-specific clades of AP1/FUL-like genes, consistent with the previously defined FUL1-FUL4 clades (Wu et al. 2017). Further analysis of the genomic collinearity in grass genomes indicates that the duplicated paralogs are mostly located in syntenic blocks ( fig. 5b), suggesting that the four FUL clades are a result of genome duplications. The grass family was reported to experience three polyploidizations referred to as tau (s) and sigma (r) at the MRCA of commelinid monocots and Poales, respectively, and rho (q) associated with the early Poaceae history (Paterson et al. 2004(Paterson et al. , 2009D'Hont et al. 2012;Jiao et al. 2014;Ming et al. 2015). Here, we propose that the s, r, and q WGDs together generated the four FUL clades: 1) the oldest duplication leading to FUL1-FUL2 (blue) and FUL3-FUL4 (yellow) occurred at the MRCA of commelinids including Poales, Zingiberales, and Arecales (circle 1, fig. 5a), corresponding to the position of s; 2) the event giving rise to FUL3 and FUL4 is across all members of Poales (circle 2) coincides with r; and 3) the duplication generating FUL1 and FUL2 (circle 3) affects most Poaceae subfamilies, except the early divergent subfamily Anomochlooideae, very close to q, whose precise position is still unclear.   MBE the representatives from the early divergent tribes (Nardeae, Lygeeae, Duthieeae, and Phaenospermateae); these tribes were estimated to diverge during the cooling period in the Eocene ( fig. 3), suggesting the possibility that the early Pooideae might already have had the potential to adapt to seasonal climate by vernalization. However, FUL4 homologs were only found in Poeae members (transcriptome data) and B. distachyon but not in other Pooideae. The expression levels of the FUL4 ortholog in the rice (OsMADS20) and B. distachyon (BdMADS31) are very low (Li et al. 2016;Wu et al. 2017), suggesting that the nondetection of FUL4 homologs in our transcriptome data sets might be due to their low expression levels. However, the absence of FUL4 homologs even in genome sequences here suggests that FUL4 might not be required for vernalization and might have undergone pseudogenization, consistent with the previous result that the FUL-like motif was lost in the B. distachyon FUL4 homolog (Li et al. 2016).

Multiple Tandem Duplications of CBF Genes at Different Stages of Pooideae Evolution
C-repeat binding factors (CBF) are transcription factors that belong to the APETALA2/Ethylene Responsive Factor (AP2/ ERF) protein superfamily. CBF can bind to the C-repeat/dehydration responsive (CRT/DRE) element of Cold-regulated (COR) and Dehydrin (DHN) genes and enhance the tolerance to low temperature, drought, and salt stresses (Stockinger et al. 1997;Haake et al. 2002;Nakashima et al. 2009). The numbers of CBF genes vary greatly across different lineages: from six in Arabidopsis thaliana and ten in rice  to 20 in barley (Skinner et al. 2005) and 37 in wheat (Guo et al. 2019). The copy number variation of CBF in barley population was found to associate with the level of freezing tolerance (Knox et al. 2010;Würschum et al. 2017). In addition, Sandve and Fjellheim (2010) suggested that the amplifications of CBF and FT genes might have played a role in helping the core Pooideae ancestor adapt to the super-cooling climate during the E-O transition, although this proposal was based on an analysis using only eight core Pooideae species. The recent advances in genome sequencing of Pooideae Wang et al. 2020;Zhang, Liu, et al. 2020) provide a great opportunity to discern the evolutionary pattern of CBF copy number across temperate grasses, as clues to the molecular basis of adaptation to cold climates. Here, we performed phylogenetic analysis of CBF2A, which promotes COR transcription and enhances freezing tolerance (Jeknic et al. 2014), and its homologs from 15 grass species and two outgroups (supplementary table S14, Supplementary Material online). The resulting gene tree ( fig. 6 and [8]). In contrast, the tropical/subtropical rice and bamboo and their relatives have low copy numbers of CBF genes, suggesting that increased CBF copy numbers might be associated with the plant distribution to regions of temperate climates, as proposed before for temperate cereals (Knox et al. 2010).
To gain further insights into the evolution of CBF genes in Triticeae, we examined the chromosomal location of CBF duplicates and found that they are located in a gene cluster on the same chromosome, with fewer than three non-CBF genes between any two closest CBF genes (supplementary fig. S32, Supplementary Material online). Specifically, CBF genes in Triticeae, Poeae, and Brachypodieae genomes are in tandem arrays on chromosome 5, 3, and 4, respectively. Even the relatively few rice CBF genes also show a tandem arrangement on chromosome 9. The genomic positions and the phylogeny ( fig. 6a and supplementary fig. S32, Supplementary Material online) both support the view that tandem duplication has been the main mechanism for generating multiple copies of CBF genes in Pooideae grasses, with duplications at different stages of the Pooideae history (see fig. 6b for an illustration of the model with an emphasis for genes in Triticeaee; the genes marked with green, yellow, and red colors matching those in the CBF phylogeny are shown in supplementary fig. S31, Supplementary Material online, right). In the model ( fig. 6b), the BOP ancestor had three CBFs (represented by I, II, and III), which are maintained in rice; in addition, CBF-I was duplicated (CBF-Ix and CBF-Iy) but CBF-II and CBF-III were likely lost in Phyllostachys heterocycla (Bambusoideae). Furthermore, CBF-II experienced significant expansion in Pooideae ( fig. 6 and supplementary fig.  S31, Supplementary Material online), with some lineagespecific duplication of CBF-III (supplementary fig. S31, Supplementary Material online), resulting in at least six copies (including four CBF-II copies) in the ancestor of Brachypodieae and the core Pooideae. Further duplications of CBF-IIc in Brachypodieae increased CBF-IIc to six copies, with maintenance of CBF-IIe and CBF-III and loss of CBF-I, CBF-IIab, and CBF-IId. In the ancestor of the core Pooideae, CBF-I and CBF-II copy numbers increased to three and five copies, respectively. The expansion of CBF genes continued in Triticeae, ultimately resulting in 18 copies in the Triticeae ancestor and $20 copies in each of the three wheat subgenomes and related species, such as T. aestivum, T. dicoccoides, Aegilops tauschii, and Thinopyrum elongatum ( fig. 6b and  supplementary fig. S31, Supplementary Material online).
Following duplication, functionally redundant genes are often lost (Lynch and Conery 2000). The analyses here revealed that CBF duplicates have experienced lineage-Phylotranscriptomics and Phylogeny of Pooideae . doi:10.1093/molbev/msac026 MBE specific losses. For example, barley likely lost six copies (boxes with X in fig. 6b) from the ancestor of Triticeae, whereas wheat (including three subgenomes) and its relatives show different patterns of gene loss (gaps in gene arrays corresponding to missing numbers) across species/chromosomes, suggesting differential evolution of CBF in Pooideae especially within Triticeae. The complex CBF evolutionary pattern is also reflected by multiple gene rearrangements in Triticeae (blue lines in fig. 6b). Altogether, CBF genes in temperate grasses have gone through stepwise expansion along with different levels of gene loss and rearrangement, possibly implying that an integrated regulatory network in response to chilling has evolved progressively during Pooideae history from a more limited ancestral regulatory machinery in the BOP ancestor.

Discussion
Here, we used over 1,200 nuclear genes from 157 transcriptomes to reconstruct a highly resolved Pooideae phylogeny, with strongly supported monophyly for 14 of 15 tribes and all subtribes with two or more sampled taxa, but polyphyly of Diarrheneae. Furthermore, using the well-resolved phylogenetic relationships among Pooideae lineages, we estimated their origins and divergence times, reconstructed the evolutionary history of several traits, investigated diversification dynamics, detected numerous duplications, and analyzed molecular evolution of FUL and CBF homologs, providing clues to the success of Pooideae in adaptation to temperate environments.

Pooideae Nuclear Phylogeny Reveals Deep Cytonuclear Discordances and Suggests Possible Hybridization
In previous plastid phylogenies, the relationships were uncertain for several early divergent Pooideae tribes (e.g., Duthieeae, Phaenospermateae, Stipeae, and Meliceae) and for Poeae subtribes (Schneider et al. 2009(Schneider et al. , 2011GPWGII 2012;Soreng et al. 2017). The well-resolved phylogeny here using many more nuclear gene than possible with plastid genes reveals the order for the early divergent tribes ( fig. 1, supplementary figs. S1 and S5-S9, Supplementary Material online) and groups the Poeae subtribes into two new clades: Poeae nuclear group 1 (PNG1) and Poeae nuclear group 2 (PNG2) ( fig. 2 and supplementary figs. S2 and S5-S9, Supplementary Material online) instead of the two recognized from plastid evidence previously. PNG1 includes all eight subtribes of PCG1 plus seven subtribes of PCG2, whereas PNG2 contains the remaining nine subtribes of PCG2. The discordances between nuclear and plastid phylogenies imply a complex history across Pooideae, with possible reticulation events and/or incomplete lineage sorting, as described/discussed below and in the supplementary note, Supplementary Material online. Diarrheneae has five species in two genera with distribution in eastern Asia (Neomolinia, three species) and North America (Diarrhena, two species) ; this classification was further supported by morphological and molecular evidence (Tzvelev 1989;Schneider et al. 2011). In a study of Hochbach et al. (2015), phylogenetic analyses of 53 Pooideae (including both Diarrhena species and one Neomolinia species) with one plastid and four nuclear sequences showed that the monophyly of the Diarrheneae tribe was supported by the sisterhood of the two genera with 63-100% BS from two nuclear sequences (PhyB, Topo6 exon 17-19), whereas the results using the plastid and two other nuclear sequences (Acc1, Topo6 exon 8-13) failed to group Diarrhena and Neomolinia together (Hochbach et al. 2015). In this study, Diarrhena is maximally supported as sister to Stipeae þ Ampelodesmeae, whereas Neomolinia is sister Brachypodieae þ core Pooideae. We further examined all 1,234 nuclear gene trees and found the polyphyly of Diarrheneae in 71% of trees. Therefore, a taxonomic revision of Diarrheneae would be plausible. The discordance between nuclear and plastid possibly implies that the ancestor of Diarrhena might have involved an ancient hybridization between Neomolinia and the MRCA of the supertribe Stipodae; this hypothesis is consistent with the observation that a large fraction of nuclear genes supports a close relationship to Stipodae, whereas in plastid phylogenies and those of a small portion of nuclear genes Diarrhena is close to Neomolinia. Another possible explanation for the discordance between some of the gene phylogenies and the species topology is incomplete lineage sorting. In addition, the placement here of Phaenospermateae as sister to the supertribe Melicodae (Meliceae þ Brylkinieae) with 100% BS in all five trees is in agreement with previous analyses using nuclear sequences (79-91% BS) (Hochbach et al. 2015). However, in plastid studies Phaenospermateae was consistently grouped with Duthieeae with 75-95% BS (Schneider et al. 2009;Hochbach et al. 2015).
Previous phylogenetic analyses of subtribes within Poeae using plastid sequences often resulted in two clades, PCG1 and PCG2, but with sparse sampling at the subtribe level. On the other hand, analyses using nuclear sequences revealed that some subtribes of PCG2 (usually Airinae, Holcinae, Sesleriinae, Scolochoinae) were nested in PCG1 (Schneider et al. 2009;Tkach et al. 2020), but with uncertain relationships. Here, our phylogenetic trees consistently resolved the subtribes into two new clades (PNG1-PNG2, fig. 2 and supplementary figs. S2 and S5-S9, Supplementary Material online), with the subtribes Holcinae and Sesleriinae being nested among the subtribes of the previously defined PCG1, in partial agreement to the previous results from nuclear sequences. An examination of the concordance numbers of the 1,234 single genes trees to the alternative hypotheses ([PNG1 and PNG2] or [PCG1 and PCG2]) of the Poeae phylogeny shows higher fractions with concordance to the PNG1þPNG2 topology (supplementary fig. S33, Supplementary Material online, for details, see the supplementary note, Supplementary Material online). In addition, we evaluated the branch support by ASTRAL-pro with the function of "newick annotations" (Zhang, Scornavacca, et al. 2020), with estimation for both the quartet support (QS) and local posterior probability (LPP) for each branch. The QS represents the percentage of quartets in gene trees that agree with branch topology, whereas the LPP is computed using a transformation of the QS, with a probability value of 0.33 considered sufficient support for a branch and higher score indicating more consistent single gene trees. Our analyses of the single trees of five OGs un- These results all support the nuclear two-clade topology over the chloroplast two-clade topology.

GD via Different Mechanisms Likely Contributed to Pooideae Adaptation to Open Habitat
Among the CGDs identified here, CGD5 (circle five, fig. 4) at the MRCA of the core Pooideae includes over 4,500 GDs and coincides with a transition from a closed habitat to an open one (figs. 3 and 7), with subsequent changes in characters ( figs. 3 and 7). Among the duplicated genes in CGD5, members of gene families involved in stress response (e.g., CBF, OsWAK11, CYP7, DCL2, MPL1) and reproductive development (e.g., OsMADS1, ENL1) were particularly retained (supplementary table S11, Supplementary Material online). Specifically, molecular phylogenetic analyses of CBF homologs ( fig. 6 and supplementary fig. S31, Supplementary Material online) detected three duplications at this node. Innovations in the architecture of spikelet, including dorsal awns (at the MRCA of Poeae), hairy ovaries (at the MRCA of Triticodae), Phylotranscriptomics and Phylogeny of Pooideae . doi:10.1093/molbev/msac026 MBE and inflorescence (at the MRCA of Triticeae) (figs. 3 and 7; supplementary figs. S19-S21, Supplementary Material online) occurred after CGD5. The awn increases the chance that the spikelet is caught on the fur of passing animals (Elbaum et al. 2007), hence facilitating long-distance seed dispersal; the awn can also enhance the penetration of spikelet/seed into the soil (Elbaum et al. 2007), potentially contributing to drought tolerance (Peleg et al. 2010). The transition of the awn location from sinus or apical to dorsal might have led to higher tolerance to environmental stresses associated with open habitat, such as drought, and also promoted seed dispersal. Another change involves the ovary indumentum, from glabrous to hairy. The innovation with hairs might be associated with the resistance to water deficit and insect resistance at the MRCA of Triticodae. These changes (including molecular and phenotypic characters) could have been helpful for plant adaptation to open habitat. The correlation between GD and trait transition is statistically supported ( fig. 7 and supplementary fig. S29, Supplementary Material online), implying that GDs might have played a role in character transitions. Therefore, the increased genetic materials and associated variation due to the CGD5 might have led to molecular changes such as greater gene functions for cold-response and key shifts in morphology ( fig. 7b), and the physiological and developmental innovations then contributed to plant adaptation to new or stressful environments. In addition, it is possible that other factor(s) might have driven both the retention of duplicated genes and morphological changes. Moreover, the open habitat might have subsequently facilitated the radiation and expansion of Pooideae ( fig. 7b), such as the rapid radiation of subtribes in the core Pooideae in the late Oligocene to early Miocene ( fig. 3 and supplementary figs. S13 and S14 and table S5, Supplementary Material online). This scenario explains, at least in part, how the Pooideae "migrated" out of the canopy environment of forests to open habitats and became a highly successful lineage in the temperate flora.
Several studies suggested that GDs, especially those of WGD, have an impact on diversification (Tank et al. 2015;Landis et al. 2018). Here, we used two different methods (MuSSE and the test proposed by Tank et al. 2015) to test this hypothesis. The results support the idea that the lineages that have experienced more duplications or retained more gene duplicates often exhibit higher diversification rates (supplementary fig. S28, Supplementary Material online). Other alternative hypotheses are still possible; for example, some lineages might have been through a period of small population sizes, possibly due to relatively high diversification rates, and experienced fixation of some gene duplicates that are not necessarily functionally related to diversification. In addition, the correlation between GDs and diversification rate upshifts is statistically supported ( fig. 7b and supplementary fig. S29, Supplementary Material online). For instance, over 4,500 GDs of CGD5 that is linked to the large clade of the core Pooideae (Poeae þ Triticodae) with over 80% of Pooideae species. Not surprisingly, we detected a significant upshift of diversification at the same node/stem (figs. 4 and 7). In addition, we found evidence for WGDs associated with several genera (figs. 3 and 4) including the large genera Poa ($550 spp.) and Agrostis ($224 spp.), implying that WGD might have also contributed to species richness of groups with relatively recent origins. From these results, we propose that GDs (including CGD and WGD), at least in part, have had a positive impact on functional innovation and diversification ( fig. 7b), thereby enhancing the adaptation of Pooideae, although other hypotheses are also possible.

Climate Change Might Have Been a Driver for the Adaptive Evolution in Pooideae
Climate has had profound impacts on the origin, distribution, diversity, and even extinction of organisms (Jacobs et al. 1999;Hooker et al. 2004;Svenning et al. 2015). Since the Cretaceous, the Earth has experienced several climate upheavals including the hot periods of the Paleocene-Eocene Thermal Maximum (PETM), the Early Eocene Climatic Optimum (EECO), and the Mid-Eocene Climatic Optimum (MECO), as well as the abrupt cooling during the E-O transition (Zachos et al. 2008); especially, the E-O transition was the beginning of the last Cenozoic ice age (Kvasov and Verbitsky 1981;Goldner et al. 2014) and was suggested to have greatly affected the origin and expansion of temperate biomes (Favre et al. 2016;Meseguer et al. 2018).
Pooideae have been a prominent part of temperate ecosystems and their diversification was suggested to be driven by paleoclimate under global cooling (Schubert et al. 2019). Our results suggest that the cooling climates have been correlated with GD, diversification, and radiation across Pooideae. In addition, the origin of Pooideae and those of early divergent lineages were in closed habitats (supplementary fig. S17, Supplementary Material online) from the late Cretaceous to the late Eocene ($68-35 Ma), with relatively warm climates including the PETM and EECO and a long subsequent cooling period (figs. 3 and 7). After the EECO, there was an overall cooling trend except the abrupt warming at the MECO, with further sudden reduction of the global temperature at the E-O transition (Zachos et al. 2008). The global cooling led to the contraction of tropical forests to low latitudes (Prothero 1994;Donoghue and Edwards 2014), resulting in many new/open habitats. The changed climate and newly available habitats presented both challenges and opportunities for Pooideae plants; indeed, Pooideae experienced successive divergences of early divergent lineages (tribes) (figs. 3 and 7) and an upshift of diversification rate (figs. 3 and 7) during this cooling period, suggesting that the changing climate might have been a contributing factor to Pooideae divergence. This is also close to the period when GDs (CGD5) and the transition to open habitat occurred in the ancestor of the core Pooideae (figs. 3, 4, and 7). The CGD5 duplicates included CBF genes (figs. 6 and 7; supplementary fig. S31 and table S11, Supplementary Material online) and the increased CBF copies possibly enhanced the ability of Pooideae species to tolerate freezing during global cooling especially the E-O transition in open habitats. At the same time, some genes associated with adaptation to seasonal cooling might have already existed in the ancestor of Pooideae, as is the case of the vernalization-related AP1/ FUL-like genes, which were duplicated due to WGD shared by all grasses and have been retained in the early divergent tribes ( fig. 5 and supplementary fig. S30, Supplementary Material online). Another example is the CBF copies that were already present at the MRCA of Brachypodieae and the core Pooideae ( fig. 6 and supplementary fig. S31, Supplementary Material online). Hence, early Pooideae might already have the potential to adapt to the temperate regions before the E-O transition, which in turn could explain the survival of early Pooideae during gradual temperature decline from the late Eocene to the E-O transition.
In the late Oligocene to early Miocene, the frequent continental plate movement and orogeny continued to create new habitats, and also resulted in seasonal climates and increased aridity. This period also witnessed the radiation of most Poeae subtribes (figs. 3 and 7; supplementary figs. S13-S15 and table S5, Supplementary Material online) and divergence of other species-rich tribes (e.g., Triticeae and Stipeae). After the Mid Miocene Climatic Optimum (MMCO), there was another rapid cooling period to a lower temperature than ever before; this was also a period of multiple divergence events for many Pooideae genera (figs. 3 and 7; supplementary figs. S13-S15 and table S5, Supplementary Material online). Furthermore, the WGD events within genera (WGD2-WGD5, fig. 4) were estimated to have occurred during this aridity/cooling period, suggesting that both the WGD and climate might have accelerated speciation and radiation in Pooideae. In addition, recent dramatic cooling since $7 Ma Phylotranscriptomics and Phylogeny of Pooideae . doi:10.1093/molbev/msac026 MBE corresponded to the appearance of species adapted to cold habitats in multiple genera ( fig. 7 and supplementary figs. S13-S15 and table S5, Supplementary Material online), suggesting that cooling/cold climates might again have been an important external factor for species evolution and dispersal. Notably, the dating results can be impacted by the use of different gene sequences, phylogenies, estimation methods, fossils, and their calibrations; therefore, we need to be cautious about the ages of divergences of Pooideae lineages and the related hypothesis for links to global temperature changes. Even though the dating results using treePL here were quite consistent among the analyses here, it is possible that substantially different ages for the Pooideae lineages might be obtained using other methods and/or gene sequences, leading to alternative explanation for Pooideae cold adaptation.
The coincidence of expansion of Pooideae with cooling climate has also been suggested in previous studies (Jacobs et al. 1999;Strömberg 2011;Schubert et al. 2019). Our results here further demonstrate that in the history of Pooideae, environmental changes/stresses (such as climate change or a new type of habitat) might magnify the advantages provided by gene/genome duplications, which increased genetic materials and functional divergence and supported innovation of key characters as suggested previously (Soltis and Soltis 2016;Van de Peer et al. 2017). The combination of internal genetic events (gene/genome duplication) and external environmental stimuli (climate change) resulted in the radiation of Pooideae into multiple successful lineages with wide expansion to new territories in temperate regions.

Candidate Ortholog Identification and Gene Set Filtering
Orthologous groups were generated from gene sequences of ten representative Poaceae species (see the supplementary note, Supplementary Material online) by OrthoMCL v1.4 (Li et al. 2003). The 1,234 OGs with low-copy nuclear genes were retained as seeds and their sequences were used to search for homologous sequences in all species by HaMStR v13.2 (Ebersberger et al. 2009) with a cutoff e-value of 10 À20 . The resulting nucleotide sequences were aligned with MUSCLE v3.8 (Edgar 2004)

Phylogenetic Analysis
We used the coalescent method to reconstruct the phylogeny of Pooideae. For the low-copy genes, the phylogeny of each of the 1,234 OGs was reconstructed by RAxML v7.2 (Stamatakis 2006) with 100 replicates under GTRGAMMA model. Then, ASTRAL v5.6 (Mirarab et al. 2014) was used to reconstruct Pooideae phylogeny for each of the five OG sets with the 100 replicates from RAxML to obtain the BS values of all nodes. OGs with multiple copies in single species (from the WGD analysis) were further filtered into 802, 480, and 181 OGs (see the supplementary note, Supplementary Material online), which were used to reconstruct species trees by using ASTRAL-Pro (Zhang, Scornavacca, et al. 2020). The phylogenetic reconstruction method for wheat relatives is described in the supplementary note, Supplementary Material online.

Statistic Tests of Topological Reliability
To investigate the possibility of different topologies of Poeae, we further used ASTRAL-Pro (Zhang, Scornavacca, et al. 2020) to annotate our three-clade Pooideae phylogeny by Zhang et al. . doi:10.1093/molbev/msac026 MBE evaluating its supporting values of LPP (-t 4) and QS (-t 8) for different topologies (the main, the first, and the second alternatives). Besides, the concordance of 1,234 single gene trees onto the different topologies of the Poeae were counted by PhyParts (Smith et al. 2015).

Dating Divergence Times with Fossil Calibrations
For molecular clock estimates, the input species tree with branch lengths was generated by RAxML v7.2 using the concatenation of 373 OGs from 190 species (supplementary table S1, Supplementary Material online) while fixing the topology as shown in supplementary figure S13, Supplementary Material online. We used a penalized likelihood method in treePL (Smith and O'Meara 2012) to estimate the divergence times based on the ML tree as a phylogenetic reference. Three different calibration strategies were employed (see the supplementary note, Supplementary Material online). We first ran a preliminary analysis with the prime option to determine the best optimization parameters (opt, optad, and optcvad). Then, we performed cross-validation to determine the best smoothing value. Lastly, we performed the formal analysis with the parameters above with the fossil constraints to obtain the time tree of Pooideae. For the sensitivity test of the calibration fossils, we removed each of 15 calibrations one at a time and estimated the age with the same procedure by treePL. The average age of each node from all of the sensitivity tests was compared with above results from the three calibrations (supplementary table S5

Ancestral State Reconstruction of Characters
We selected six characters to reconstruct the ancestral states of Pooideae. Information of character states was retrieved from The Grass Genera of the World, Flora of China, and Scientific Database of China Plant Species (for more details, see the supplementary note, Supplementary Material online). The states of characters are listed in supplementary table S7, Supplementary Material online. We used the rayDISC function in the corHMM package (Beaulieu et al. 2013) in R for ancestral character reconstruction based on the ML method. Ancestral states at internal nodes were estimated by marginal probabilities. We performed the analysis with nine different models and root.p arguments, and selected the best-fitting results with the lowest AIC(c) scores and the highest weight (supplementary fig. S22, Supplementary Material online) as suggested by the author of the program.

Diversification Analysis
BAMM (Bayesian analysis of macroevolutionary mixtures) ) and MEDUSA (Modeling Evolutionary Diversification Using Stepwise AIC) (Alfaro et al. 2009) are two widely used programs to detect variation in rates across phylogeny, even though the former has received criticism (Moore et al. 2016). Here, we used the two programs to estimate diversification dynamics based on the time tree in this study with sampling fraction information (supplementary table S8, Supplementary Material online). We calculated the sampling fraction at the tribe level based on species information of Soreng et al. (2017). The detailed description of the implementation BAMM and MEDUSA is described in the supplementary note, Supplementary Material online.
Gene Tree Mapping, Synonymous Substitution Rate (Ks) Estimation, and Synteny Analysis To identify GD clusters as possible evidence for WGDs in Pooideae, a total of 164 data sets including eight outgroups were used to assign genes to gene families essentially following previous methods ) (see the supplementary note, Supplementary Material online, for details). Finally, 45,722 of aligned homolog group were used to reconstruct a final ML gene tree by RAxML v7.2 (Stamatakis 2006) with GTRCAT model. We performed tree reconciliation analysis by comparing gene trees with a reference species tree (supplementary fig. S23, Supplementary Material online) as described in previous studies Xiang et al. 2017;Ren et al. 2018;Leebens-Mack et al. 2019). The GD was counted at each of the MRCA node by Tree2GD (https://github.com/Dee-chen/Tree2gd; last accessed July 25, 2021) with !50% BS of the divergence of two subclades. The synonymous substitution per site (Ks) for paralogs in a specific taxon was calculated as described previously (Leebens-Mack et al. 2019;) (see the supplementary note, Supplementary Material online). Moreover, gene collinearity in assembled genome was analyzed by MCScanX , then the intersection between syntenic gene pairs and duplicates derived from the phylogenetic method of GD cluster is used as evidence of genome doubling.

Annotation of Expanded Gene Families and GO Enrichment Analysis
Among genes that experienced GD during Pooideae evolution, we selected the genes (gene families when considering homologs from multiple species) with species coverage !20 species and !5 species in both of the duplicated subclades, to represent genes that have expanded in a wide range of species. Then, the gene sequence(s) from any of the three Pooideae genomes (T. aestivum, Hordeum vulgare, or B. distachyon) were used to identify homologs in Oryza sativa or Arabidopsis thailiana, with their annotated functions listed in supplementary table S11, Supplementary Material online. We further carried out GO enrichment analysis using AnnotationHub (Morgan 2019) and clusterProfiler (Yu et al. 2012) packages of R to find the enriched functional categories for duplicates of all candidate WGDs and CGD5.

Molecular Evolution of AP1/FUL and CBF Homologs
To reconstruct the evolutionary history of AP1/FUL homologs, we selected a total of 46 species (supplementary table S13, Supplementary Material online). We retrieved candidate AP1/FUL homologs of all representatives by using HMMER Phylotranscriptomics and Phylogeny of Pooideae . doi:10.1093/molbev/msac026 MBE v3.1 (Mistry et al. 2013) with the entry domain (PF01486.17) and default settings. To identify CBF homologs, with are members of a subclade of the AP2/ERF superfamily containing the AP2 domain, we performed a BLASTP search (e-value cut off ¼ 10 À5 ) to favor genes that are more closely related to CBF. The treatments of the retrieved protein sequences are described in the supplementary note, Supplementary Material online. In addition, we used MCScanX ) and location of gene on the chromosome to confirm the type of duplication.

MuSSE Analysis for GD-Related Diversification Rates
To investigate potential association between the GDs (CGD and WGD) and diversification rate, we performed MuSSE analysis (Multiple State Speciation and Extinction) (Maddison et al. 2007) to investigate the effect of times of duplication on diversification rate in the resolved phylogeny (FitzJohn et al. 2009). The R package of diversitree v0.9.11 (FitzJohn 2012) was used for the analysis. The numbers of times each species has experienced GDs detected in this study are used as traits. Transition of character states were constrained to change stepwise as 1<>2<>3<>4-8. We selected the best-fitting model from four different models (i.e., free lambda and mu, fixed lambda, fixed mu, and fixed lambda and mu for each state) by the lowest AIC value and the highest AIC weight (supplementary table S12, Supplementary Material online). Then we ran the MCMC chain for 5,000 generations and sampled every 100 generations with an exponential prior as 1/(2r), where r was the lambda of the species with the original state.

Correlations between Diversification Rate Upshifts, Trait Transitions, and GDs (CGD and WGD)
A statistical test for the correlation of diversification rate shifts, trait transition, and GDs (CGD and WGD) was performed following the method of Tank et al. (2015) (also see Landis et al. [2018]) as the coincidence between two of the three features (diversification rate shifts vs. trait transition, trait transition vs. GDs, and GDs vs. diversification rate shifts). For example, for association between GDs and diversification rate shifts, the number of coincidence of observed data (e.g., GDs and diversification rate upshifts observed in this study) were compared with two alternative null hypotheses: 1) coincidences between the observed upshifts and randomly distributed GDs (r.wgd); 2) coincidences between the observed GDs and randomly distributed upshifts of diversification rates (r.shift). Each test with 1,000 random distributions was repeated 100 times to produce the result of the distribution of P values.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.