Adaptive evolution driving the young duplications in six Rosaceae species

Abstract

Background

In plant genomes, high proportions of duplicate copies reveals that gene duplications play an important role in the evolutionary processes of plant species. A series of gene families under positive selection after recent duplication events in plant genomes indicated the evolution of duplicates driven by adaptive evolution. However, the genome-wide evolutionary features of young duplicate genes among closely related species are rarely reported.

Results

In this study, we conducted a systematic survey of young duplicate genes at genome-wide levels among six Rosaceae species, whose whole-genome sequencing data were successively released in recent years. A total of 35,936 gene families were detected among the six species, in which 60.25% were generated by young duplications. The 21,650 young duplicate gene families could be divided into two expansion types based on their duplication patterns, species-specific and lineage-specific expansions. Our results showed the species-specific expansions advantaging over the lineage-specific expansions. In the two types of expansions, high-frequency duplicate domains exhibited functional preference in response to environmental stresses.

Conclusions

The functional preference of the young duplicate genes in both the expansion types showed that they were inclined to respond to abiotic or biotic stimuli. Moreover, young duplicate genes under positive selection in both species-specific and lineage-specific expansions suggested that they were generated to adapt to the environmental factors in Rosaceae species.

Background

Gene duplications contribute to the generation of new genetic materials and novel gene functions, which drive the evolution and divergence of genomes and genetic systems [1, 2]. In plant genomes, the frequent occurrence of whole-genome duplications, segmental duplications, and polyploidizations results in masses of duplication loci [3, 4]. The whole-genome duplication (WGD), a sort of gene duplications sharply accelerates the scale of chromosome or the whole genome, but followed by a series of gene loss, gene conversion and so on [5, 6]. For tandem duplication, it might be caused by unequal crossing over leading to the progeny duplicates located adjacently to each other in a cluster intra-chromosome [6, 7]. The tandemly duplicate copies exhibit a coordinated expression mode and increase the divergence distance among themselves [7]. The transposon-related duplication or tansponson-mediated duplication is replicative transposition involved with transposable elements [6]. For example, in Oryza sativa (rice) and Arabidopsis, approximately 15–62% and 90%, respectively, of the gene loci are estimated to arise from gene duplication [8,9,10].

The large-scale existence of duplicate genes implies the retention and evolution of duplicates among plant genomes [5]. However, duplicate genes face three long-term fates: nonfunctionalization (or pseudogenization), characterized by one of the copies losing its function; neofunctionalization reflected by one of the copies gaining a novel function; and subfunctionalization exhibited by duplicate copies inheriting parts of the original gene function [5]. Nonfunctionalization/pseudogenization is the most widespread fate of the duplicate copies. However, neofunctionalization is the preservation mechanism to retain them, which is reflected by the positive selection during or after duplicate fixation [1, 11].

The signatures of positive selection acting on duplicate genes, commonly indicating that the duplicates are subject to adaptive evolution, were previously reported in plant genomes. For example, in Arabidopsis thaliana, the imprinted gene MEDEA (MEA) undergoes positive Darwinian selection along with neofunctionalization after duplication [12]; similarly, in Arabidopsis and a few grass species, the centromere protein C (CENP-C) genes with complex duplicate regions are under positive selective pressure [13]; the chalcone synthase (CHS) genes undergo positive selection in Dendranthema genomes [14]. Furthermore, a group of young duplicate genes that underwent adaptive evolution were detected in plants, such as the extremely expanded nucleotide-binding site leucine-rich repeat (NBS-LRR) genes of Vitis vinifera, Populus trichocarpa, and the Rosaceae species [15,16,17]. However, the evolutionary characteristics of young duplicate genes have been rarely reported at genome-wide levels among closely related plant species.

The whole-genome sequencing of Fragaria vesca, Malus x domestica, Pyrus communis, Prunus persica, Rosa chinensis, and Rubus occidentalis provides us an opportunity to investigate the evolution of the recent duplicate genes among the six Rosaceae genomes. The Rosaceae is a large family possessing high economical values, composed by four subfamilies, Spiraeoideae, Rosoideae, Maloideae, and Prunoideae. The six species involved in three subfamilies of Rosaceae, covering different evolutionary distances, containing Rosoideae (F. vesca, R. chinensis and R. occidentalis), Maloideae (M. x domestica and P. communis) and Prunoideae (P. persica). The origination of Rosaceae family is predicted during the Late Cretaceous [18]. Nine ancestral chromosomes existing in the ancestral Rosaceae genome, modern Rosaceae genomes are generated after a series of chromosome fission, fusion, and duplications during the evolutionary processes of Rosaceae family [19]. Especially, the genomes of M. x domestica and P. communis underwent a common recent WGD, but no similar large-scale duplication was reported in the other four genomes [20,21,22,23,24,25,26]. In our study, a genome-wide identification and genetic evolution analysis of young duplications were performed among the six diploid genomes. Our results demonstrated that the young duplicates underwent adaptive evolution for acclimatization in the six species.

Results

