In situ dissecting the evolution of gene duplication with different histone modification patterns based on high-throughput data analysis in Arabidopsis thaliana

Identification of duplicate genes with different modes in A. thaliana

The MCScanX-transposed software is commonly used to identify gene duplications that occurred within different areas based on the MCScanX algorithm which can detect multiple gene collinearity within and between related genomes (Wang et al., 2013). We applied MCScanX-transposed software with default parameters to identify duplicated genes with different modes, including WGD, proximal, tandem and transposed duplications in the A. thaliana genome. The unidentified genes were considered as singleton genes.

Histone profiling

Public ChIP-chip data (Dataset S2) were originally mapped to the TAIR10 version of the A. thaliana reference genome. A total of 10 distinct HM marks including H2Bub, H3K4me1, H3K4me2, H3K4me3, H3K9ac, H3K9me2, H3K27me1, H3K27me 3, H3K36me3 and H3K14ac were analyzed (Additional file 4). For each of the ten HM marks in the open reading frame regions (ORF) and upstream regulatory (promoter) regions (500 bp upstream of the transcription start site), we calculated the log2-transformed average enrichment ratio with nucleosome normalizing. The ratios were then transformed to z scores using the formula Z_X = (χ − μ)/δ, where χ is the ratio value for a gene, μ is the mean ratio of all genes, and δ is the standard deviation of this ratio for all genes. The HM divergence of a duplicate gene pair was measured by D_HM = 1− r, where r is the Pearson correlation coefficient of the HM profiles for the duplicated gene pair.

Randomized controlled analysis

Randomization tests examine data in a nonparametric way, without the need to make assumptions about object populations. The most important advantage of randomization tests is that they allow users to specify any test statistic of interest and remain statistically valid. To test whether the duplication group of interest is significantly different from others, the randomization test is applied.

For each test, 10,000 independent randomizations are performed for the controlled analysis. Note that, (i) to obtain the randomized control of HM divergence between duplicate pairs, we randomly generated the same number of singleton pairs. (ii) To test the chromosome bias on the HM divergence between duplicate pairs, randomized singleton pairs were chosen where both copies were located on the different chromosomes. Then, the Wilcoxon rank-sum test was used to evaluate the significance of distribution between the randomized control and observed data sets.

Statistical method

To compare the chromosome bias of HM divergence in ORF and promoter regions of different modes of gene duplication, we calculated the D_HM and compared the mean values of D_HM in each mode. Wilcoxon rank-sum test was applied to calculate the significance. The same method was applied to assess whether the patterns of expression divergence of duplicated genes were consistent across different perturbation conditions or not.

K_S and K_A calculation

Coding sequence divergence was measured by the rate of synonymous substitutions (K_S) and nonsynonymous substitutions (K_A). K_S and K_A between duplicated gene pairs were estimated using MCScanX-transposed software with default parameters.

Gene expression analysis

A total of 88 RNA-seq samples (Dataset S1) from 16 publications were divided into three sample conditions based on different perturbation conditions: (i) the normal development group, named as “Normal” containing 18 wild-type samples; (ii) the group of samples with environmental stresses, named as “Stress” containing 20 samples; (iii) the group of samples with mutations on chromatin modifier (CM), named as “CM_mu” containing 50 samples. These RNA-seq datasets were re-analyzed according to the protocol proposed by Trapnell et al. (2012). This protocol used TopHat2, which is a fast splice junction mapper for RNA-Seq reads (Kim et al., 2013) and Cufflinks2, which could assemble transcripts, estimate their abundances, and test for differential expression based on aligned RNA-Seq reads (Trapnell et al., 2010). Here, the raw datasets were downloaded from the SRA database according to the corresponding SRA ID listed in Additional file 1. Then the reads were aligned to TAIR10 using TopHat2. Finally, the gene expressed abundances were calculated using Cufflinks2.

Functional analysis

