Both mechanism and age of duplications contribute to biased gene retention patterns in plants

Origin and turnover rate of paralogs

For each of the 25 study species, we calculated Ks time divergence between pairs of paralogs and used a synteny-based approach to categorize members of all gene families as derived from WGD or tandem duplications. Duplicate pairs whose origins were uncertain based on available data were classified as “undefined”. Across the 25 target genomes, the majority of paralogs detected had Ks ≤ 2 (Table 1) including 79% of paralogs in A. thaliana, 86% in Glycine max, and 92% in Malus domestica. Paralogs with Ks > 2 were excluded from our analyses due to concerns that Ks saturation could impair reliable inferences [31]. Most species displayed clear prominent peaks in their Ks age histograms, which is illustrated by histograms for five species with contrasting histories of polyploidy (Fig. 1). Histograms for the remaining 20 species are depicted elsewhere (Additional file 1: Figure S1). In the K-S goodness of fit test, all histograms for all species except Carica deviated significantly (P < 0.05) from the null model of constant duplicate gene birth and death (Additional file 1: Table S1). SiZer maps identified a significantly increasing gradient in the Ks age histograms of WGD-derived paralogs of most species, which provides support for polyploid signals being well distinguished from background duplications.

Specie	Chr	Initial PCG	Duplicates	Number of duplicates by duplication type			Number of duplicates by Ks range
				WGD	Tandem	Undefined	0 < Ks ≤ 0.5	0.5 < Ks ≤ 1	1 < Ks ≤ 1.5	1.5 < Ks ≤ 2	Ks > 2
Arabidopsis lyrata	16	32670	6378	3442	1816	1120	2251	2216	966	945	228
Arabidopsis thaliana	10	33602	6194	2740	1232	2222	1657	2407	1183	947	222
Amborella trichopoda	26	26460	3322	15	998	2309	1861	427	402	632	137
Brachypodium distachyon	10	26678	3573	1025	1768	780	981	1024	835	733	175
Carica papaya	18	28072	1915	24	455	1436	603	210	402	700	126
Citrullus lanatus	22	23438	2806	385	1015	1406	691	435	751	930	227
Eucalyptus grandis	22	36449	11120	390	6424	4306	8106	925	1029	1060	240
Fragaria vesca	14	34809	3974	1021	1606	1347	1500	979	684	811	184
Glycine max	40	46509	15242	9721	2087	3434	11697	1961	790	794	185
Helianthus annuus	34	44144	9925	57	995	8873	4651	2805	1551	918	108
Lotus japonicus	12	26818	2682	184	627	1871	1159	774	415	334	51
Malus domestica	34	63515	15551	2761	1308	11482	13084	1258	683	526	107
Manihot esculenta	36	30800	7134	2530	703	3901	4915	837	716	666	110
Medicago truncatula	16	57587	5098	1083	2419	1596	2902	1262	543	391	115
Oryza sativa ssp. indica	24	48788	8349	1957	2869	3523	3361	2169	1665	1154	317
Oryza sativa ssp. japonica	24	59430	5559	1482	2173	1904	1928	1584	1233	814	183
Physcomitrella patens	54	36137	3769	306	202	3261	637	1848	883	401	99
Populus trichocarpa	36	41521	9721	5609	1988	2124	7572	738	704	707	147
Ricinus communis	20	31221	2558	155	614	1789	628	435	683	812	176
Sorghum bicolor	20	34686	4267	1048	1698	1521	1468	1061	993	745	186
Solanum lycopersicum	24	34432	7100	1234	2561	3305	3184	2287	872	757	209
Selaginella moellendorffii	16–27	22285	1885	351	608	926	1457	129	102	197	66
Theobroma cacao	20	46269	3488	722	1553	1213	1199	601	822	866	201
Vitis vinifera	38	26644	4536	528	1935	2073	1918	852	1042	724	128
Zea mays	20	39597	6336	590	1396	4350	3792	1095	813	636	153