Young duplicate genes in the six Rosaceae species

A total of 35,936 gene families were explored across the six Rosaceae species containing 21,650 young duplicate gene families, which indicated that young duplications occurred in 60.25% of the total gene families (Table 1 and Additional file 1: Table S1). Species-specific and lineage-specific expansions were detected in these young duplicate gene families based on their duplication patterns. The total family number of species-specific expansions (14,988) outdistanced that of lineage-specific expansions (6662). In species-specific expansions, distinct family numbers were found among the six species, such as the most gene families in M. x domestica (6184), moderate number in P. communis (3122), and the least gene families (791) in R. occidentalis. Interestingly, in lineage-specific expansions, there was an extremely high value (6105) in the lineages of M. x domestica and P. communis, probably because of the close phylogenetic relationship between the two species and the common recent WGD shaping and increasing their genomes [21, 22]. Except in the lineages of M. x domestica and P. communis, a broad range of family numbers (1 to 149) were detected in lineage-specific expansions. The second largest gene number (149) was observed in the lineages of F. vesca and R. chinensis, which may be attributed to their close evolutionary relationship. The similar phenomenons were also found in the lineages of M. x domestica, P. communis and P. persica or R. chinensis and R. occidentalis (Table 1 and Additional file 1: Table S1).

Table 1 Number of young duplicate gene families for two types of expansions

For the families belonging to lineage-specific expansions, it is worth mentioning that seven young duplicate gene families included 156 gene members from the lineages of all six species. That is, each of the six species has two or more gene members in each of the seven gene families. To detect the species-specific duplication events in these families, two or more genes from one species clustered together in a clade (bootstrap values > 50) were marked as species-specific duplication events in the seven NJ trees (Additional file 5: Fig. S1). There were 9, 8, 3, 5, 7, 3, and 1 species-specific duplication events involving 15, 15, 6, 10, 11, 6, and 2 genes in family679, family730, family1336, family2291, family4459, family4952, and family5347, respectively (Additional file 5: Fig. S1). The results demonstrated that 65 genes (65/156 = 41.67%) were involved in species-specific duplications among the seven young duplicate gene families.

Duplication types of the young duplicate genes

The young duplicate genes could be classified into three duplication types, containing tandem duplication, transposon-related duplication and WGD, at the genome-wide level among the six Rosaceae species. In species-specific and lineage-specific expansions, young duplicate genes were involved in all the three duplication modes, but distinct gene numbers and percentages were displayed in different duplication types in the six species (Table 2). For example, there were relatively lower gene numbers and proportions in the three duplications types among species-specific and lineage-specific expansions in F. vesca and R. occidentalis.

Table 2 Gene numbers and percentages of young duplicate genes from three duplication types in species-specific and lineage-specific expansions

In species-specific expansions, the gene numbers of tandem duplications were much higher than those of the other two duplication types in F. vesca, M. x domestica, P. persica, R. chinensis and R. occidentalis. Accordingly, the highest percentages of the young duplicate genes came from tandem duplications were also detected in the five species. It was indicated that tandem duplications played important roles in the young duplications after the speciation of the five plants. Especially, 37.54% of the young duplicate genes were produced by tandem duplications in P. persica, representing the highest percentage compared with the proportions of this duplication type in the other species. However, in P. communis, the largest gene number and proportion were discovered in WGDs.

In lineage-specific expansions, young duplicate genes partly changed the distributions in the three duplication types compared with those of species-specific expansions. The largest gene numbers were detected in tandem duplications of F. vesca and R. chinensis, and their related proportions were 13.55 and 12.90% in the two species, respectively. More young duplicate genes were derived from WGDs in M. x domestica, P. communis and P. persica, and from transposon-related duplications in R. occidentalis. It is worth noting that relatively large percentages of young duplications were belonging to the WGDs in M. x domestica (39.99%) and P. communis (44.85%), but lower percentages of tandem duplicated genes in the two species (15.17% in M. x domestica and 12.07% in P. communis). The results illustrated WGDs driven the expansions of young duplicate genes in M. x domestica, P. communis and P. persica before the species differentiation and divergence. Therefore, all of these demonstrated that tandem duplications and WGDs might be the major force promoting the occurrence of young duplicate genes in the six Rosaceae species.

Domain preference of the young duplicate genes

The protein domains of the young duplicates were explored in the species-specific and lineage-specific expansions to uncover the functional preference of the duplicate genes among the six Rosaceae species.

A total of 2117 different domains were detected in the species-specific expansions among the six species (Additional file 2: Table S2). It is worth mentioning that 43.50% (921/2117) of the domains appeared in only one species, indicating that approximately one half of the protein domains were uniquely encoded by species-specific duplicate genes in the six species. On the contrary, only 5.15% (109/2117) of the domains occurred simultaneously in all the six species. Interestingly, the low-frequency domains were relatively low in number in all the species, while some of the high-frequency domains were high in number in the related species (Fig. 1). For example, there were many high-count domains, especially the domains of PPR, LRR, Pkinase, p450, and NB-ARC, shared by the species-specific duplicates of the six Rosaceae species.