To understand the biological function of duplicate genes, the Database for Annotation, Visualization and Integrated Discovery (DAVID) tool with version 6.7 (Huang, Sherman & Lempicki, 2009) was applied to identify the enriched Gene Ontology (GO) functional terms based on each mode of the gene duplication. Further, REVIGO, which could summarize the long lists of GO terms by removing redundant terms (Supek et al., 2011), was used to remove redundant GO terms and visualize the GO enrichment results. All genes and its annotated standard GO terms of arabidopsis are listed in Additional file 5.

Identifying different modes of duplicate genes

All A. thaliana genes were classified by four modes of gene duplication including WGD, proximal, tandem and transposed duplication, and singleton genes (“Materials and Methods”). According to the definition, 8,111 duplicate gene pairs were identified including 4,293 (53%) WGD gene pairs, 2,130 (26%) tandem duplicate gene pairs, 784 (10%) proximal duplicate gene pairs and 904 (11%) transposed duplicate gene pairs (Dataset S2). As shown in Fig. 1A, duplicate genes are scattered across the whole genome, whereas only transposed duplicate genes tend to be located in transposable elements (TEs)-enriched region.

We then calculated the enrichment levels of each HM mark for each gene within its ORF and promoter region (“Materials and Methods”, Figs. S1 and S2). Comparing with singletons, WGD genes were more likely to be targeted by active marks (including H3K4me1/2/3, H3K9ac, H3K14ac and H3K36me3) and H3K27me1/3, instead of H3K9me2 in both ORF and promoter regions (Wilcoxon rank-sum test: p-value < 0.05 for all cases). Additionally, transposed duplicate genes show a higher level of H3K9me2 than other modes of duplicate genes and singletons (Wilcoxon rank-sum test: p-value < 0.05 for all cases).

Then, we calculated the pairwise Pearson correlations of HM marks for paralog pairs. The result showed that the correlation coefficients between the duplicate gene pairs for each mark are positive, except for the coefficients between WGD pairs for H3K9me2 (Table 1). Moreover, these correlation coefficients in promoter regions are larger than in ORF regions. On the other hand, some histone modification marks showed co-occurrence between duplicate gene pairs, and the overall patterns of histone modification marks between duplicate pairs are consistent (Pearson correlation coefficient: 0.418 < r < 0.991 and p-value < 1.0E−4, Mantel test (Mantel, 1967), Figs. 1B–1F and Fig. S3).

Divergence of histone modification profiles among different modes of duplicate genes

We further examined the divergence of HM profile between paralog pairs and singletons. The divergence of HM profiles was calculated (“Materials and Methods”) in both ORF and promoter regions, which were denoted by scores: D_HM–O and D_HM–P, respectively. To compare the HM divergence between paralog pairs and randomized singleton pairs, we randomly selected 8,111 pairs from singletons with 10,000 repeats. We found that both D_HM–O and D_HM–P are lower between duplicate pairs than that of randomized pairs (Wilcoxon rank-sum test: p-value < 2.2E−16 for both cases), which is consistent with the previous study in yeast (Zou et al., 2012).

To investigate the HM status of different gene duplication regions, we compared the effects of HM profile on different gene regions of each gene duplication mode. The results showed that the divergence of HM profile both in ORF and promoter regions are distinct among the different modes of gene duplication (Kruskal-Wallis Test, p-value < 2.2E−16; Figs. 2A and 2B). Notably, the transposed duplicates have a higher divergence score than other duplicate modes, meanwhile, tandem and proximal duplicates have lower divergence score.

The physical distance between single-gene duplications followed a trend: transposed duplicates > proximal duplicates > tandem duplicates. We assumed that position effects may have an influence on the divergence of HM profiles of genes with different duplication modes. To test this hypothesis, we selected all A. thaliana gene pairs on the same chromosomes and computed the correlation between HM divergences and physical distances of gene pairs both in ORF and promoter regions. Notably, the HM divergence of gene pairs on the same chromosome is not significantly correlated with physical distance (Pearson correlation coefficient: r = 0.066 for ORF regions and r = 0.089 for promoters). However, gene duplicate pairs on the same chromosome share a more similar HM profile than that on different chromosomes in ORF and promoter regions (Wilcoxon rank-sum test: p-value < 2.2E−16 for both cases, Figs. 2C and 2D).

