We inferred simultaneously duplicated genes produced by the ceWGT in 28 representative plants selected from 28 plant families and constructed a phylogenetic tree by integrating prior information on species relationships ( Figure 1 ). First, using the most updated gene collinearity inference toolkit WGDI (Sun et al., 2022), we extracted the collinear gene pairs within each genome of the studied species. In these plant genomes, we detected 479–5,024 homologous/paralogous blocks containing 5,081–51,349 paralogous genes in collinear positions along the compared chromosomes ( Supplementary Table S2 ). The mean number of duplicated blocks is 1,844, and the mean number of duplicated genes is 22,466. Among these plants, Boehmeria nivea has the fewest colinear blocks (479) and collinear genes (5,081), while Hevea brasiliensis has the highest number of collinear blocks (5,024) with 51,349 collinear genes. More collinear blocks show more chromosomal rearrangement or larger-scale DNA fractionation suggesting divergent stability among considered genomes.

Second, we estimated the synonymous nucleotide substitutions (Ks) between duplicates genes. This effort confirmed that each considered plant has been affected by another polyploidization event after the ceWGT and showed that the latter polyploidization often occurred tens of millions of years after the ceWGT ( Figure 2 ). Considering difficulties to separate paralogous genes produced by more than two polyploidization events or by two events close in time, this study prioritizes species that each have experienced the ceWGT and an additional WGD that occurred much later than the former ( Figures 2 , 3 ). This strategy facilitates a precise discrimination of paralogous genes produced by the considered polyploidization events and estimation of Ks values corresponding to each event. Besides, a check of homologous gene dot plots helped separate the ceWGT-derived duplicates from those produced by additional polyploidization ( Figure 3 ). For instance, in addition to the shared ceWGT, initially revealed in grapevine, mango has undergone an additional WGD event.

Actually, the ceWGT event resulted in the production of 284–4,142 paralogous gene blocks, and 1,746–41,113 paralogous genes in these plant genomes. The most paralogous genes were found in Fagus sylvatica, the fewest in Lactuca sativa, while the average is 9,903.

Owing to long tails often present in the Ks distributions, we extracted the median Ks value of each inferred collinear gene blocks, which often provide statistically stable inference, and then characterized the Ks distributions of the medians from all collinear blocks from each genome. The obtained Ks distributions were then fitted according to the normal distribution to extract the peak value of Ks ( Supplementary Table S3 ).

In our analysis, we observed divergent Ks distribution patterns among the plants under consideration. Specifically, for the shared ceWGT event, the maximum Ks value was detected in Simmondsia chinensis, reaching up to 2.0408, whereas the minimum was found in Nyssa sinensis, as low as 1.2145. The average Ks value across the shared ceWGT events was 1.5989. In contrast, for pitaya (Hylocereus undatus), the Ks distribution exhibited a peak at 2.0142, which is approximately two-thirds larger than the peak observed in Nyssa sinensis. We divided the studied plants into two groups: one group having been affected by a single round of polyploidization (the ceWGT) and the other one by additional round of polyploidization (a WGT or a WGD), which were denoted as P1R and P2R, respectively ( Figure 4 ). The P1R group involves 12 plants, and the P2R includes 16 ones. For the P1R plants (including wild strawberry, Acer yangbiense, starfruit, cocoa tree, and grapevine), the grapevine Ks values peaked at 1.2564, much smaller than that in strawberry (Ks = 1.6221, 29.1%) and that in jujube (Ks = 1.5890, 26.5%). For species in the P2R group, an even more significant variation was observed in their Ks distribution. For instance, the ceWGT Ks value in dragon fruit peaked at 2.0142, being 51.13% higher than that in tea tree ( Figure 2 ). A divergent Ks distribution shows a rather divergent evolutionary rate difference among plants, to be further discussed below.

The above characterization of the ceWGT revealed thousands of simultaneously duplicated genes in each species, which enables us to check whether these paralogous genes have gathered divergent levels of synonymous nucleotide substitutions, to further show whether the studied plants have evolved at different paces.