Fig. 1

Top 20 protein domains of the young duplicate genes in species-specific expansions. The x-axis means the numbers of different domains. The y-axis means the domains taking the top 20 places of domain numbers. a: F. vesca, b: M. x domestica, c: P. communis, d: P. persica, e: R. chinensis and f: R. occidentalis

Although the numbers of domains found in lineage-specific expansions (2000) and in species-specific expansions were more or less equal, the domain frequency detected in both type of expansions was distinctly different. Clearly, only 5.20% of the protein domains (104/2000) were discovered in one species, such as B-lectin, Vicilin, and Trigger, demonstrating that a small amount of lineage-specific duplicate genes had exclusive domains in some species (Additional file 2: Table S2). In addition, 22.95% of the protein domains (459/2000) were found to co-occur in all the six species, with 7.56% (151/2000), 4.85% (97/2000), 25.30% (506/2000), and 34.15% (683/2000) of them appearing simultaneously in the lineages of five, four, three, and two species, respectively. Similar to the high-frequency domains of species-specific expansions, the domains of lineage-specific expansions also exhibited high occurrence in all the six species and also possessed a large number of copies in them, containing the Pkinase, PPR, LRR, p450, WD40, and Ribosomal, etc. (Fig. 2). Therefore, it may be concluded that the high-frequency duplicate domains in species-specific and lineage-specific expansions, involved in growth and development (Ribosomal, Ank, and Peptidase) or response to environmental stresses (PPR, NB-ARC, LRR, and Pkinase), might play a key role in the evolutionary processes of the six Rosaceae species.

Fig. 2

Top 20 protein domains of the young duplicate genes in lineage-specific expansions. The x-axis means the numbers of different domains. The y-axis means the domains taking the top 20 places of domain numbers. a: F. vesca, b: M. x domestica, c: P. communis, d: P. persica, e: R. chinensis and f: R. occidentalis

Duplication time of the young duplicate genes

The Ks values are molecular scales of duplication time and the divergence time. To further detect the timing of young duplication events in the six Rosaceae species, Ks values were calculated in both species-specific and lineage-specific duplicate gene families.

In species-specific expansions, the average Ks values of the orthologs were higher than those of the paralogs only in P. communis, R. chinensis and R. occidentalis (Table 3). However, the Ks values of paralogs obviously peaked at the range of 0 to 0.1 with extremely high frequency and slowly decreased from 0.1 to 1 in all species, except P. communis, in which the Ks values peaked at the range from 0.1 to 0.2 (Fig. 3). These results illustrated that a considerable portion of the young duplicate genes were generated at the very recent times. In the lineage-specific expansions, the orthologs had larger Ks values than paralogs, which suggested that species divergence was followed by duplication events. In addition, the Ks values distributed differently with lower frequency from 0 to 1 compared with those in species-specific expansions. For example, the peak values of Ks were in the range of 0–0.1 in F. vesca, M. x domestica, P. persica and R. occidentalis, 0.1 to 0.2 in P. communis and R. chinensis (Fig. 3). Although the peaks were still detected at 0 to 0.1 in four species, no extreme advantage in Ks frequency compared with those at the range of 0.1–0.2 or 0.2–0.3. The observation proved that, in the period of the recent time, much more species-specific duplicate genes were produced than the lineage-specific ones. Moreover, the appreciable clustering of the Ks values around 0.2 in M. x domestica and P. communis of lineage-specific expansions was consistent with the recent WGD in the two species [21, 22].

Table 3 Average Ks values and Pi values of young duplicate gene families for two types of expansions

Fig. 3

The Ks values of paralogs of young duplicate gene families in the two types of expansions. The x-axis means the range of Ks values from 0 to 1, and the range was divided into ten parts in unit of 0.1. The y-axis represents the occurrence of Ks value in each unit. a: F. vesca, b: M. x domestica, c: P. communis, d: P. persica, e: R. chinensis and f: R. occidentalis

The nucleotide diversity of the young duplicate genes

To deeper explore the evolutionary differences between paralogs and orthologs, we calculated the nucleotide diversity values (Pi values) among species-specific and lineage-specific duplicate genes (Table 3).

In species-specific expansions, the paralogs had larger average Pi values than the orthologs in each of the six species. Moreover, t-test analysis were also operated between the Pi values of paralogs and orthologs, showing Pi values of paralogs were significantly higher than those of orthologs in each of the six species (P < 0.01). The results manifested that copies derived from species-specific duplications (paralogs) might undergo a relative faster sequence divergence leading to the larger diversities among paralogs than orthologs in the six species.

However, the opposite results of paralogs with lower average Pi values than the orthologs were found in lineage-specific expansions of the studies species, except P. communis. The paralogs have significantly smaller Pi values than the orthologs in F. vesca, M. x domestica, P. persica, R. chinensis and R. occidentalis (t-test, P < 0.01). It could be inferred that the ancestor copies inherited from ancestor species to the studies species (orthologs) might be driven by a faster divergence speed after the lineage-specific duplications in the five species.

