Comparative genomics and transcriptomics analysis reveals evolution patterns of selection in the Salix phylogeny

Transcript sequences of 10 Salicaceae species

SSRs identified in 10 Salicaceae species

Orthologue identification and functional characterization between 10 Salicaceae species

All of the pairwise orthologues were identified by comparative analysis between the 10 Salicaceae species (Table 4). The results showed that S. purpurea had the maximum average number (6597) of orthologous genes, whereas S. babylonica had the minimum average number (707). The highest number of orthologous genes (9713) was found between S. purpurea and S. suchowensis, while the lowest number (681) was found between S. babylonica and S. fargesii. 238 single copy orthologues were found in all 10 Salicaceae species (Fig. 2). The orthologues were annotated with GO terms (Additional file 1: Table S3). Taking P. trichocarpa as an out-group species, the phylogenetic tree of Salix was constructed based on combined 238 orthologous transcripts using Maximum Likelihood (ML) method (Fig. 4).

	S. purpurea	S. suchowensis	S. sachalinensis	S. dasyclados	S. viminalis	S. eriocephala	S. fargesii	S. matsudana	S. babylonica
S. purpurea
S. suchowensis	9713/0.02
S. sachalinensis	6319/0.02	5035/0.02
S. dasyclados	6598/0.02	5221/0.02	5147/0.02
S. viminalis	6898/0.03	5429/0.03	4812/0.03	4950/0.02
S. eriocephala	6585/0.03	5238/0.03	5132/0.03	5301/0.02	4923/0.02
S. fargesii	6258/0.03	4911/0.03	4669/0.02	4765/0.02	4638/0.02	4756/0.02
S. matsudana	7934/0.04	6251/0.04	5217/0.04	5427/0.03	5322/0.04	5411/0.04	5078/0.04
S. babylonica	892/0.05	696/0.05	708/0.04	714/0.04	688/0.03	713/0.04	681/0.03	726/0.02
P. trichocarpa	8179/0.11	6246/0.12	4165/0.1	4325/0.11	4519/0.11	4400/0.11	4107/0.1	5317/0.11	549/0.11

Phylogenetic analysis and divergence time

The genetic distance of species was related to synonymous mutation rate calculated by orthologous genes, so the synonymous mutation rates of all pairs of orthologues were estimated in 10 Salicaceae species (Table 4). Between different branches (Fig. 3), S. purpurea has the Ks peak (0.02) with S. suchowensis, S. sachalinensis and S. dasyclados, 0.03 Ks peak with S. viminalis, S. eriocephala and S. fargesii, 0.04 Ks peak with S. matsudana, 0.05 Ks peak with S. babylonica, and the maximum Ks peak 0.11 with out-group P. trichocarpa. Between different genera, most of Salix species has the Ks peak 0.11, whereas S. fargesii was found the minimum Ks peak 0.10 with P. trichocarpa (Table 4). It is suggested S. fargesii was a relatively ancient species compared to others.

Using P. trichocarpa as an out-group species, the phylogenetic tree of Salix was derived with the pairwise Ks values of the orthologous transcripts as a distance metric based on the neighbour-joining (NJ) method (Fig. 4). In the phylogenetic tree, the average Ks value is 0.11 between Genus Salix and Genus Salix (Calculated by Table 4), and which is nearly consistent with the value of 0.12 in previous studies [18]. Based on existing fossil evidence, the divergence time of genera Salix and Populus was about 48 million years old in middle Eocene sediments [16, 17]. With this time as the separation of the two lineages and K = 0.11, the rate of substitution (r) was calculated to about 1.14 * 10^− 9 per site and year (T = K/2r), and which is very close to previous value of 1.28 * 10^− 9 [18].

Using the fossil calibrations (48 Mya) of genera Salix and Populus [16, 17], the divergence times were estimated based on the 238 single copy genes and pairwise Ks distance metric of the orthologous transcripts (Table 4). The divergence of subgenus Vetrix and Salix occurred at about 17.6–16.0 Mya in the Salix phylogeny, and S. fargesii diverged at about 10.9–10.6 Mya with other species of subgenus Vetrix (Fig. 4). There were still some inconsistencies on the divergence time of subgenus Vetrix and Salix based on nuclear and plastome genes in previous studies [13, 15]. The time of 17.6–16.0 Mya between subgenus Vetrix and Salix supports the value of 16.9 Mya estimated by complete plastome genomes.