Chr Number of Chromosomes, Initial PCG Initial number of Protein-coding gene sequences

The Ks age histograms of WGD-derived paralogs (Fig. 1a, depicted in black) were clearly distinct from those of the tandem-derived paralogs (Fig. 1b, depicted in gray). While tandem histograms exhibited a descending slope (similar to a half-parabola) for most of the species, WGD-derived paralog histogams had peaks that overlapped with peaks from histograms of all paralogs (Fig. 1a, depicted in brown). SiZer maps also confirmed the presence of peaks for WGD-derived paralogs histograms (Fig. 1d).

Because of our focus on transcription factor paralogs, their Ks age histograms are shown (Fig. 1 and Additional file 1: Figure S1; depicted in yellow) along with the Ks histograms of WGD- and tandem-derived paralogs. The SiZer maps (Fig. 1e) showed increasing gradients for transcription factor paralogs that overlapped with the slopes of WGD-derived paralogs.

Biased retention of paralogs after large- and small-scale duplications

To assess the universality of the GBH across land plants, we identified the most strongly overrepresented GO functional categories in both predicted WGD- and tandem-derived paralogs in these 25 genomes. We found that WGD- and tandem-derived paralogs did not share the top 10 most frequently overrepresented GO categories (Fig. 2a and c). While the most overrepresented categories of WGD-derived paralogs fell under macromolecular complexes (GO:0032991), internal to cell (GO:0005622), and cytoplasm (GO:0005737) functional GO categories; those of tandem-derived paralogs grouped into programmed cell death (GO:0012501), defense response (GO:0006952), and apoptotic process (GO:0006915) GO categories.

In six species, WGD-derived paralogs were not enriched for the overrepresented GO categories found in the remaining plant species. Five of them—Cariaca, Ricinus, Populus, Selaginella, and Physcomitrella—have few WGD-derived paralogs predicted by DAGchainer (Table 1), consistent with possible under-estimation or misidentification of WGD-derived paralogs in these species (see Discussion below). For another five taxa—A. thaliana, Medicago, both Oryza subspecies, and Populus—WGD-derived transcription factor paralogs were overrepresented (Fig. 2b). Surprisingly, WGD-derived transcription factor paralogs were not significantly overrepresented in Arabidopsis lyrata, which shares the same WGD events as A. thaliana, although there was a trend in the expected direction.

Unexpectedly, transcription factor activity (GO:0003700) WGD-derived paralogs were not significantly overrepresented in 20 plant species, ten of which exhibit evidence of recent WGDs in their evolutionary history, with a significantly increasing gradient in SiZer (Fig. 1 and Additional file 1: Figure S1) within Ks range < 1 (and consistent with previous reports—see below). Finally, results from our analyses of tandem duplications showed tandem-derived transcription factor paralogs were significantly underrepresented across the 25 focal genomes.

Biased retention toward ancient transcription factors

We analyzed the biased retention of transcription factor paralogs based on Ks time divergence as opposed to mechanism of duplication. This was accomplished by mapping known WGD events onto a phylogeny for the 25 species targeted by this study (Fig. 3, Additional file 1: Table S2).

In general, transcription factor (GO:003700) paralogs tend to be overrepresented amongst ancient (Ks > 1) duplication regardless of mechanism of duplication (Fig. 3). Eleven of the 25 focal species exhibited significant enrichment at Ks range > 1.5 (Fig. 3, pentagons), but no such retention bias at lower Ks ranges (≤1.5). When we compared transcription factor paralog enrichment at Ks > 1.0 versus < 1.0, 17 species showed significant enrichment for the older transcription factor paralogs (Fig. 3, stars). For four of these, A. thaliana, Medicago, and the two Oryza subspecies, the overrepresented transcription factors originated from WGD events (Fig. 2b). However, for the remaining 13 species, the ancient paralogs are not obviously associated with a WGD event. Although A. thaliana, Oryza sativa ssp. indica, and Solanum exhibited significant signals of polyploidy in the Ks range < 1 (Fig. 1; Additional file 1: Figure S1), their transcription factor paralogs were only significantly overrepresented in the Ks range > 1 (Fig. 3f).