Selective pressure on young duplicate genes

The ratio of nonsynonymous to synonymous substitution (Ka/Ks) is an important indicator of the functional constraints on genes. Therefore, the Ka/Ks ratios of paralogs and orthologs were examined in all species-specific and lineage-specific duplicate gene families.

In both species-specific and lineage-specific expansions, most of the gene pairs with Ka/Ks ratios smaller than 1 illustrated that a majority of the young duplicate genes were subject to purifying selection among the six Rosaceae species. Nevertheless, a fraction of the gene pairs showed Ka/Ks ratios greater than 1, suggesting that they underwent positive selection (Fig. 4). In species-specific expansions, the paralogs had greater median and average values compared with the orthologs in the six Rosaceae species. Moreover, the Ka/Ks ratios exhibited highly significant differences between paralogs and orthologs in each species (t-test, P < 0.01), demonstrating that paralogs had significantly larger Ka/Ks values than orthologs in species-specific young duplicate gene families among the six Rosaceae species. A similar phenomenon was observed in lineage-specific expansions, the paralogs had highly significantly greater Ka/Ks values than the orthologs in all the six species (t-test, P < 0.01). These results indicated that paralogs were driven by weaker functional constraints and had faster evolutionary rates than orthologs in the young duplicate gene families of the six Rosaceae species.

Fig. 4

The Ka/Ks ratios of young duplicate genes in the two types of expansions. The box plots are exhibiting the distributions of Ka/Ks values among paralogs and orthologs in the two expansion types. The small square and the line in the box represent average and median values of the Ka/Ks values, respectively. a: F. vesca, b: M. x domestica, c: P. communis, d: P. persica, e: R. chinensis and f: R. occidentalis

Furthermore, the phenomena were more directly displayed by the linear analysis of Ka/Ks ratios between paralogs and orthologs from the same young duplicate family in species-specific and lineage-specific expansions (Additional file 6: Fig. S2). The paralogs had higher Ka/Ks values than the orthologs of the same family and are represented by the corresponding dots above the trend lines (blue lines: slope equal to 1). Therefore, the farther the dots were from the trend lines, the faster did the evolutionary rates occur in the related family. The protein domains of these families were examined, and it was found that some of them were connected with in response to biotic or abiotic stresses, including PPR, FAR1 and UDPGT (Additional file 3: Table S3).

Chromosomal location of young duplicate genes

The physical location of the young duplicate genes, in both the species-specific and lineage-specific expansions, was uneven on the chromosomes in the six Rosaceae species. Accordingly, the trends of the gene densities were basically consistent with those of gene numbers in the six species.

In species-specific expansions, four distribution patterns of the duplicate genes on the chromosomes in the six species were noticed (Additional files 7, 8, 9, 10, 11 and 12: Figs. S3–S8). In the first pattern, more duplicate genes preferred to distribute themselves in the regions near the two telomeres on each chromosome, such as in chromosomes 3, 5, and 6 of F. vesca and chromosomes 2, 3, 7, 8, 9, 10, 11, 14, 15, and 17 of M. x domestica. In the second pattern, the species-specific duplicate genes exhibited peak distributions in the neighborhood of one of the telomeres on each chromosome, such as in chromosomes 1, 4, 5, 6, 13, and 16 of M. x domestica. In the third pattern, there was a relative mean distribution of the young duplicate genes on the chromosomes, such as on all chromosomes of P. communis and R. occidentalis, chromosomes 3, 4, 5, and 6 of P. persica, and chromosomes 2, 4, and 5 of R. chinensis. In the fourth pattern, the peaks of the duplicate genes were located on the non-telomere regions of the chromosomes, such as in chromosomes 1 and 2 of P. persica and chromosomes 1, 6, and 7 of R. chinensis. Most of the chromosomes belonged to the first, third and forth patterns showed the gene densities consistent with the distribution pattern of the young duplicate numbers. That is, the chromosomal regions with gene number peaks generally presented gene density peaks, while the non-peak regions were also shared by the similar distributions between the gene numbers and gene densities, especially the related chromosomes from F. vesca and R. chinensis. Although similar distribution trends were revealed between the gene numbers and gene densities, there were some differences at some positions of the chromosomes. For instance, no number peak on the right telomere region of chromosome 4 from P. persica, but a density peak of young duplicate genes displayed here. This result demonstrated relatively lower genes distributed at this telomere region. The similar phenomenon was also exhibited in chromosome 1 of R. occidentalis.