To explore whether polyploidization contributes to the elevation of plant evolutionary rates, we compared the two groups to find whether there is difference between them. The coefficient of variation is 0.161 for the P1R group, and 0.149 for P2R, indicating a higher degree of Ks dispersion in the former. This reflects a greater heterogeneity in the evolutionary rates among P1R species. The P1R group has a much lower average Ks (1.557 ± 0.251) than that (1.631 ± 0.243) of the P2R group (t-test p-value = 0.001) showing an elevated average evolutionary rate of the latter. This shows that the P2R paralogous genes produced by the ceWGT have gathered 4.75% more nucleotide substitutions than those of the P1R. Grossly, combining the above findings, no matter what the plants are or where they are originated, we showed that an additional polyploidization can significantly increase the evolutionary rates of plants. Moreover, this justify that the occurrence of polyploidization elevates a species’ evolutionary pace and results in evolutionary effects in thousands of millions of years.

Divergent evolutionary paces make it difficult to perform a reasonable dating of the evolutionary events, such as ancient polyploidization or speciation of plants, under the traditional evolutionary models. These models often based their inference on the same or similar evolutionary rate(s) and other parameters. In fact, if we directly infer the occurrence time of the ceWGT, using the Ks distribution from different plants and assuming the same evolutionary rates (often a Ks rate 7.0 × 10⁻⁹ synonymous substitutions per site per year), we would have much diverged estimation in different plants, varying from 90 to 164 millions of years. Besides, they often used the same gene families, which evolve at paces affected by the plants from which they are derived, and divergent plant evolutionary paces were not well considered. Additionally, scientists often did not have knowledge whether the previously used (duplicated) genes had the same origination, let along whether they were produced simultaneously.

Here, using genome-scale data (thousands of paralogs from each plant, simultaneously duplicated by the ceWGT), we managed to perform a correction to the Ks between paralogs from the studied plants. Grapevine was used as a reference, in that it preserved much of the genome structure of the ancestral core eudicot genome, showing a stable genome sequence during more than 100 million years of evolution. The correction was made by a linear transformation of the mean Ks values of the paralogs from the other plants to that in grapevine (see Methods for details; Figure 5 ; Supplementary Table S4 ). Actually, the Ks between the grapevine paralogs has a mean value of 1.2564, which is among the smallest ones in all plants, again showing the stable nature of the grapevine genome. After performing the correction, all of the Ks distribution peaks from different species were aligned to the same location as that of grapevine’s ( Figure 5B ).

Supposing that the ceWGT having occurred approximately 115 to 130 million years ago (Yang et al., 2020b), as previously inferred, during the Paleogene period (approximately 23 to 66 million years ago), there was a significant scale of collective polyploidization events. In this period, species, such as M. indica, H. brasiliensis, C. illinoinensis, and N. sinensis, and their respective relatives, experienced WGDs.

After correcting for the evolutionary rates, we re-estimated the timing of the most recent polyploidization events in a range of plant species ( Figure 5C ; Supplementary Table S5 ). For example, the results indicate that the recent polyploidization event in lettuce (Lactuca sativa) was estimated to occurred at approximately 59 to 67 million years ago (Mya), while pitaya (Hylocereus undatus) did so around 47 to 54 Mya. In addition, mango (Mangifera indica) had a specific polyploidization event occurring ~22 to 25 Mya.

Utilizing a stratification method based on the synonymous substitution rate (Ks), we systematically ranked the paralogous genes according to their Ks values in each of the selected species from 14 different orders. Considering the genes that fall within the top and bottom 20% of the Ks distribution, we identified the fast- and slowly evolving genes in each species.

Often, the largest proportion of genes in both groups were involved in pathways related to “Genetic information processing” or “Protein families of genetic information processing.” An only exception was found with M. indica, in which the largest proportion of the fast-evolving genes were involved in pathways relating to “Protein families of signaling and cellular processes.” Interestingly, we found that the fast-evolving genes in all species and the slowly evolving ones in 13 species were not involved in the pathways relating to “Metabolism of cofactors and vitamins” (the exceptional one: I. triloba) or “Nucleotide metabolism” (the exceptional one: C. canephora). Similarly, the fast-evolving genes were not involved in pathways relating to “Metabolism of other amino acids” in all selected species, while the slowly evolving genes in four species were involved to an unneglectable percentage (4.44%–7.14%) ( Figure 6 ).