In genomes of only four taxa (Carica, Malus, Manihot, and Populus) were recent transcription factor paralogs overrepresented, and only for Populus were WGD-derived transcription factor paralogs significantly overrepresented (Fig. 2b).

In addition to analyzing the retention of transcription factor paralogs, we submitted our data to enrichment analysis aiming to find additional GO categories that could have experienced biased retention patterns. A number of GO categories, including those involved in transcription, regulation, transport, and response to stimulus were frequently overrepresented among ancient paralogs (Ks > 1) and not exclusively associated to WGDs (Additional file 1: Figure S2). While three of these functional GO categories—cell periphery (GO:0071944), plasma membrane (GO:0005886), and response to abiotic stimulus (GO:0009628)—were overrepresented among WGD-derived paralogs; two categories—response to stimulus (GO:0050896) and catalytic activity (GO:0003824)—were overrepresented among tandem-derived paralogs.

Tandem paralogs have faster turnover rate

Our synteny-based approach identified pairs of WGD-derived genes similar to those that have been reported in previous studies. In A. thaliana, for example, circa 80% of the 2740 duplicate gene pairs we classified as WGD-derived are in common with the list of polyploidy-derived paralogs published by [23]. Differences among studies may be due to new gene annotation tools that recently became available. In some instances, the number of paralogs predicted as having their origin in WGD events can be underestimated due to widespread genomic changes (e.g., gene loss and/or chromosomal rearrangements) after polyploidization events [19]. Such processes are particularly problematic for ancient polyploidization events, which may explain the low number of WGD paralogs we predicted in the basal plants, Amborella and Physcomitrella, as well as for Lotus, Carica, and Ricinus (Table 1). On the other hand, our approach indicates the presence of a small number WGD-derived paralogs in Selaginella, which is not thought to have a WGD in its evolutionary history (Table 1). This result could be evidence for Selaginella as ancient polyploidy. Alternatively, it suggests that Selaginella has had an ancient large segmental duplication or some fraction of the identified WGD derived paralogs are false positives. However, we selected WGD pairs using a syteny based approach, which is the most conservative method presently available.

Tandem paralogs were similarly identified based on the genomic coordinates of genes. In Eucalyptus, 32% of its 36,449 protein-coding genes originated via tandem events, which is the largest proportion of tandem-derived paralogs amongst the 25 plant species we investigated. Physcomitrella exhibited the smallest proportion (~1%) of tandem-derived paralogs. These findings are very similar to those previously reported for Eucalyptus [32] and Physcomitrella [33], respectively.

We identified peaks in the Ks age histograms; based on SiZer maps, these peaks likely result from WGDs (Fig. 1 and Additional file 1: Figure S1). Previous studies have also identified these WGD events using data that span across several families [34], or from a given plant species [9, 10, 32]. In the Ks histogram of A. thaliana, for example, there were two prominent peaks (Fig. 1), which coincided with the α and β polyploid events reported by early investigations [34–36]. In our analysis, the tail of the most recent duplication masked the second peak; thus, a single, significantly increasing slope was identified by SiZer. In A. lyrata, SiZer identified two significant peaks as expected given the recent history of polyploidy in Arabidopsis [36].

Differences in the Ks age histograms from WGD- and tandem-derived paralogs indicates that the turnover rate of tandem paralogs is faster than that of WGD paralogs, as previously suggested by others [7, 14, 29, 33]. The pattern we uncovered suggests lower turnover rates of transcription factor paralogs than those observed for tandem paralogs. Furthermore, it appears that the origin and biased retention of transcription factor paralogs are not restricted to large-scale duplication events.

Patterns of transcription factor retention following WGDs