Evolutionary pattern of Salix spp. genes

Ka/Ks rate of orthologous genes could reflect the evolution pattern of species. Ka/Ks > 1 indicates that the gene has involved in positive selection during evolution.

In subgenus Vetrix except S. fargesii (Table 5), S. dasyclados was identified 454, 306, 267 and 289 positive selection genes with the spepies of subgenus Vetrix, S. purpurea, S. suchowensis, S. sachalinensis and S. eriocephala. Between them, cold- stress genes were found to be annotated to NP_190879.1, AAM23265.1, NP_849749.1 and AAN77157.1 (Additional file 1: Table S4), which producing the proteins of P-loop NTPases [35], L-asparaginase [36], HOS10 with Myb damain [37] and thylakoid-bound ascorbate peroxidase [38]. 254 positive selection genes were identified between S. sachalinensis and S. viminalis, and one light-stress gene (NP_565524.1) was found by producing the SEP protein [39].

In subgenus Salix and S. fargesii (Table 5), 257 and 36 positive selection genes were identified between S. matsudana and S. fargesii, S. matsudana and S. babylonica. Universal-stress genes (NP_001132550.1 and NP_001132238.1) were wildly found between them by producing universal stress protein. Universal stress protein could induce by many environmental stressors such as nutrient starvation, drought, extreme temperatures, high salinity, and the presence of uncouplers, antibiotics and metals [40–42].

Salix phylogeny derived by comparative genomics and transcriptomics

In previous studies, Salix phylogenies were usually derived by several nuclear and plastid markers [7–15]. Different markers always obtained different phylogenetic trees. Comparative genomics and transcriptomics could make use of more and more nuclear sequences. Phylograms were derived using two methods in this work. 238 single copy genes were strictly selected to construct the phylogenetic tree by maximum- likelihood method, which used most of the sites. Another method is based on the neighbor-joining method using the pairwise Ks values of the orthologous transcripts as a distance metric, which used most of the orthologous. The divergence times of subgenus Vetrix and Salix estimated by two methods were consistent with the value by complete plastome genomes [15]. It is improved that enough single copy nuclear sequences should obtain similar results with enough plastid sequences in the Salix phylogeny.

Paleoclimate change in the divergence of Salix phylogeny

The divergence time of genera Salix and Populus was about 48 Mya at the period of Paleogene (66-23Mya) [43–45]. During the Paleogene, the global climate went against the hot and humid conditions of the late Mesozoic era and began a cooling and drying trend [23]. As the Earth cooled, tropical plants were restricted to equatorial regions and became less numerous. Deciduous plants became more common which could survive through the seasonal climates, during which Salix and Populus diverged.

Miocene (23 - 5Mya) is the main period in the divergence of Salix phylogeny (Fig. 4). During the period, there is evidence of a warm period from 21 Mya to 14 Mya named as the Middle Miocene Climate Transition (MMCT) [23], and the rare pleasant environment might cause the species diversity. The divergence time of subgenus Vetrix and Salix was about 17.6–16.0 Mya corresponding to the period of MMCT. Then global temperatures took a drop and some species were extinct by 14Mya [46–48], so the north subgenus Vetrix needed to migrate or adapt in order to survive. One group with S. fargesii diverged from subgenus Vetrix and migrated to south. The resident group of Vetrix had to adapt the cold and drought climate. By 8 Mya, the climate sharply cooled and formed the Quaternary Ice Age (2.6–0.1Mya) [49]. The climate change from MMCT to Quaternary Ice Age should play an important role in the divergence of Salix phylogeny.

Universal- stress genes and migration of S. fargesii

The divergence time of S. fargesii (section Psilostigmatae) was about 10.9–10.6 Mya after the MMCT (21–14 Mya), and in which period the climate changed to be cooling. Wind and animal pollination had been proved to play an important role in the spread of willows [50–52]. Section Psilostigmatae are mainly distributed along the Changjiang (Yangtze) [53], Laantsang and Nujiang river of China (Fig. 1). The river provided the feasibility for the migration of the willow catkins by animal or other pollination, which is consistent with the distribution of Section Psilostigmatae. In previous studies, it was shown that the evolutionary history of the salix has involved multiple reticulation events that may mainly be due to hybridization [13]. Migration of S. fargesii maybe provided the possibility for the hybridization of Salix.