Similarly, the four distribution patterns of the duplicate genes were also detected in lineage-specific expansions (Additional files 7, 8, 9, 10, 11 and 12: Figs. S3–S8). In the first pattern, a clear clustering around the two telomeres was found in chromosomes 2, 3, 7, 8, 9, 10, 11, 14, 15, and 17 of M. x domestica and chromosomes 3 and 6 of P. persica. In the second pattern, it was demonstrated by chromosomes 1, 4, 5, 6, 13, and 16 in M. x domestica and chromosomes 2, 4, 5, 7, and 16 in P. persica. The distribution tendency of the third pattern was revealed by all the chromosomes of F. vesca, P. communis, and R. occidentalis, and chromosomes 1, 2, 3, and 6 of R. chinensis. The fourth pattern was represented by chromosome 8 in P. persica and chromosomes 4, 5, and 7 in R. chinensis. Moreover, the positions with lower gene numbers but higher gene densities was also discovered in different chromosomes, such as chromosome 6 in M. x domestica and chromosome 12, 15 of P. communis. Interestingly, for the chromosomal locations of young duplicates in lineage-specific expansions, the extremely similar distributions between the gene numbers and gene densities were also discovered in all the chromosomes of F. vesca and R. chinensis. Therefore, the results indicated that the distribution patterns of young duplicate genes on the same chromosomes between the species-specific and lineage-specific expansions were similar on some of the chromosomes among the six species.

To further verify the relationship of gene locations between the two expansion types, we calculated the correlation coefficients of the young duplicate gene numbers between species-specific and lineage-specific expansions on the same chromosome in the six species (Additional file 4: Table S4). Most of the correlation coefficients of M. x domestica and some of them of P. communis were greater than 0.6, which indicated a strong correlation between the locations of species-specific and lineage-specific duplicate genes in the corresponding chromosomes.

Collinearity between young duplicate genes

For exploring the evolutionary history of the young duplicate genes, the collinearity block analysis was performed among the young duplicate genes of lineage-specific expansions co-occurred in the six species. The family members and their neighbouring genes predicted in the collinearity blocks were showed across the corresponding chromosomal regions (Additional file 13: Fig. S9). Although the collinearity blocks located in different chromosomes among the six species, the young duplicate genes of the four families with collinearity relationships were frequently found distributed in the chromosome 3 in F. vesca, chromosome 3, 5, 10 and 11 in M. x domestica and P. communis, chromosome 4 and 6 in P. persica, chromosome 5 in R. chinensis, and chromosome 3 in R. occidentalis, respectively. These indicated that these chromosomal regions might originate from the ancestral chromosomes of the common ancestor species of Rosaceae. In family679, one F. vesca member, two M. x domestica members, one P. communis member, one P. persica member, one R. chinensis member, and one R. occidentalis member were predicted locating in the same collinearity block, and the other shared collinearity block contained one, two, two, one, one and one members from F. vesca, M. x domestica, P. communis, P. persica, R. chinensis, and R. occidentalis. Similarly, there were also one or two members from each species involved in the same collinearity block in the other three families. Interestingly, some members and their neighbouring genes still represented the consistent transcriptional directions after the long evolutionary processes, such as the members from F. vesca (Fv3g35380) and R. chinensis (RC5G0545100) in familiy679. The family members along with the related neighbouring genes exhibited these ancestral regions, which were conserved and retained in the chromosomes of the six modern Rosaceae species after experiencing series of evolutionary events, including the young duplications.

Discussion

Rosaceae evolution history exerting influences on the young duplicate genes

To retrospect the history of the Rosaceae family, it could be better discovered the characteristics of young duplicate genes reflected by the evolutionary processes of the six species.

In the ancestor genome of Rosaceae family, nine protochromosomes represent the ancestral Rosacese karyotype, which experienced varying scales of fusions and fissions at chromosome level to successfully develop into the ancestral karyotype of subfamilies, such as ancestral Rosoideae karyotype and ancestral Prunoideae karyotype [19]. There are different evolutionary processes happen in ancestral Maloideae karyotype compared with the former two karyotypes, that is, duplications and large scales of chromosome fusions and fissions contribute to the generating of ancestral Maloideae one [19]. A recent WGD shared by the modern Rosaceae species of Malus and Pyrus is also reported, which might be one of the main factor leading to the large numbers of young duplicate families in the lineage of M. x domestica and P. communis (Table 1) [21, 22, 26]. In addition, this recent WGD might also explain the result, for which the young duplicate genes belonging to WGDs accounted for larger proportions in M. x domestica and P. communis from lineage-specific expansions (Table 2). On the contrary, the other four Rosaceae genomes lack the similar duplication observed in Malus and Pyrus species, but some fragmentary triplicated sub-genomic blocks are in P. persica genome [19, 20, 23, 25]. Therefore, these phenomena supported the relatively lower young duplicates and lower numbers of WGDs in F. vesca, R. chinensis and R. occidentalis, and a little higher number of WGDs in P. persica from lineage-specific expansions (Table 2). Besides, the duplication time of the young duplicate genes in M. x domestica and P. communis from the lineage-specific expansions (Fig. 3) partly coincided with the shared recent WGD in Malus and Pyrus species (Ks peak around 0.2) [21]. These demonstrated that the recent WGD might partly contribute to the generating of lineage-specific young duplicate genes in the two species.