Consistent with the expectations of the GBH, WGD- and tandem-derived paralogs did not share the top 10 most frequently overrepresented GO categories. Six species—Malus, Cariaca, Ricinus, Populus, Selaginella, and Physcomitrella—were exceptions and did not share the most frequent GO categories, which is consistent with the possible under-estimation or misidentification of WGD-derived paralogs in these species. In Malus, for example, the GO categories that were overrepresented include: plasma membrane (GO:0005886), response to abiotic stimulus (GO:0009628), response to biotic stimulus (GO:0009607), and response to endogenous stimulus (GO:0009719). Analyses of an EST library of Malus domestica also found that these categories were overrepresented [37]. Consistent with the GBH, we did not found tandem-derived transcription factor paralogs overrepresented in any of 25 focal genomes.

Other findings were inconsistent with the predictions of the GBH. In plants, the genome of A. thaliana has been frequently used to support dosage-constraints of transcription factors [23, 24]. Unexpectedly, our findings reveal transcription factor activity (GO:0003700) WGD-derived paralogs to be significantly overrepresented in only five plant taxa—A. thaliana, Medicago, the two Oryza subspecies, and Populus. Ten of the 20 remaining study species exhibited evidence for recent WGDs. Other studies have also reported a downward bias in the retention of transcription factor paralogs following WGD events. In Compositae paleologs, for example, it has been observed that genes involved with structural components or cellular organization were significantly overrepresented; whereas transcription factors were significantly underrepresented [10]. These authors argued that patterns of intrinsic selection on different gene categories may vary across higher taxonomic categories. The fate of paralogs originated by either WGD or small-scale events would depend on intrinsic properties, such as gene function and the environment in which the new polyploid was born [21].

Age of duplications contribute to biased gene retention

Regardless the mechanism of duplication, we showed that ancient paralogs of transcription factors were preferentially retained over paralogs of more recent origin. In agreement to our findings, a previous study in A. thaliana reported that genes involved in transcriptional regulation showed greater retention after the later (β) genome duplication than after the youngest (α) duplication [24]. Likewise, transcription factors not directly associated with WGDs were overrepresented among genes of ancient origin in A. thaliana [13]. Again, our results indicate that out of 25 plant species with very different histories of polyploidy, such as A. thaliana which has two recent WGD events [35] and Vitis which has no known recent WGD events [9], 17 share this pattern of biased retention of ancient transcription factor paralogs. Although transcription factor paralogs with recent origin were over-represented in four species (Carica, Malus, Manihot, and Populus), we could only clearly determine that those of Populus were WGD-derived paralogs. The over-representation of young (Ks < 0.5) transcription factor paralogs in Carica is intriguing, given that no WGD events likely took place in its recent evolutionary history [38] and that DAGchainer only predicted tandem-derived transcription factor paralogs for Carica within the Ks range ≤ 1.0 (Additional file 1: Table S2). Given that Carica lacks recent WGD events [38] and we did not identified transcription factors paralogs originated from WGD events within Ks < 1, the many transcription factor paralogs of Carica appear to derive from small-scale duplications within its genome.

Our findings differ from a recent study of core gene families in 37 angiosperm genomes [13], which reported remarkable consistency in the rate at which genes return to a single copy state, as well as in the gene families that are retained as multi-copy. The findings were related to differences in gene function and the authors concluded that similar selection pressures within and between lineages are largely responsible for the non-random patterns observed, at least for core genes [13]. The apparent differences between the two studies derive partly from the fact that core gene families represent a fairly small fraction (13%) of all gene families and that single copy genes in were included their analysis, which drive many of the reported patterns. In contrast, we restricted our analyses to duplicate genes.

Data collection and selection of paralogs