Notably, duplicate pairs tend to be located on the same chromosome comparing with randomized singleton pairs (Chi-squared test: χ² = 152.0918, df = 1, p-value < 2.2E−16), which is termed as “chromosome bias”. This falls in line with the research in yeast (Zou et al., 2012). We then wanted to know whether chromosome bias gives a possible explanation for less difference of HM divergence between the duplicate gene pairs than that between randomized singleton pairs. To rule out the possibility, the same number of duplicate gene pairs and random pairs were selected in which both copies located on different chromosomes, with a result that HM divergence was lower between paralog pairs than that between randomized singleton pairs (Fig. 2E, Student’s t-test: p-value < 0.001 for both ORF and promoter regions).

Histone modification divergence of paralog pairs along with evolutionary time

To explore how HM profile evolves following gene duplication, the relationship of the HM profile of duplicate genes and the divergence of their coding sequence was characterized. In our study, coding sequence divergence was measured by synonymous substitutions (K_S) rate and nonsynonymous substitutions (K_A) rate. The duplicate pairs with K_A < 0.5 and K_S < 2.0 were selected for further analysis. We related the divergence of HM profile of duplicate genes to K_A (Figs. 3A and 3B) and K_S (Figs. 3C and 3D) using a scatterplot. We observed that D_HM–O and $D_{HM-P}$ are weakly positively correlated with K_S and K_A (Pearson correlation coefficient: 0.06 < r < 0.19; Fig. 3E). As shown in Fig. 3E, HM profile of duplicate genes in ORF diverges much quicker than that in promoters regions (Wilcoxon rank-sum test: p-value < 6.808E−10 for all gene modes), and $D_{HM-O}$was more correlated with K_S and K_A (Pearson correlation: r = 0.18 and 0.19) than $D_{HM-P}$ (Pearson correlation coefficient: r = 0.06 and 0.12). However, there is no significant difference in HM profile divergence between in ORF and promoter regions for randomized singleton genes (Wilcoxon rank-sum test: p-value = 0.8768). The result illustrated that the divergence of HM in ORF regions is more similar to the coding sequences divergence between paralog pairs.

We further explored the association of duplication modes with coding sequence divergence D_HM–O of WGDs tend to have a lower correlation with K_S and K_A than single-gene duplications (Fig. 3F; Fig. S4). Moreover, at similar K_S or K_A levels, transposed duplicate pairs were more likely to have higher HM divergence than the tandem and proximal duplicate gene pairs (Figs. 3A–3D). The different extent of HM divergence may be explained by the hypothesis that collinear duplicates tend to have similar chromatin environments, while transposed duplications re-locate to new chromosomal positions that usually have different chromatin environment (Wang et al., 2013).

Relationship between epigenetic divergence and gene expression divergence between duplication gene pairs

Our recent study revealed that epigenetic marks play a key role in determining transcriptional outcomes (Wang et al., 2015b). We thus asked whether the divergence of epigenetic patterns contributes to the expression divergence of paralogs. Correlation analysis was performed for HM divergence ( $D_{HM-P}$ and $D_{HM-O}$) and expression divergence (E) between paralog pairs, and we found that there is indeed a relation between them (Pearson correlation coefficient: r = 0.265 for $D_{HM-O}$ and r = 0.113 for $D_{HM-P}$, Figs. 4A and 4B). However, it is a weaker correlation between expression divergence and epigenetic divergence in promoter regions than that in ORF regions. Furthermore, the correlation differs among different modes of gene duplication (Pearson correlation coefficient: 0.239 < r < 0.368 for $D_{HM-O}$ and 0.096 < r < 0.241 for $D_{HM-P}$, Fig. S5). For proximal and transposed duplicates, there is a stronger correlation between HM divergence and expression divergence.