Although different evolutionary pathways shaped distinct modern Rosaceae genomes, some shared ancestral regions/genes were still detected in their existing chromosomes among the six Rosaceae species (Fig. S9). These further indicated the studied species originated from the common ancestor, and still retained some ancestral genetic traits in their modern genomes [27].

Functional preference of young duplicates adaptation to environmental stresses

Duplicated genes were previously represented in response to different ecological conditions in various plants [1]. For example, the gene duplicates were adaptation to extreme temperature in Arabidopsis and Brassicaceae, salt stress in Citrus and Sorghum bicolor, drought in wild tomato, weedicides in Medicago sativa, Glycine max and Nicotiana tabacum, nutrient limitation in Poaceae, and pathogen infections in Arabidopsis and O. sativa [28,29,30,31,32,33,34]. In our study, the domains of PPR, NB-ARC, LRR, Pkianse, Ribosomal, Ank and Peptidase exhibited high frequency in species-specific and lineage-specific duplications, which were involved in the adaptation to the changing environmental forces in plants.

The pentatricopeptide repeat (PPR) protein is composed of the PPR motif, a degenerate 35-amino-acid repeat tandemly arranged for 2 to 26 copies [35]. In plants, the PPR gene family is one of the largest gene families because of its extreme expansions in copy number, which makes it broadly function as gene-specific transcriptional regulators in the mitochondria and chloroplasts in response to ecological factors [36]. A mitochondrial PPR protein, PPR40, has created a close connection between oxidative respiration and stress/hormonal responses in A. thaliana [37]. In addition, some of the PPR genes, such as the restorer-like genes, are subject to diversifying selection, associating with the “arm-race” co-evolution of nuclear and cytoplasmic genomes in plants [38,39,40]. The PPR genes going through duplication and divergence have some characteristics similar to those of R genes, and they may play a vital role in adaptive evolution by the modification of cytoplasmic proteins in plants [39, 40].

The NBS-LRR gene family is one of the largest disease resistance gene families in plants, which encodes the NBS-LRR proteins mainly composed of the NB-ARC and LRR domains [41]. Plant NBS-LRR proteins, ancient in origin and abundant in copy numbers, are involved in the broad-spectrum resistance against various pathogens including fungi, viruses, bacteria, oomycetes, insects, and nematodes [41, 42]. The presence of a Pi2/9 locus containing tandemly arranged NBS-LRR genes with a variety of copy numbers and high levels of divergence in the rice genome suggests that the genes of this locus adaptively cope with diverse pathogens [43]. In the Pi2/9 locus, genes Pi2 and Pi9 exhibit broad-spectrum resistance to rice blast [44, 45]. Similarly, many of the plant NBS-LRR genes are in the phase of developing a rapid evolutionary pattern driven by diversifying selection, which is accomplished by numerous duplications, frequent sequence exchanges, and pseudogene losses [46]. These evolutionary processes facilitate the fast divergence of NBS-LRR genes to adapt to the rapidly changing pathogens [46, 47].

The protein kinase (Pkinase) family is one of the largest superfamilies in eukaryotic organisms, composed of a mass of homologous proteins with common conserved structural motifs [48]. However, the members of the Pkinase family exhibited an abundant variation in their structures, regulation patterns, and substrate specificities and functions involved in many cellular processes such as development, differentiation, proliferation, and defense responses against abiotic or biotic stresses [48, 49]. Especially, many of the leucine-rich repeat receptor-like kinase (LRR-RLK) genes are part of the disease resistance gene family in plants [16]. The LRR-RLK genes involved in interactions with the environment are mainly detected in lineage-specific duplicate clades driven by positive selection in plant genomes [50, 51]. For example, an LRR-RLK gene, Xa21, from wild rice, O. longistaminata, showed broad-spectrum resistance to bacterial blight (Xanthomonas oryzae pv. oryzae) [52, 53]. Among the Xa21 homologous loci in gramineous genomes, the frequent occurrence of duplication, gene translocation, and copy number variation demonstrated that Xa21 homologs undergo rich diversification in order to adapt to the specific fast-changing pathogens over long evolutionary timescales [54].

Adaptive evolution of young duplicate genes

Gene families are driven by adaptive evolution following gene duplications, which is represented by duplicate genes under positive selection [1, 55]. In this study, a part of gene pairs from species-specific and lineage-specific expansions with Ka/Ks values greater than 1 indicated that the genes were driven by adaptive evolution in the six Rosaceae species. For example, NBS-LRR genes, a large family of disease resistance genes, always experience large-scale gene duplications and rapid adaptive evolution in many plant genomes, including several Rosaceae species [16, 17]. Based on the divergence analysis of the Ka/Ks (or dN/dS) ratios, it has been previously reported that many gene families in plants are driven by adaptive evolution. For example, the monosaccharide transporter (MST) gene family, an ancient origin and largely duplicate family, reveals evidence of positive selection in rice and Arabidopsis [56]; the chalcone synthase (CHS) multigene family, which plays an important role in flavonoid biosynthesis and the divergence of its three subfamilies (SF1, SF2, and SF3), has been under positive selective pressure in the Dendranthema species [14]; and the S-phase kinase-associated protein 1-like (Skp1-like) and F-box families also show a similar pattern of adaptive evolution in A. thaliana [57, 58].

