we obtained the genome sequences of a. thaliana ( arabidopsis genome initiative 2000 ) , rice ( o. 
 sativa , international rice genome sequencing project 2005 ) , and soybean ( g. 
 max , schmutz et al . 2010 ) from the plant genome database ( duvick et al . 2008 ) and the plant duplication database ( lee et al . 2013 ) along with their outgroups a. 
 lyrata ( hu et al . 2011 ) , sorghum bicolor ( paterson et al . 2009 ) , and phaseolus vulgaris ( schmutz et al . 2014 ) , respectively . polymorphism data were obtained from 80 a. 
 thaliana accessions ( cao et al . 2011 ) . to identify coding single nucleotide polymorphism ( snp ) information for rice we used the rice haplotype map project database ( second generation , http://www.ncgr.ac.cn/ricehap2/index.html ) and soybean we used snp information deposited in snpdb ( sherry et al . gene expression data were obtained for the three plants species from the plexdb ( dash et al . in particular , we focused on development stages preceding gametogenesis , during gametogenesis and embryogenic developments ( table 1 and supplementary file s1 , supplementary material online ) . 2003 ) and invariant set ( is ) methods were performed with the affy bioconductor package to normalize all data sets simultaneously . scatterplots of expression between replicates showed better results for rma normalization ( data not shown ) . therefore , unless stated otherwise , expression data shown in this study are based on a normalization across experiments using rma with log2 transformation . since different laboratory conditions can affect expression patterns ( massonnet et al . 2010 ) , we controlled for these effects in the a. thaliana data ( schmid et al . 2005 ) by removing data sets that were obtained from plants with different growth conditions before rna extraction ( supplementary file s1 , supplementary material online ) . to check whether the differences in expression between experimental conditions were negligible compared to the differences between stages , we generated scatterplots for the mature pollen stage ( supplementary fig . s7 , supplementary material online ) that was common to different experiments ( honys and twell 2004 ; schmid et al . 2005 ; borges et al . scatterplots showed an expression profile that was similar between experiments with rma normalization over all experiments and when normalized independently ( supplementary fig . s7b and c , supplementary material online ) and also showed more variation between expression levels when compared with nonnormalized and is normalized expression ( supplementary fig . scatterplots between nonnormalized experiments and between is normalized experiments showed less variation in expression levels , but in general , the correlations between expression levels from different experiments were highly independent from the normalization method . for rice and soybean , all experiments were kept for normalization . gene expression data for p. 
 patens for mature gametophyte , early- and midsporophyte ( odonoghue et al . 2013 ) were downloaded from geo ( gse32928 ) and the array and genome annotation ( v1.6 ) was obtained from www.cosmoss.org/physcome_project/wiki/downloads . in this data set , two samples per chip are hybridized , each with a different fluorescent dye ( green cy3 and red cy5 ) . a transcriptome index is the average of an evolutionary parameter like gene age ( tai ) , divergence ( tdi ) , and diversity ( tpi ) that is weighted by the expression level of each gene . confidence intervals were obtained by bootstrapping , using 100 sets of genes for each experimental stage . for estimates of gene age we used one - way blast ( default parameters ) hits against a set of genomes that are assigned to certain phylostrata and the blast hit to the most distant phylostratum defines the gene age ( alb and castresana 2007 ) . the oldest genes have a gene age value of 1 and the highest gene age value was assigned to genes that are specific to a given species ( youngest genes ) . for a. thaliana we classified 13 phylostrata , 9 for rice , 15 for soybean , and 5 for p. 
 patens . altogether we used 40 plant genomes , details about the hierarchical order , the genomes assigned to each phylostratum and number of genes with assigned gene age can be found in supplementary figure s8 , supplementary material online . for each species the largest age category was gene age of value 1 . to calculate a per gene estimate of divergence , we calculated dn / ds using pairwise alignments of homologous genes identified by inparanoid from the whole - genome comparison with its respective outgroup ( remm et al . we obtained per gene estimates of dn / ds ( = ka / ks ) estimates for genes specific to species pairs with the kaks_calculator ( zhang et al . tpis=i=1n(pn / n/((ps+1)/s)))eisi=1neis , where s is the developmental stage , n the number of genes , eis the expression intensity of gene i in developmental stage s , pn and ps the numbers of nonsynonymous and synonymous polymorphisms , respectively , and n and s are the numbers of nonsynonymous and synonymous sites , respectively . we used the ratio of nonsynonymous per site polymorphisms to synonymous per site polymorphism to estimate the distribution of fitness effects . higher values of pn / ps reflect an excess of slightly deleterious mutations ( keightley and eyre - walker 2007 ) . for technical reasons we used ps+1 rather than ps as suggested by stoletzki and eyre - walker ( 2011 ) because some genes have no synonymous polymorphisms and therefore would need to be excluded from the analysis which is biased ( stoletzki and eyre - walker 2011 ) . for compactness we refer to the term pn / n/((ps+1)/s ) as pn / ps throughout the manuscript . we tested whether transcriptome indices are different between stages by bootstrapping 100 samples of each index per stage and then performing a two - sample t - test to test for the differences in the means of bootstrapped values . if not noted otherwise , only the highest p value in the comparison of stages is reported . we calculated the weighted median transcriptome index of an evolutionary parameter x and assumed that i=1nei=1 . the weighted median of the evolutionary index is then xf with f such that 
 i < fei  < 1/2  and  i > fei  1/2 
  the standardized transcriptome index that does not consider genes with a nonsignificant expression ( supplementary fig . s4 , supplementary material online ) was calculated as follows : 
 t(x)is=i=1nxieisi=1neisx , 
 where x is the arithmetic mean of x1,  ,xn and n the number of significantly expressed genes . we further obtained per gene neutrality index ( ni ) for a. thaliana as follows : 
 ni = dspndnps 
 where ps=(ps+1)/s . the number of protein interactions for a. thaliana were obtained from the arabidopsis interactome database ( ftp://ftp.arabidopsis.org/home/tair/proteins/protein_interaction_data/interactome2.0/ ) . we obtained the genome sequences of a. thaliana ( arabidopsis genome initiative 2000 ) , rice ( o. 
 sativa , international rice genome sequencing project 2005 ) , and soybean ( g. 
 max , schmutz et al . 2010 ) from the plant genome database ( duvick et al . 2008 ) and the plant duplication database ( lee et al . 2013 ) along with their outgroups a. 
 lyrata ( hu et al . 2011 ) , sorghum bicolor ( paterson et al . 2009 ) , and phaseolus vulgaris ( schmutz et al . 2014 ) , respectively . polymorphism data were obtained from 80 a. 
 thaliana accessions ( cao et al . 2011 ) . to identify coding single nucleotide polymorphism ( snp ) information for rice we used the rice haplotype map project database ( second generation , http://www.ncgr.ac.cn/ricehap2/index.html ) and soybean we used snp information deposited in snpdb ( sherry et al . gene expression data were obtained for the three plants species from the plexdb ( dash et al . in particular , we focused on development stages preceding gametogenesis , during gametogenesis and embryogenic developments ( table 1 and supplementary file s1 , supplementary material online ) . 2003 ) and invariant set ( is ) methods were performed with the affy bioconductor package to normalize all data sets simultaneously . scatterplots of expression between replicates showed better results for rma normalization ( data not shown ) . therefore , unless stated otherwise , expression data shown in this study are based on a normalization across experiments using rma with log2 transformation . since different laboratory conditions can affect expression patterns ( massonnet et al . 2010 ) , we controlled for these effects in the a. thaliana data ( schmid et al . 2005 ) by removing data sets that were obtained from plants with different growth conditions before rna extraction ( supplementary file s1 , supplementary material online ) . to check whether the differences in expression between experimental conditions were negligible compared to the differences between stages , we generated scatterplots for the mature pollen stage ( supplementary fig . s7 , supplementary material online ) that was common to different experiments ( honys and twell 2004 ; schmid et al . 2005 ; borges et al . scatterplots showed an expression profile that was similar between experiments with rma normalization over all experiments and when normalized independently ( supplementary fig . s7b and c , supplementary material online ) and also showed more variation between expression levels when compared with nonnormalized and is normalized expression ( supplementary fig . scatterplots between nonnormalized experiments and between is normalized experiments showed less variation in expression levels , but in general , the correlations between expression levels from different experiments were highly independent from the normalization method . for rice and soybean , gene expression data for p. 
 patens for mature gametophyte , early- and midsporophyte ( odonoghue et al . 2013 ) were downloaded from geo ( gse32928 ) and the array and genome annotation ( v1.6 ) was obtained from www.cosmoss.org/physcome_project/wiki/downloads . in this data set , two samples per chip are hybridized , each with a different fluorescent dye ( green cy3 and red cy5 ) . a transcriptome index is the average of an evolutionary parameter like gene age ( tai ) , divergence ( tdi ) , and diversity ( tpi ) that is weighted by the expression level of each gene . confidence intervals were obtained by bootstrapping , using 100 sets of genes for each experimental stage . for estimates of gene age , we followed the procedure of quint et al . we used one - way blast ( default parameters ) hits against a set of genomes that are assigned to certain phylostrata and the blast hit to the most distant phylostratum defines the gene age ( alb and castresana 2007 ) . the oldest genes have a gene age value of 1 and the highest gene age value was assigned to genes that are specific to a given species ( youngest genes ) . for a. thaliana we classified 13 phylostrata , 9 for rice , 15 for soybean , and 5 for p. 
 patens . altogether we used 40 plant genomes , details about the hierarchical order , the genomes assigned to each phylostratum and number of genes with assigned gene age can be found in supplementary figure s8 , supplementary material online . for each species the largest age category was gene age of value 1 . to calculate a per gene estimate of divergence , we calculated dn / ds using pairwise alignments of homologous genes identified by inparanoid from the whole - genome comparison with its respective outgroup ( remm et al . we obtained per gene estimates of dn / ds ( = ka / ks ) estimates for genes specific to species pairs with the kaks_calculator ( zhang et al . tpis=i=1n(pn / n/((ps+1)/s)))eisi=1neis , where s is the developmental stage , n the number of genes , eis the expression intensity of gene i in developmental stage s , pn and ps the numbers of nonsynonymous and synonymous polymorphisms , respectively , and n and s are the numbers of nonsynonymous and synonymous sites , respectively . we used the ratio of nonsynonymous per site polymorphisms to synonymous per site polymorphism to estimate the distribution of fitness effects . higher values of pn / ps reflect an excess of slightly deleterious mutations ( keightley and eyre - walker 2007 ) . for technical reasons we used ps+1 rather than ps as suggested by stoletzki and eyre - walker ( 2011 ) because some genes have no synonymous polymorphisms and therefore would need to be excluded from the analysis which is biased ( stoletzki and eyre - walker 2011 ) . for compactness we refer to the term pn / n/((ps+1)/s ) as pn / ps throughout the manuscript . we tested whether transcriptome indices are different between stages by bootstrapping 100 samples of each index per stage and then performing a two - sample t - test to test for the differences in the means of bootstrapped values . if not noted otherwise , only the highest p value in the comparison of stages is reported . we calculated the weighted median transcriptome index of an evolutionary parameter x and assumed that i=1nei=1 . the weighted median of the evolutionary index is then xf with f such that 
 i < fei  < 1/2  and  i > fei  1/2 
  the standardized transcriptome index that does not consider genes with a nonsignificant expression ( supplementary fig . s4 , supplementary material online ) was calculated as follows : 
 t(x)is=i=1nxieisi=1neisx , 
 where x is the arithmetic mean of x1,  ,xn and n the number of significantly expressed genes . we further obtained per gene neutrality index ( ni ) for a. thaliana as follows : 
 ni = dspndnps 
 where ps=(ps+1)/s . the number of protein interactions for a. thaliana were obtained from the arabidopsis interactome database ( ftp://ftp.arabidopsis.org/home/tair/proteins/protein_interaction_data/interactome2.0/ ) .
reproductive traits in plants tend to evolve rapidly due to various causes that include plant - pollinator coevolution and pollen competition , but the genomic basis of reproductive trait evolution is still largely unknown . to characterize evolutionary patterns of genome wide gene expression in reproductive tissues in the gametophyte and to compare them to developmental stages of the sporophyte , we analyzed evolutionary conservation and genetic diversity of protein - coding genes using microarray - based transcriptome data from three plant species , arabidopsis thaliana , rice ( oryza sativa ) , and soybean ( glycine max ) . in all three species a significant shift in gene expression occurs during gametogenesis in which genes of younger evolutionary age and higher genetic diversity contribute significantly more to the transcriptome than in other stages . 
 we refer to this phenomenon as 
 evolutionary bulge  during plant reproductive development because it differentiates the gametophyte from the sporophyte . 
 we show that multiple , not mutually exclusive , causes may explain the bulge pattern , most prominently reduced tissue complexity of the gametophyte , a varying extent of selection on reproductive traits during gametogenesis as well as differences between male and female tissues . 
 this highlights the importance of plant reproduction for understanding evolutionary forces determining the relationship of genomic and phenotypic variation in plants .