High rates of evolution preceded shifts to sex-biased gene expression in Leucadendron, the most sexually dimorphic angiosperms

Abstract

Differences between males and females are usually more subtle in dioecious plants than animals, but strong sexual dimorphism has evolved convergently in the South African Cape plant genus Leucadendron. Such sexual dimorphism in leaf size is expected largely to be due to differential gene expression between the sexes. We compared patterns of gene expression in leaves among 10 Leucadendron species across the genus. Surprisingly, we found no positive association between sexual dimorphism in morphology and the number or the percentage of sex-biased genes (SBGs). Sex bias in most SBGs evolved recently and was species specific. We compared rates of evolutionary change in expression for genes that were sex biased in one species but unbiased in others and found that SBGs evolved faster in expression than unbiased genes. This greater rate of expression evolution of SBGs, also documented in animals, might suggest the possible role of sexual selection in the evolution of gene expression. However, our comparative analysis clearly indicates that the more rapid rate of expression evolution of SBGs predated the origin of bias, and shifts towards bias were depleted in signatures of adaptation. Our results are thus more consistent with the view that sex bias is simply freer to evolve in genes less subject to constraints in expression level.

Keywords: Leucadendron; chromosomes; evolution; evolutionary biology; gene expression; plant; proteaceae; sex-biased gene expression; sexual dimorphism.

Figure 1.. Strong morphological sex differences of the leaves evolved repeatedly in the genus Leucadendron, yet sex-biased expression affects only few genes, and the transcriptomes of males are not similar between species, nor are transcriptomes of females similar between species.

(A) Typical male (left) and female (right) shoot tips of L. rubrum, a wind-pollinated species with strong sexual dimorphism in leaves, stems, and inflorescence. (B) Left: supermatrix species tree; scale bar indicates expected number of substitutions per site. All branches showed full Shimodaira–Hasegawa-like support. Middle: schematic outlines of male (blue) and female (red) leaves, drawn to scale from a single example photograph; background shading indicates lower sexual dimorphism as light grey and higher sexual dimorphism as dark grey, as classified by the data of Tonnabel et al., 2014. Right: bar plots showing the percentage of male-biased (blue) and female-biased (red) genes among all expressed genes.

Figure 2.. Percent of sex-biased genes in leaves as a function of morphological sexual dimorphism (ratio of female over male leaf area) in 10 species of Leucadendron.

Male-biased expression is in blue, and female-biased expression is in red points resp. lines. The underlying data are found in Supplementary file 1 – Table S3.

Figure 2—figure supplement 1.. Sex-biased gene expression as a function of morphological sexual dimorphism in the leaves of 10 species of Leucadendron, shown in 6 alternative scatter plots for all combinations of 2 sexual dimorphism measures and 3 statistics of sex-biased gene expression.

Male-biased expression is in blue, and female-biased expression is in red points resp. lines. The underlying data are found in Supplementary file 1 – Table S3.

Figure 2—figure supplement 2.. Heatmap and cluster dendrogram for biological processes putatively regulated by sex-biased genes in 10 Leucadendron species.

Values are the proportions of sex-biased genes annotated with a given biological process (functional categorization according to TAIR database). The cluster dendrogram is based on Pearson’s correlation coefficient. Species names are abbreviated as follows: L. brunioides, brun; L. dubium, dubi; L. ericifolium, eric; L. linifolium, lini; L. muirii, muir; L. olens, olen; L. platyspermum, plat; L. rubrum, rubr; L. spissifolium, spis; L. xanthoconus, xant.

Figure 3.. Gene expression heatmap and hierarchical clustering dendrogram for the sexes of 10 Leucadendron species and the hermaphroditic Leucospermum outgroup, for sex-biased genes only (650 genes).

Gene expression values (columns) are the mean log₂(TMM − TPM + 1) per species and sex. The clustering of groups (rows) is based on distances calculated as 1 − Pearson’s correlation coefficient.

Figure 4.. Summary of evolutionary histories inferred for sex-biased gene expression (SBGE) in Leucadendron.

Left: species tree annotated with inferred numbers of sex-biased genes at ancestral nodes (bold), and gains of sex-biased genes annotated on each branch; no losses were inferred. Right: table marking the sex-bias status for the 63 genes that showed sex bias in more than 1 species (either shared bias in the same direction, or divergent sex bias). For details, see Supplementary file 1 – Table S4.

Figure 5.. Molecular sequence evolution measured as omega (dN/dS) for sex-biased and unbiased genes of Leucadendron.