Moreover, paralogs had faster evolutionary rates than orthologs with Ka/Ks ratios greater than 1 among the young duplicate gene families, which could be inferred based on the fact that the former might undergo a higher frequency of adaptive evolution than the latter in the six Rosaceae species [59]. Instances contain genes involved in response to ecological conditions of different stresses under adaptive evolution among the species-specific duplicate families, including LEA genes in M. x domestica, FAR-RED IMPAIRED RESPONSE1 (FAR1) genes in P. communis, and the glucuronyltransferase family (UDPGT domain, UGT) in R. chinensis (Additional file 3: Table S3). In our study, species-specific expansions involved more gene members and gene families compared with lineage-specific expansions among the six Rosaceae species. The results demonstrated that the species-specific duplications mainly contributed to the recent expansions of the six species. It was previously reported that species-specific duplications could result in species-specific functions and traits in a given species, helping the species to maintain better growth and development and exploit stronger adaptability to the constantly changing environment [16, 60].

In apple trees, a period of exposure to low temperature was necessary to enable dormancy induction and growth resumption [61]. The LEA genes encode late embryogenesis abundant proteins in higher plants, many of which are seen to play a role in environmental stress response, such as drought, high salinity, and low temperature [62,63,64]. For example, over-expression of the OsLEA3–2 gene can enhance the salt and drought tolerance of transgenic Arabidopsis and rice during growth periods [65]. In woody plants, including apples, peaches, and grapes, the LEA genes are the major ones among the consistently up-regulated genes induced by the conditions of photoperiod and low temperature [66]. Therefore, it may be inferred that the apple-specific duplication of LEA genes under positive selection might be closely associated with adaption to the physiological characteristic of dormancy.

Species-specific duplications of FAR1 were detected with the feature of adaptive evolution in the pear genome, which are photophilous fruit trees requiring approximately 1600–1700 h of sunlight annually. FAR1 and FAR-RED ELONGATED HYPOCOTYL3 (FHY3) are light-signaling transcription factors, which regulate light-induced inositol biosynthesis, oxidative stress response, starch synthesis, and starch granule formation [67, 68]. The FAR1 and FHY3 proteins promote the tolerance to oxidative stress and prevent the salicylic-acid-dependent cell death by activating the myo-inositol-1-phosphate synthase1 gene (MIPS1) to adapt to the light changes in Arabidopsis [67]; they also modulate light-induced starch synthesis by transcriptional activation of the starch-debranching enzyme ISOAMYLASE2 (ISA2) in response to light and sugars in A. thaliana [68]. Thus, we speculated that pears generate more FAR1 genes to adapt to a given light environment for more sugar synthesis and accumulation for fruit ripening. Besides, half of the FAR1 genes reveal up-regulation during salinity stress treatments in Aegilops tauschii, suggesting that they may play certain roles in salt stress response [69]. Pears have a relatively strong adaptability to ecological conditions, containing the tolerance to cold, drought, and salt stresses [70], indicating that the evolution of FAR1 genes might support the species to adapt to different ecological conditions.

The UGT genes comprise the largest gene family in plants, responding to a variety of biotic and abiotic stresses by glycosylating various phytohormones and other metabolites [71, 72]. UGT73B3 and UGT73B5 in tomato and AtSGT1 and UGT76B1 in Arabidopsis are involved in the biotic stress response to the infection caused by Pseudomonas syringae with the crosstalk of salicylic acid and jasmonic acid [73,74,75]. Similarly, the tobacco glycosyltransferase gene, TOGT1, plays an important role in enhancing the resistance to tobacco mosaic virus (TMV) [76]. In Arabidopsis, UGT genes, UGT85U1/2 and UGT85V1, in response to abiotic stresses, are involved in improving the tolerance to salt and oxidative stress [77]. Therefore, all the above-mentioned examples suggest that duplications and adaptive evolution of UGT genes in Chinese rose (R. chinensis) might be driven by specific ecological conditions, including biotic or abiotic factors.

Conclusions

In summary, this study uncovered the evolutionary patterns of young duplicate genes among the six Rosaceae species. All results showed that functional bias of the young duplicate genes was related to environmental stimuli response. The young duplicate genes were driven by adaptive evolution to adapt to the ecological factors in the six Rosaceae species. This work provides a systematic ways to study young duplications at genome-wide levels among closely related species.

Methods

Identification of young duplicate gene families