Our result shows that universal- stress genes were identified to be involved in positive selection between S. fargesii and subgenus Salix. It is suggested that selective evolution of universal- stress gene should play an import role in migrating to south for S. fargesii. When S. fargesii moved to south, universal- stress gene can help to produce the abiotic resistances to adapt the complex environment of south areas.

Cold-, light- stress genes and the north resident group of sub genus Vetrix

After the MMCT (21–14 Mya), the global climate came back the cooling and drying trend [23]. In previous studies, it was shown that the evolutionary history of the Salix maybe affected by the profound climatic cooling during the Tertiary [13]. The resident group of section Vetrix had to adapt to the cooling stage especially the high-latitude. Cold-, light- stress genes were widely identified to be involved to positive evolution among S. purpurea, S. suchowensis, S. sachalinensis, S. eriocephala and S. dasyclados. It is suggested that the cool and dry climate had played an important role in the speciation of north group of Section Vetrix.

Data sources

In order to discover the evolutionary pattern of orthologues, the cDNAs and transcripts of 9 Salix and one Populus (out-group) were downloaded from the public databases (Table 1). The cDNAs of P. trichocarpa (v3.1), Salix purpurea (v1.0) and S. suchowensis were directly derived from the JGI [5] and willow genome project of NJFU ( http://bio.njfu.edu.cn/ss_wrky ). Transcriptome sequencing of S. sachalinensis, S. dasyclados, S. viminalis, S. eriocephala, S. matsudana, S. babylonica and S. fargesii were obtained from the SRA database of NCBI. Geographic distributions of section Psilostigmatae was draw by ArcGIS based on Flora of China ( http://frps.iplant.cn ) (Additional file 1: Table S2).

Data filtering and de novo assembly

SRA datasets with FASTQ format were filtered to remove raw reads of low quality. Transcriptome assembly was achieved using the short-read assembly program Trinity [54]. The assembled sequences (> = 300 bp) were combined and clustered with CD-HIT (version 4.0) [55, 56]. Sequences with similarity > 95% were divided into one class, and the longest sequence of each class was treated as a unigene during later processing.

Identification of SSRs in 10 Salicaceae species

Putative SSRs in Unigenes and cDNAs were identified by MISA software. The options of Di- to hexa-nucleotide SSRs were set to 6 (for di-), 5 (for tri- and tetra-), and 4 (for penta- and hexa-), and all SSRs were characterized in 10 Salicaceae species.

Identification of orthologues among 10 Salicaceae species

OrthoMCL software [57] was used to cluster the transcribed sequences. Based on the proteins of Salix purpurea as reference, one-to-one sequences of each group were then filtered to use in subsequent analyses. The annotations obtained from Nr were processed through the BLAST2GO program [58] to get the relevant GO terms. Heatmap of orthologues was draw by R language.

Estimation of synonymous substitution and non-synonymous substitution rates

In order to remove the unigenes without open reading frames, pair-wise orthologues were searched against plant protein sequences of GenBank with BLASTX tool. The method has been used in previous studies [59]. Clustalw software [60] was used to align the filtered pair-wise orthologues, and the output files were formatted to NUC format for subsequent analysis. The rates of synonymous substitutions (Ka) and non-synonymous substitutions (Ks) were estimated with PAML software [61].

Phylogenetic analysis

There were still some inconsistencies on phylogenetic relationship in previous studies. Phylograms were derived using two methods in this work. Single copy genes by orthoMCL were aligned by Muscle [62] and formated by Gblock [63], maximum- likelihood method was used to build the phylogenetic tree by MEGA6 [64] (bootstrap is 1000 and Kimura 2-parameter model). Another method is based on the neighbor-joining method of MEGA6 [64] using the pairwise Ks values of the orthologous transcripts as a distance metric (Table 4). Populus trichocarpa was used as an out-group to root trees.