Taking grapevine as an example ( Figure 7 ), according to the annotation analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG), a Chi-square test showed that genes evolving at different rates were significantly divergently related to biological pathways (χ² = 43.976, p-value = 0.0001) ( Supplementary Table S6 ). Actually, for the fast-evolving group in grapevines, genes involved in the “metabolism” category accounted for ~40% of the total, the most of which were divided into two groups: “carbohydrate metabolism” and “lipid metabolism.” Approximately 36% of the fast-evolving were related to “genetic information processing” and 9% to “environmental information processing.” As to the slowly evolving group, 50% of them belonged to the “metabolism” category, 10% more than those in the fast-evolving group. Among the “metabolism” category, the slowly evolving genes were mainly related to “biosynthesis of other secondary metabolites” (16%) and “carbohydrate metabolism” (9%), while genes related to “genetic information processing,” “cellular processes,” and “organismal systems” accounted for 33%, 10%, and 7%, respectively.

We collected 28 high-quality, chromosome-level core eudicot plant genomes (mainly from NCBI and PHYTOZOME). The data materials mainly include genome annotation files (General feature format, GFF), gene translated protein files (Peptide, PEP), and coding sequence files (CDS) ( Supplementary Table S1 ). Python scripts (https://github.com/SunPengChuan/wgdi) were used to process data format to facilitate subsequent research.

The first step is to perform multiple sequence alignment. According to the sequence alignment tool BLAST (Altschul et al., 1990), select the -blastp module to perform homology alignment within and between genomes of the selected species’ genomic protein sequences. The E-value threshold of the output result is set to 1e⁻⁵ to accommodate the duplicated genes produced by paleopolyploidization ∼10 Mya, and the output file format (-outfmt) is set to 6 during the specific operation.

Based on the previously obtained homologous sequence alignment result file (Blast file), combined with the genome annotation file (Gff file) and chromosome length file (Lens file), we used the -d module implemented in WGDI (Sun et al., 2022) to draw homologous gene dot plots. The dots of different colors (red, blue, and gray) in the dot plots represent the level of similarity of their gene pairs. According to the dot plots of homologous genes, homologous collinearity was inferred within each genome and between genomes. Next, we used the -icl module to perform collinearity analysis and obtained the collinearity regions, described by scores, statistical significance, collinear gene numbers, etc.

According to the collinear genes obtained previously, the Ks values were calculated by combining the cds and pep files. Here, the -Ks module in WGDI was used. This module used Muscle software (Rc, 2021) to perform protein binding based on the protein sequence, used pal2pal.pl to convert the protein binding into codon binding based on the CDS sequence, and finally calculated Ks using yn00 from PAML (Yang, 2007).

We performed the Ks correction among species by aligning the ceWGD peaks of Ks values between duplicated genes in each genome ( Figure 5A ) to the peak location in grapevine.

Using the aforementioned method for calculating Ks values, we calculated the correction coefficients for duplicated genes within individual genomes. We used K s ' to denote the Ks values after correction, which was obtained by multiplying the original Ks values by the coefficient (Zhuang et al., 2019; Song et al., 2021b).

For any species q, its correction coefficient is defined as q, its correction coefficient is defined as C q = K s ' / K s = P V i / P q = k ; thus, K s q ' = K s ×   C q = K s ×   ( P V i / P q ) = K s × k . Among them, P q represents the peak value of the Ks distribution of species q.

For example, the correction coefficient of duplicated genes in diploid wheat C. sinensis was defined as C C a = K s ' / K s = P V i / P C a = 1.0608 ; therefore, K s ' = K s ×   C C a = K s ×   ( P V i / P C a ) = K s × 1.0608 .