The whole genome sequences and annotations of the six Rosaceae species, F. vesca (v4.0.a1), M. x domestica (v1.0), P. communis (v1.0), P. persica (v2.0.a1), R. chinensis (v1.0) and R. occidentalis (v1.0), were downloaded from the Genome Database for Rosaceae (GDR, https://www.rosaceae.org/) [20,21,22,23,24,25, 78,79,80]. An all-vs.-all BLASTN search was performed among the nucleotide sequences (CDSs) of the six species (e-value< 1.0e-40) by using the local BLAST+ program. Based on the BLAST results, a PERL script was performed to identify gene families among the six species, for which the genes with coverage exceeded 60% were considered to constitute a gene family. Following this, these gene families that had at least two members from the same one species and those with the highest identity more than 90% were defined as young duplicate gene families by a PERL script [81].

Two types of expansions, species-specific and lineage-specific expansions, were detected across all young duplicate gene families. Species-specific expansion was characterized by young duplications occurring in only one species; that is, the corresponding species comprised ≥2 members and the other species comprised ≤1 member in each family. Lineage-specific expansion was characterized by the occurrence of young duplications in two or more species. That is, when the lineage-specific expansion was found in < 6 species (2 or 3 or 4 or 5), the involved species had ≥2 members and the other species had ≤1 member in each family; when the lineage-specific expansion was discovered in six species, all the studied species owned ≥2 members in each gene family.

Identification of paralogs and orthologs

In each young duplicate gene family, their proteins sequences were used to run the script with default parameters by using OrthoFinder, and the members were classified into different orthogroups [82, 83]. In each orthogroup, genes derived from the same one species are defined as paralogs, and genes from different species are orthologs based on the illustration of OrthoFinder [82,83,84]. Taken a young duplicate family as an example, the family10 exhibiting the species-specific expansion occurred in F. vesca genome. There were four gene members in this family, two of them from F. vesca (Fv4g31010 and Fv4g31020), one from R. chinensis (RC7G0418500) and the other one from R. occidentalis (T00039). The four gene members were clustered in an orthogroup according to the Orthofinder identification. Therefore, the two genes from F. vesca were considered as paralogs; and the genes from different species, such as Fv4g31010 vs. RC7G0418500, were defined as orthologs. Similarly, the young duplicate families from the lineage-specific expansions also possessed the paralog and ortholog relationships between their gene members.

Phylogenetic tree of young duplicate genes

For the accuracy of sequence alignments, CDSs of the seven young duplicate gene families involved in the six species were firstly translated into proteins, which experienced the alignments by using the MUSCLE program with default options in MEGA7 [85]. Then, the protein alignments were used to guide the alignments of nucleotide sequences, which were employed to constructed the phylogenetic trees using the neighbor-joining (NJ) method with 1000 replicates based on the p-distance model in MEGA7 [85].

Duplication types of young duplicate genes

In each of the six species, the whole-genome protein sequences were considered as both database and query sequences to proceed a BLASTP search by the local BLAST+ program (e-value< 1.0e-40). Moreover, chromosomal locations of whole-genome genes were achieved from the annotation files of the six species. Based on the BLAST results and the physical positions, different modes of duplication gene pairs were classified by using the comparative genomic tool DupGen_finder [86]. Afterwards, the duplication types of young duplicate genes were screened out from the genome-wide classification results in each species.

Domain preference of young duplicate genes

The protein sequences of all young duplicate genes were obtained by a PERL script, used to study the structural domains of them by the local Pfam database with e-value 1.0 (http://pfam.xfam.org/). Then, the domain numbers of the young duplicate genes were counted by Excel according to the Pfam results.

Calculation of nonsynonymous and synonymous values of young duplicate genes

Likewise, the CDS alignments of each young duplicate orthogroup were got from the guidance of protein alignments using the ClustalW2.1 program [87]. The nonsynonymous substitutions (Ka), synonymous substitutions (Ks), and the ratio of the nonsynonymous to synonymous substitutions (Ka/Ks) were calculated using MEGA7 based on the CDS alignments in each orthogroup [85]. For Ks ratios, the values less than 1 were retained considering the saturation of nucleotide substitutions. Afterwards, the counting statistics of Ks values were performed in units of 0.1 in paralogs and orthologs among the two types of expansions by Excel, respectively. Finally, the bar graphs of Ks values were generated by Origin. For Ka/Ks ratios, the values between paralogs and orthologs were sorted by Excel and the graphs were also drawn by Origin.

Calculation of nucleotide diversity values of young duplicate genes

Based on the alignments of young duplicates in each orthogroup, the nucleotide diversity values (Pi) were calculated by MEGA7 [85]. Afterwards, the Pi values of paralogs and orthologs were sorted by Excel.

Chromosomal location of young duplicate genes

The physical positions of young duplicate genes and whole-genome genes were available from the annotation information of the six species. In all species, each chromosome of them was divided into a quantity of windows (1 Mb/1 Window), and the gene numbers and gene densities of species-specific and lineage-specific expansions were calculated in windows on corresponding chromosome of each species, respectively.

Syntenic analysis of young duplicate genes between species

The syntenic relationships between the genes of the six species were detected at genome-wide levels by using TBtools [88] and MCScanX [89]. After obtaining the syntenic blocks between the studied species, the related young duplicate genes from the lineage-specific expansions were picked out to perform the visualization of microsyntenic relationships.