A previous study has suggested that histone modification enzymes and environmental stresses influence the expression evolution of duplicate genes in yeast (Zou et al., 2012). To test whether the expression divergence patterns of duplicate genes were consistent across different perturbation conditions in A. thaliana, we divided expression profiles into three types: Normal, Stress and CM_mu. The expression divergence between duplicate pairs under these three types of conditions were denoted as E_Normal, E_Stress and $E_{CM\mathrm{\_}mu}$, respectively. Our results showed that the E_{CM_nu} is significantly greater than $E_{Normal}$ and E_Stress (Fig. 4C, Wilcoxon rank-sum test: p-value < 2.2E−16). In contrast, there is no significant difference in expression divergence of duplicate genes in “Normal” and “Stress” conditions (Fig. 4C, Wilcoxon rank-sum test: p-value = 0.6097). As shown in Fig. 4C, different perturbation conditions have a similar effect on WGD and tandem duplicated genes with all duplicate genes but have little or no effect on the proximal and transposed duplicate genes. Taken together, epigenetic related enzymes (chromatin modifiers) may affect the expression evolution of duplicate genes, especially for WGD and tandem duplicated genes.

Functional divergence and epigenetic divergence between duplicate pairs

The previous study suggested that duplicate pairs with different biological functions differ significantly in ORF HM and promoter HM profiles in yeast (Zou et al., 2012). Here, we estimated how different biological functions influence the epigenetic divergence between duplicate pairs in A. thaliana. The biological functional divergence is classified into five models based on the biological functions of the duplicate genes (Fig. 5A). For each duplicate gene pair (using A and B as an example), A and B share the same GO terms (Type I); A and B partially share GO terms (Type II); A and B do not share GO terms, but they both have associated GO terms (Type III), or just A has associated GO terms (Type IV), or both A and B do not have associated GO terms (Type V). Therefore, the extent of functional divergence showed the following trend: Type I ≈ Type V > Type II > Type III ≈ Type IV.

The relationship between functional divergence and duplication modes was investigated. We found that the WGD gene pairs were mostly enriched in Type II, and WGD gene pairs have higher functional divergence than single-gene pairs (Fisher test: p-value < 1.439E−14, Fig. 5B). The single-gene duplication shows the following functional divergence trend: Transposed > Proximal > Tandem (Fig. 5B).

We then investigated the relationship between functional divergence and HM divergence. Duplicate genes from Type I and Type V have lower HM divergence levels than Type II, III and IV (Figs. 5C and 5D; Wilcoxon rank-sum test: p-value < 0.005 for all cases). Type V duplicate genes have much lower HM divergence than Type I in both ORF and promoter regions (Wilcoxon rank-sum test: p-value < 0.005 for both cases). Besides, except for Type V, other types have higher HM divergence in promoter regions than in ORF regions. The results suggested duplicate genes with distinct biological processes undergo the different evolutionary rates of HM pattern in their ORF and promoter regions. That is, functional divergence is consistent with HM divergence.

To understand the biological function of duplicate genes, the functions of the duplicate gene with high HM divergence ( $D_{HM-O}$ > 1 or D_HM–P > 1) were investigated. Using gene set enrichment analysis, eight functional models were identified: RNA biogenesis, translational regulation, protein catabolism, protein modification, transport, defense response, metabolism and transcriptional regulation (Fig. 6). It has been reported that the duplication of light-harvesting complex (LHC) in plants could diverge and acquire novel functions, such as response to various stress conditions including high light, high salinity, and nutrient limitation (Rochaix & Bassi, 2019). The biological functions associated with nutrient stress response, magnesium ion transport and phosphate ion transport were identified, indicating the validity of the results in this study.