(A) Omega estimated over the deep evolutionary timescale between Arabidopsis thaliana and the genus Leucadendron. Violin plots (horizontal bar marks the median, dots show the mean) for the three categories of bias status, with different sets of genes in each category. Mean omega was not significantly different between unbiased and either sex-biased category (permutation tests p > 0.05). (B) Omega estimated over the more recent evolutionary timescale between species of Leucadendron. Histograms show pairwise differences in omega for genes observed under both male bias and unbiased (i.e., in different species). The mean difference in omega is indicated by a vertical line, and does not deviate from zero significantly (permutation test, p = 0.8). (C) The same as (B) but for genes observed under both female bias and unbiased, again the mean difference in omega is not deviating from zero (permutation test, p = 0.9).

Figure 5—figure supplement 1.. Distributions of fitness effects (DFEs) of sex-biased genes (SBGs; male bias, blue; female bias, red) and unbiased genes (grey) in two Leucadendron populations.

Bars show the Akaike-information criterion (AIC)-weighted point estimates under four different DFE models, and error bars show the 95% confidence intervals of 200 site-level bootstrap replicates, fitted to the four models in proportion of their AIC weight. (A) DFEs in L. dubium. The SBGs shown here were defined by differential gene expression tests in this species. (B) DFEs in the distantly related L. ericifolium. The SBGs shown are those of L. dubium, whereas L. ericifolium itself had almost no SBGs, and none overlapped with L. dubium. The unbiased genes here are largely shared with L. dubium.

Figure 6.. Sex-biased genes in Leucadendron have ancestrally and intrinsically higher rates of expression evolution.

(A) Histograms of mean absolute standardized phylogenetically independent contrasts (PICs) of gene expression for sex-biased (dark grey) and unbiased genes (light grey). The difference in means of the two categories is 2.3 (permutation, p = 2 × 10⁻⁵). Sex-biased expressions themselves were excluded when calculating the PICs. (B) Interspecific gene expression distance as a function of DNA sequence distance for 55 species pairs of dioecious Leucadendron and the hermaphrodite relative Leucospermum. Sex-biased expressions themselves were excluded when calculating gene expression distances. Gene expression distances are shown for three different categories: distance between species (mean over the sexes) for unbiased genes (black points and line), distance between males for sex-biased genes (blue points and line), and distance between females for sex-biased genes (red points and line). DNA distances > 4% pertain to Leucadendron–Leucospermum pairs. The shaded envelopes around linear regression lines represent parameter standard errors. DNA sequence distance was significantly correlated with each of the three categories of gene expression distance (Mantel tests, all p ≤ 0.0225).

Figure 6—figure supplement 1.. Demonstration that interspecific differences in expression level tend to correlate positively with the level of intraspecific expression variation, posing the risk that genes with noisy expression are mistakenly inferred to have fast rates of expression evolution.

This example scatter plot and correlation test shows all expressed genes in the species pair L. rubrum and L. dubium. The x-axis shows the intraspecific coefficient of variation over all 12 replicates per species, averaged over both species. The y-axis shows the interspecific log₂ fold-change between the mean expression levels of the two species. Similar results were obtained for all other pairs of species and also with entirely artificial, randomly simulated data (not shown). However, as described in the main text and shown in Supplementary files 1 and 2, this trend did not confound the inference of faster expression evolution for sex-biased genes in this study, because sex-biased genes showed low expression noise (i.e., expression variation within sexes).

Figure 6—figure supplement 2.. Challenging the inference of sex-biased genes by repeated permutation of the sexes.

In most species and sexes, observed counts of sex-biased genes are significantly greater than counts from permuted datasets, indicating that the observed sex-biased genes are unlikely to represent stochastic artefacts. 'Fbiased' refers to the number of female-biased genes, and 'Mbiased' denotes the number of male-biased genes. 'perm' stands for permutated datasets (light coloured bars in histograms), whereas 'obs' stands for the observed result from the unpermutated data (solid vertical lines). Sex-biased genes are defined under the conventional thresholds of 5% FDR and two-fold change between the sexes.

Figure 6—figure supplement 3.. Expression specificity (Shannon entropy) over different tissues and developmental stages for sex biased and unbiased genes in Leucadendron, using gene expression specificity data of homologous Arabidopsis thaliana genes (Klepikova et al., 2016).

Violins show the complete distributions, horizontal bars indicate the medians and dots the means. The difference of the means is 0.095 (permutation, p = 2 × 10⁻⁵).