Full genome annotations, protein-gene codes, DNA sequences, gene families, and Gene Ontology (GO) annotations from the 25 focal species were retrieved from PLAZA 2.5 and 3.0 Dicots [39], with the exception of sunflower (Helianthus annuus), as detailed in Additional file 1: Table S3. Protein-gene code files with alternative transcripts removed were used to identify paralogous gene pairs using BLASTp all-against-all, with an e-value cutoff of e⁻²⁰, with a minimum 50% identity, alignment length > 300 bp, number of mismatches < 550, and number of gap opens < 30. Self hits were removed and only paralogous gene pairs with both copies belonging to the same gene family were maintained for further analysis. For the selection of paralog pairs for H. annuus, CDS sequences and BLASTn all-against-all were used based on HA412.v1.1 version of the genome (http://www.sunflowergenome.org/).

Determining paralog duplication mechanism

The DAGchainer software package [40] was used to predict the mechanism of by which paralogs originated based on their genomic coordinates. WGD-derived paralog pairs were predicted by running DAGchainer to find syntenic/collinear regions among chromosomes, in the same species, using default parameters and ignoring tandem duplication alignments (-s and -I options). Tandem-derived paralogs were predicted by using the accessory segmental duplication tool, also made available by DAGchainer, to find collinear sets of homologous genes, with the ‘max intervening genes value’ set to 10. All the other paralog pairs, not predicted as WGD or tandem-derived, were marked as undefined (UD), as these paralogs may have been originated by either large- or small-scale duplications.

Age of duplication events

We calculated relative divergence times for each paralog pair in terms of synonymous substitutions per synonymous site (Ks). First, we aligned the nucleotide sequences of gene pairs using TranslatorX [41], based on protein alignments performed by MUSCLE v3.8.31 [42]. Divergence times (Ks) were calculated with the yn00 software from the PAML v4.1 package [43]. This method assumes the F3x4 codon frequency model and accounts for transition/transversion rate bias and codon usage bias, which is an approximation of the maximum likelihood method recommended for pairwise comparisons in the manual of PAML. Because of issues associated with Ks saturation and stochasticity [31], only paralogs with Ks ≤ 2 and Standard Error (SE) < 0.5 were used in further analyses.

Custom python scripts were used to parse the BLAST all-against-all output in order to identify the closest paralog gene pairs. First, self hits were removed. Then, paralogs were organized into a single gene list and then used to select the corresponding paralog pair(s) for each of these genes based on the following three rules: (I) if a single gene was predicted as WGD-derived by DAGchainer, keep the duplicate pair with the lowest Ks value, while still allowing pairing with tandem-derived or undefined genes; (II) if not predicted as WGD-derived, but predicted as tandem-derived, keep the gene pair with the lowest Ks value; and (III) if the single gene was not predicted as WGD- or tandem-derived, keep the undefined paralog pair with the lowest Ks value.

GO annotation and over-representation analysis

Functional Gene ontology GO terms (categories) were determined for each gene and paralog pair and then evaluated for enrichment by the ErmineJ v3.0.2 software [44]. All the three GO domains (Biological Process, Molecular Function and Cellular Component) were included in the Over-Representation Analysis (ORA), with a minimum gene set size equal to 10 and the Best Scoring gene replicate treatment. Eight different groups of paralogs were analyzed: WGD-derived, tandem-derived, and paralogs representing the following Ks ranges: (A) 0 < Ks ≤ 0.5, (B) 0.5 < Ks ≤ 1, (C) 1 < Ks ≤ 1.5, (D) 1.5 < Ks ≤ 2, (E) 0 < Ks ≤ 1 and (F) 1 < Ks ≤ 2. The GO categories were considered overrepresented if Corrected P < 0.05, as calculated by the ErmineJ software.

Statistics

SiZer maps: identifying significant peaks in Ks histograms

Significant peaks in the Ks histograms were found by SiZer [46] implemented on R software, with the following command line: SiZer.1 < − SiZer(x, y, h = c(.05,5), degree = 1, derv = 1). A SiZer map is a way of examining when the p-th derivative of a scatterplot-smoother is significantly negative, possibly zero or significantly positive across a range of smoothing bandwidths. In a SiZer map, blue indicates a significantly increasing gradient, red is a significantly decreasing gradient, purple is a non-significant gradient and gray indicates that data are too sparse for reliable estimation.