Three species, S. tetrasperma, S. babylonica, and S. oreinoma, representing two subgenera of the genus Salix, were sampled. We collected healthy, tender and fresh leaves from adult plants of target species. The voucher herbarium specimens for the three sampled Salix species are all deposited in Herbarium of Kunming Institute of Botany, Chinese Academy of Sciences (KUN) (Supplementary Table S1).

Total DNA enriched with cp DNA was extracted from 200 g of fresh leaves according to the methods of Jansen et al. (2005) and Zhang et al. (2011). Purified DNA (5 mg) was fragmented and used to construct short-insert libraries according to the manufacturer’s manual (Illumina Inc., San Diego, CA, United States). DNA from the different individuals was indexed by tags and pooled together in one lane of Illumina’s Genome Analyzer for sequencing.

We filtered out non-cp DNA reads from the raw sequences based on the known cp genome sequences. Next, the filtered reads were used to de novo assemble the cp genomes with SOAPdenovo software, which is specially designed to assemble Illumina short reads. SOAPdenovo pipeline (e.g., k = 31 bp and scaffolding contigs with a minimum size of 100 bp) can carry out accurate analyses of unexplored genomes, resolve repeat regions in contig assembly and improve gap closing, etc., in a cost effective way (Luo et al., 2012). Then, all contigs were mapped to the reference cp genome of P. trichocarpa using BLAST¹ search from NCBI with default parameters. The orders of aligned contigs were determined according to the reference genome. Gaps between the de novo contigs were replaced with consensus sequences of raw reads mapped to the reference genomes.

Based on the reference genomes, the four junctions between LSC/IRs and SSC/IRs were confirmed with PCR-based product sequencing, respectively. To avoid assembly errors and to obtain high quality complete cp genome sequences, validation of assembly was also carried out with intensive PCR-based sequencing. We designed 39 pairs of primers (see Supplementary Table S2 for detail) based on the varied regions of the eight preliminary cp genome assemblies. PCR products were sequenced using the BigDye v3.1 Terminator Kit for ABI 3730xl (Life Technologies, Carlsbad, CA, United States). Sanger sequences and assembled genomes were aligned using Geneious 7 (Kearse et al., 2012) to determine if there were any differences. The complete cp genome sequences were deposited in GenBank (Supplementary Table S1).

Genome annotation was accomplished using the Dual Organellar Genome Annotator (DOGMA) (Wyman et al., 2004) and also compared with the available complete chloroplast genome of P. trichocarpa (GenBank accession number NC009143) to annotate the genes encoding proteins, transfer RNAs (tRNAs), and ribosomal RNAs (rRNAs). All of the identified tRNA genes were further verified using the corresponding structures predicted by tRNAscan-Se v1.21 (Schattner et al., 2005).

We performed phylogenetic analyses based on different datasets, i.e., whole cp genome minus the second inverted repeat region (IRa) to avoid considering the same information twice, concatenated non-coding sequences (including intergenic spacers and introns), and concatenated exons of protein-coding genes (ndhA was excluded due to its pseudogenization in S. babylonica). Sequences were aligned with Muscle [62]. After manual correction in Geneious (Kearse et al., 2012), phylogenetic analyses were performed based on maximum likelihood (ML) criteria and Bayesian inference (BI). The ML analysis employing a GTR+G model of substitution for all datasets was run in RAxML v7.2.8 (Stamatakis, 2006), and a bootstrap analysis of 2000 replicates was performed simultaneously (option “-f a”). BI was performed with MrBayes v3.2 (Ronquist et al., 2012). Two independent Markov chain Monte Carlo (MCMC) chains were run, each with three heated and one cold chain for three million generations. Each chain started with a random tree, default priors, and sampling trees every 100 generations. The GTR+G model was used for all datasets as suggested by RaxML manual.

Topological incongruence among conflicting datasets was tested using Dendroscope v. 3.2.8 under the galled network consensus algorithm (Huson and Scornavacca, 2012). ML trees of concatenated protein coding sequences and whole genome sequence were used to infer phylogenetic networks.

We collected the coding DNA sequence (CDS) of each orthologous gene in eight Populus species and six Salix species, and aligned them with translation aligner using Geneious (Kearse et al., 2012), which could take into account frame shifts and premature stop codons, and generate codon-based alignments. We used the branch site model of PAML to compute the Ka/Ks ratio for orthologous genes in each external and internal branch of phylogeny trees that were generated based on protein coding sequence alignment with ML method. We tested two branch site models (with the parameters model = 2 and NSsites = 2): the “model 1” with both the branch site specific Ka/Ks and background Ka/Ks varying freely, and the “model 2” with the branch site specific Ka/Ks fixed at 1 and background Ka/Ks varying freely (Yang, 2007). We then performed the Likelihood Ratio Test (LRT), which tests whether the likelihood of the “model 1” is significantly different from that of the “model 2” by comparing two times the log likelihood difference. We computed p-values using a chi-square distribution with one degree of freedom (Yang, 1998).

We used the pipelines developed by UCSC genome browser to search for the homologous regions of the chloroplast genome in the nuclear genome for P. trichocarpa² and S. purpurea³, respectively, since the nuclear genome sequences are only available for these two species (Tuskan et al., 2006). We first aligned the chloroplast genome to the nuclear genome of the two reference genomes with LASTZ (Harris, 2007). We then transformed the “lav” output format of LASTZ to “axt” format using lavToAxt. Finally we chained the “axt” files using axtChain and generated chain format outputs (Kent et al., 2003; Schwartz et al., 2003).

Although we did not set the identity filter of LASTZ for blocks or HSPs (high-scoring segment pairs), we computed the identity of the final homologous sequence pairs generated by axtChain. The identity ranges from 60 to 99%. Based on the chain file, we generated the homologous regions between the chloroplast genome and nuclear genome, which could represent chloroplast-nuclear DNA transfer events.

To test the hypotheses of whole chloroplast genome horizontally transferring to nucleus in P. trichocarpa, we conducted 100,000 Monte Carlo sampling test, which is implemented by the comparison of the observed data with random samples generated in accordance with the hypothesis being tested. The significance of the Monte Carlo sampling test is determined by the rank of the test criterion of the observed data relative to the test criteria of the random samples composing the reference set (Hope, 1968). Specifically, we assumed the whole chloroplast genome was split into eight fragments, which were inferred from the alignment between the chloroplast genome and nuclear genome in species P. trichocarpa. For each simulation, we randomly sampled the insertion locations of the eight fragments in the whole P. trichocarpa genome. We then counted the number of simulations, in which the five fragments (i.g. the 1st, 3rd, 5th, 7th, and 8th fragments) were inserted into the same chromosome within the range of the length of chromosome 13 and were arranged according to their previous order in the chloroplast genome. Next, we divided this number by the number of simulations (i.g. 100,000). This mathematical derivative is treated as the P-value for inference of the testing hypothesis significance.

Using the Illumina Hiseq 2000 system, we obtain 1,392,310, 3,779,094, and 4,595,286 bp paired-end clean reads (average read length 91 bp) for S. babylonica, S. tetrasperma, and S. oreinoma, respectively. We mapped these sequence reads to the reference cp genome of S. purpurea and achieved at least 1600× (1615× for S. babylonica, 4384× for S. tetrasperma, 5330× for S. oreinoma) coverage for cp genome. Based on de novo and reference-guide assembly, we obtain the complete cp genome for S. tetrasperma. The assemblies of the other two cp genomes contain seven to eight gaps, and we filled the gaps using PCR-based sequencing.

Four junction regions of cp genomes were validated using PCR-based sequencing for each cp genome, respectively. Furthermore, in order to overcome the errors of heterogeneous indels from homopolymeric repeats (Moore et al., 2006; Yang et al., 2010), we conducted PCR-based validation to correct the errors. We designed 18 pairs of primers based on the variable regions of alignments to validate these sequences in each cp genome (Supplementary Table S2). In total, we amplified and sequenced ∼ 152 kb from all the three Salix species. Then, we compared these sequences directly to the assembled genomes and we observed no nucleotide mismatches or indels/insertions. This result confirmed the reliability of assembled chloroplast genome sequences. Finally, we obtain the complete cp genome sequences of S. babylonica, S. tetrasperma, and S. oreinoma.

The complete cp genomes of the three Salix species sequenced vary from 155,531 to 155,740 bp in size and exhibit a typical circular structure including a pair of IRs (range from 27384 to 27436 bp) and two single-copy regions (LSC, 84466–84580 bp; SSC, 15862–16323 bp). Each of the three cp genomes contains 111 unique genes (110 for S. babylonica) (Figure 1 and Table 1), including 77 unique CDSs (76 for S. babylonica because of the pseudogenization of ndhA), four unique rRNAs, and 30 unique tRNAs. Seventeen genes contain introns; 14 of them (atpF, ndhA, ndhB, petB, petD, rpl2, rpl16, rpoC1, trnA^UGC, trnG^GCC, trnI^GAU, trnK^UUU, trnL^UAA, and trnV^UAC) exhibit one intron and three of them (clpP, rps12, and ycf3) contain two introns. All the CDSs have canonical ATG start codons except ndhD, which has GTG as the start codon.

Based on cp genome alignment from 14 Salicaceae species, we identified 30 highly divergent non-coding regions for implementation of phylogenetic analysis in Salix, Populus and Salicaceae. Those molecular markers are mostly derived from intergenic or intronic non-coding regions, though four protein-coding genes including ccsA, rpl20, rps7 and ndhA/ndhE are also identified. Among the four genes, rps7 is highly diverged and lineage-specific within genus Salix. The length of rps7 is 468 bp for all Populus species except P. cathayana, whereas the three Salix species (i.g. S. oreinoma, S. suchowensis, and S. purpurea), which are closely related to subgenera Chamaetia and Vetrix, have identical rps7 sequences with a length of 270 bp. The other three Salix species (i.g. S. interior, S. babylonica, and S. tetrasperma) in subgenus Salix have identical rps7 sequences with a length of only 180 bp (Figures 2A, 3, and Supplementary Table S3).

The phylogenetic resolution of the genus Salix is proven to be extremely difficult (Skvortsov, 1999; Argus et al., 2010; Wu et al., 2015) and no agreement on Salix phylogeny is made thus far. Molecular phylogenetic analysis in Salix only succeeds in resolving relationships within the subgeneric level, but fails in the taxonomic level under subgenus and shows almost no resolution in the clade with about 73% of Salix species (i.e., subgenera Chamaetia and Vetrix) (Leskinen and Alstrom-Rapaport, 1999; Azuma et al., 2000; Chen et al., 2010; Hardig et al., 2010; Abdollahzadeh et al., 2011; Wu et al., 2015). Therefore, there is a need to identify highly divergent regions in Salix cp genomes as molecular markers.

The genus Populus is comprised of ca. 30 species, which can be divided into six major clades (Eckenwalder, 1996; Argus et al., 2010). We selected eight species from four clades in our phylogenomic analyses. Our analysis indicates that Populus is monophyletic. Phylogeny, based on whole genome and concatenated non-coding datasets, shares the same topology (Figure 2B). However, phylogenetic relationship based on concatenated protein coding dataset is different in the topological positions of P. euphratica, P. cathayana and P. balsamifera (Figure 2A), as indicated by network analysis (Figure 2C). The three New World Populus species (P. trichocarpa, P. balsamifera, P. fremontii, Clade A in Figure 2A) form a robust monophyly. This clade is sister to the remaining Populus species in the protein coding gene tree (Figure 2A). However, P. euphratica is sister to P. cathayana from the tree based on whole cp genome and non-coding regions (Figure 2B). In the tree based on protein coding regions, P. euphratica and P. cathayana fall in the clade that contains all species from Old World with low support value and short branch length (clade B in Figure 2A). This topology is absent from previous studies with more Populus species sampled (Hamzeh and Dayanandan, 2004; Wang et al., 2014), indicating that insufficient sampling might have caused this topology.

We identified at least four possible pseudogenization events in Salicaceae: (1) ndhA in S. babylonica, which contains a 1039 bp deletion including a start codon in exon 1; (2) accD in P. yunnanensis, which loses 1056 bp (72%) of gene content; (3) ycf1 in P. yunnanensis, which is deleted >5000 bp (94%) of DNA sequences; (4) psbC in P. cathayana, which misses ∼1000 bp (∼70%) of DNA sequences (Figure 2A).

We identified seven significantly positively selected genes (atpE, rps7, ycf2, ccsA, petD, psbC, and psbJ) that contain positively selected sites (Table 2). Among them, three positively selected genes (i.e., rps7, petD, and psbC) are identified in P. cathayana, ycf2 is identified in P. yunnanensis, ccsA and psbJ are identified in P. tremula. An exception, atpE, shows significant positive selection on one site in the Salix subgenus Salix clade (i.e., S. babylonica, S. tetrasperma, and S. interior) (Table 2).

We used LASTZ (Harris, 2007) to search for NUPT including DNA fragments less than 200 bp. We identified 571 and 713 NUPTs with total ungapped lengths of ca. 536 and 193 kb in P. trichocarpa and S. purpurea, respectively (Figure 4, Table 3, and Supplementary Table S7). This result is different from a previous study that used BLAST method for searching NUPT (Yoshida et al., 2013). The number of NUPTs in S. purpurea is larger than that in P. trichocarpa, but the total length of NUPTs in S. purpurea is much lower than that in P. trichocarpa (Table 3).

The fragmented assembly of S. purpurea genome (2015) might cause uncertainty for NUPT searching. Therefore, we subsequently discuss the plastid-to-nucleus DNA transfer in Salicaceae, based on the analysis from P. trichocarpa. The most abundant NUPTs in P. trichocarpa are around 70–199 bp (Figure 4). The number of NUPTs with 1000 bp longer is 77, and these large NUPTs are responsible for the majority of the total NUPT length (>75%, ca. 410 kb out of 536 kb). Furthermore, the average identity of those large NUPTs is 90.2%, and that of the six largest NUPTs with length greater than 10 kb is 99.4% (Supplementary Table S7), indicating that large NUPTs are relatively recently integrated into the nucleus (see Data Sheet S1 for detail). Additionally, all large NUPTs contain large gaps.

In our analyses, we identified five aligned regions (∼4,000–16,000 bp in length) and three/four gaps (∼7,000–15,000 bp in length) between the chloroplast and nuclear genomes in the same order as in the cp genome by LASTZ (Figure 5). We conducted MC sampling to demonstrate that it is more likely due to the insertion of the whole cp genome into the nuclear genome and generation of the gaps through insertion/translocation rather than the separate insertions of split cp genome fragments. Previously, it is also discovered that large cp fragments can be inserted into the nuclear genome based on the similarity between cp genome sequence and nuclear genome sequence using the sequence alignment tools (e.g., BLAST), and sequencing the insertion junction product (Yuan et al., 2002; Yoshida et al., 2013). Here, we used LASTZ as the alignment tool, which is a pairwise aligner for aligning DNA sequences and is originally designed to align the sequences in the size of human chromosomes and from different species (Harris, 2007). Therefore, it is more suitable to detect similarity in a large scale, such as between the cp genome and the nuclear genome. After alignment, we found that three large chunks of cp genomic regions are highly similar to three regions in nuclear genome, and more importantly, these regions are in the same order between the cp genome and nuclear genome.

Therefore, we conducted MC sampling to test the hypotheses of “whole cp genome insertion” or “split cp genome insertion,” which has not been done for detecting a whole cp genome transfer event before. We determine that the probability of those gaps produced by “split cp genome insertion” is less than 1e-05 and is significantly low by chance. Consequently, our simulation analysis supports that the whole chloroplast was transferred to the nuclear genome first and the gaps were created by subsequent insertion/translocation. Lastly, we found that the nucleus genes of NUPTs in P. trichocarpa, each of which matches multiple cp genes, contain multiple cp genes mostly in their original order in the cp genome (Supplementary Table S8). Taken together, all observations revealed by our analyses clearly favor the “bulk” and “insertion” DNA hypothesis of organelle-to-nucleus DNA horizontal transfer. Our approaches can serve as an example to analyze large-scale transfer/insertion events in other species.

We compared the three sequenced Salix cp genomes in this study with previous published complete cp genomes of 11 Salicaceae species, including eight Populus and three Salix species (Supplementary Table S2). We found that the structure and synteny of the cp genomes of the 14 Salicaceae species are highly conserved (Supplementary Figure S1). Also, the lengths of various segments, or parts (IRs, LSC, SSC, coding sequence and non-coding sequence), of the cp genomes are also quite conserved and vary in a small range (Supplementary Figure S2). We observed that four genes (i.g. infA, rps16, rpl32, and ycf68) were lost from Salicaceae cp genomes compared with other angiosperms, among which, infA and rps32 were reported to be transferred into the nucleus (Millen et al., 2001; Ueda et al., 2007). The inverted repeat of the Salicaceae cp genome results in the complete duplication of rps19, rpl2, rpl23, ycf2, ycf15, ndhB and rps7, as well as exons 1 and 2 of rps12, all four rRNA genes (4.5S, 5S, 16S, and 23S) and seven tRNA genes (trnI^CAU, trnL^CAA, trnV^GAC, trnI^GAU, trnA^UGC, trnR^ACG, and trnN^GUU). Both IRs of the Salicaceae cp genome run 1,695–1,907 bp into ycf1, which differs largely in genera Salix and Populus: 1,747 bp in all Salix species and 1,695–1,907 bp in Populus species, which result in most proportions of fluctuations in IR length. Additionally, IRb runs 50–51 bp into rpl22 (Supplementary Figure S3).

Similar to other angiosperms, the IR (pairwise identity 98.1% for both IRa and IRb) regions are more conserved in the nine species than the LSC (pairwise identity 94.0%) and SSC (pairwise identity 94.9%) regions. Differences are observed in the cp genome in the 14 Salicaceae species including genome size, gene losses, the pseudogenization of protein-coding genes, and IR expansion and contraction (Table 1 and Supplementary Figure S2).

Phylogenetic analysis, based on the complete cp genomes and non-coding and protein-coding datasets, indicates that the phylogeny of Salicaceae s.str, Populus and Salix s.l. are all resolved as monophyly, which is largely consistent with previous studies (Davis et al., 2005; Wurdack and Davis, 2009; Chen et al., 2010; Wang et al., 2014; Wu et al., 2015) (Figure 2 and Supplementary Table S4). However, we identified the incongruence of phylogenies in Salicaceae that are inferred from these three datasets (i.e., whole genome and concatenated coding sequence between non-coding sequences of the cp genome). This incongruence might be caused by homoplasy, which could result from convergence, parallelism or reversal. Because the cp genome is inherited maternally as a single unit, the observed phylogenetic tree incongruences in our result are most unlikely caused by lineage sorting or hybridization/introgression, which are generally used to explain the conflict signals among characters of plant taxa including Salicaceae (Wang et al., 2014; Wu et al., 2015).

Our phylogenetic tree, based on the whole cp genome, coding-sequence and non-coding sequences of cp genome, shows that all conflicting branches are short (Figures 2A,B), which could be caused by fast species radiation or short stem-lineages (Wagele and Mayer, 2007). In this case, apomorphies evolved in stem-lineage may be rare, and chance similarities (non-phylogenetic signal produced by nucleotide substitution processes) that evolved later can accumulate and dominate in the form of phylogenetic-signal-like (i.e., false and misleading phylogenetic signals) patterns. It is difficult to distinguish these kinds of homoplasies from apomorphies, and consequently, this might lead to the wrong phylogenetic tree (Wagele and Mayer, 2007). Moreover, multiple nucleotide substitutions along long branches may destroy synapomorphies, resulting in accumulation of homoplasies along long branches and attracting distantly related clades to be clustered together in a topology (i.e., long-branch attraction) (Wagele and Mayer, 2007). Therefore, the conflict topologies in our phylogenomic trees most likely result from incomplete sampling or homoplasy. This is in line with previous cp DNA markers based on phylogenomic study of bamboo, which comes to the conclusion that homoplasy should be in account for conflicting phylogenetic signals between cp genome datasets (Ma et al., 2014). Moreover, our analyses show that the phylogenetic tree based on genomic data with large number of informative characters (e.g., cp genomes) should be used in caution when the examined taxa are quickly radiated or bear short branches.

Phylogenomic analysis indicates that the species of Salix were clustered in a robust monophyletic clade. Two robust subclades were resolved within Salix. One subclade includes the species of subgenera Chametia and Vetrix. Two species (S. suchowensis and S. purpurea) from subgenus Vetrix cluster as a monophyletic group, which is sister to a group including a species, S. oreinoma from subgenus Chamaetia. The other subclade is comprised of species from subgenus Salix, the New World species, S. interior, which is sister to a group containing two Old World species, S. babylonica and S. tetrasperma (Figure 2). The relationships resolved in this study are in line with previous phylogenetic studies of Salix (Azuma et al., 2000; Chen et al., 2010; Hardig et al., 2010; Wu et al., 2015), but we provide evidence from a phylogenomic perspective.

In the evolutionary process of a certain lineage of an organism, changing environments (e.g., climate changing) impose selective pressures and result in adaptive evolution. Identification of genes involved in this process (i.e., positively selected genes) is central to understanding the evolutionary pattern of organisms (Yang, 1998), and pinpointing specific targets for adaptive studies.

Among seven significantly positively selected genes in Salicaceae cp genomes, gene rps7 shows a frame shift mutation. A single nucleotide insertion at site 102 causes a frame shift mutation for gene rps7 and results in a large variation of an amino acid sequence that is shown under positive selection (Figure 3). This frame shift mutation further introduces an early stop codon and shortens its amino acid sequence about 42% of the length compared to other Populus species. The rps7 gene encodes the ribosome S7 protein, also known as ribosomal protein S7 (uS7), which is crucial for the assembly and stability of the ribosome. S7 protein is an important part of the translation process, and is universally present in the small subunit of prokaryotic and eukaryotic ribosomes, and might play either a general or a specific regulatory role in translation initiation in the chloroplast (Fargo et al., 2001). The rps7 gene of P. cathayana is positively selected, although it contains a premature stop codon caused by a frame shift mutation compared with other Populus species surveyed in this study (Figure 3). Despite the premature stop codon, this truncated rps7 gene in P. cathayana has a conserved domain for the uS7 superfamily, as revealed by NCBI-CDD blast results (Supplementary Figure S4). This reveals that this shortened rps7 gene may still function normally. Alternatively, even if the truncated cp rps7 cannot maintain its original function, it is possible that an intact copy of cp rps7 has been transferred into the nucleus and can properly function in the nucleus genome. P. cathayana does not have complete genome sequences, so we searched the cp rps7 gene in nucleus genome of P. trichocarpa, which has fine quality whole nuclear genome sequences. We observe that the Potri.013G138900.1 gene in chromosome 13 is identical to cp rps7 in length and has 100% identity for both gene and coding sequences. This implies that rps7 may have functionally transferred into the nucleus and the cp copy of rps7 might be freely subject to natural selection.

Similarly, the protein coding sequences of petD and psbC genes in P. cathayana are shortened by about 29 and 70%, respectively. Compared to other Salicaceae species, both mutations are caused by a single nucleotide insertion. The petD gene encodes subunit IV of the cytochrome b₆/f complex. It is required for photosynthetic electron transport and hence, supports photosynthetic growth. The mutations in the 5′ UTR or initiation codon can affect its function (Chen et al., 1993; Sakamoto et al., 1994; Sturm et al., 1994). The psbC gene encodes one of the components of the core complex of photosystem II. It binds chlorophyll and helps catalyze the primary light-induced photochemical processes of PSII (Rochaix et al., 1989; Cai et al., 2010). However, the effects of shortened coding sequences of petD and psbC on their function remains unknown, especially for psbC, which is shortened about 70% of the length. Our analysis finds that petD and psbC have been transferred into the coding region of the nuclear genome for three and eight times in P. trichocarpa, respectively (Supplementary Tables S5, S6).

Three of the last few amino acids of ycf2 in P. yunnanensis are detected to be under significant positive selection. These amino acid changes are caused by a frame shift (a single nucleotide insertion) 6 bp ahead of them. Consequently, ycf2 in P. yunnanensis is 21 bp shorter than in other Salicaceae species. The ycf2 gene is a putative ATPase with unknown function. This gene exists in many plants, including non-photosynthetic plants. Previous experiments in tobacco indicate that it plays an essential role in cell survival in the tobacco chloroplast (Drescher et al., 2000). Our analysis indicates that ycf2 is transferred into the nuclear genome 18 times and all targeted in protein-coding regions in P. trichocarpa (Supplementary Tables S5, S6).

We identified two positively selected genes, ccsA and psbJ in P. tremula. Species of P. tremula mainly distribute in cool temperate regions of Europe and Asia. The protein encoded by ccsA is a component of the cytochrome c synthase complex of the membrane-bound System II, and is required during the biogenesis of c-type cytochromes at the step of heme attachment (Xie et al., 1998). The psbJ encodes one of the components of the core complex of photosystem II. Experiments in tobacco indicate that plants with a mutated psbJ gene are unable to grow photoautotrophically (Hager et al., 2002). Our analysis shows that ccsA has been transferred to the coding region of the nuclear genome once, but psbJ is not transferred in P. trichocarpa (Supplementary Tables S5, S6).

All the positively selected genes locate in terminal species, except atpE, which shows significant positive selection on one site in the Salix subgenus, Salix clade (i.e., S. babylonica, S. tetrasperma, and S. interior) (Table 2). The atpE gene encodes the ε subunit of CF1 of the H⁺-translocating ATP synthase, and functions in part to prevent wasteful ATP hydrolysis by the enzyme (Cruz et al., 1997). The common ancestors of this clade might be adapted to temperate conditions and diverge in about early Oligocene (Wu et al., 2015). The positively selected aptE gene might be related to the adaptive evolution of the common ancestor group of Salix subgenus Salix.

As mentioned above, most of the positively selected protein-coding genes are transferred into the nucleus in P. trichocarpa. Positive selection is generally regarded as evidence of adaptive evolution, and these positive selected genes might have driven the successful adaptation of the selected taxa and/or lineages. However, the specific functions or effects of these positively selected cp genes in the target species remain unknown, and structural biological studies are needed to clarify the implication of these findings.

Most of the chloroplast genes are transferred to nucleus and then deleted from the plastome. However, a transferred cp gene is not readily expressed in the nucleus, nor would it be able to give rise to a product equipped with the capability of returning to the chloroplast and ousting its progenitor chloroplast gene. The NUPT events have happened repeatedly and still are active during endosymbiotic evolution (Timmis et al., 2004; Cullis et al., 2009; Yoshida et al., 2013). Although some previous studies investigated the patterns of genomic integration of NUPTs in plant species, the mechanisms of cp-to-nucleus gene transfer are not well understood (Timmis et al., 2004; Lane, 2011; Yoshida et al., 2013).

The number of NUPTs we identified in P. trichocarpa is larger than that from a previous study (Yoshida et al., 2013), which used NCBI-BLASTN for NUPT identification instead of LASTZ, as in our study. However, the total length of NUPTs is similar. There is no clear explanation for the variation of the number and amount of NUPTs in plant species, but it may correlate with genome complexity, proportion of repetitive elements, and/or other factors (Yoshida et al., 2013).

Consistent with Yoshida et al. (2013), our analysis shows that large NUPTs are young. Therefore these large NUPTs may experience rapid recombination, insertion/deletion, and fragmentation; supporting previous results that the plant nuclear genome is in equilibrium between frequent integration and rapid elimination of plastid DNA (Matsuo et al., 2005).

There are two main hypotheses to explain the mechanism of DNA transferring from organelle-to-nucleus: the “bulk DNA” (DNA-mediated) view and “cDNA intermediate” (RNA-mediated) view (Timmis et al., 2004). Given the presence of large NUPTs containing both non-coding and coding regions, our results and most previous studies support “DNA-mediated” transfer hypothesis, (Matsuo et al., 2005; Michalovova et al., 2013; Yoshida et al., 2013). Despite the frequent occurrence of large NUPTs, e.g., nearly the entire chloroplast fragment in rice chromosome 10 (Yuan et al., 2002), whether the whole plastid genome can be transferred to the nucleus genome remains largely unknown or not specifically verified yet. We observe that the entire P. trichocarpa chloroplast genome is aligned to a region on chromosome 13 of the P. trichocarpa nuclear genome (Data Sheet S1). Together with MC sampling, our study clearly provide the evidence that the whole chloroplast genome horizontally transfers to the nuclear genome.

As discussed above, the whole chloroplast genome can be transferred to the nucleus, so the missing and pseudogenized genes in the Salicaceae chloroplast genome, e.g., infA, ndhA, rpl16, rpl32, could be transferred to the nucleus and may function properly. Similar cases have been reported previously, involving rpl32 in Populus (Ueda et al., 2007) and infA in some angiosperm species (Millen et al., 2001). Our analyses and other studies show that organellar DNA horizontal transfer in plants is frequent. These transfers are thought to play an important role in gene and genome evolution in plants, and functional transfer of the chloroplast genes may facilitate the regulation of gene expression (Cullis et al., 2009; Yoshida et al., 2013). However, functional gene horizontal transferring from organelle-to-nucleus is rare as the transferred coding sequences must acquire gene promoter and terminator sequences for proper transcription in the nucleus, and must also acquire transit peptides necessary for reimporting the protein back into organelle (Timmis et al., 2004; Stegemann and Bock, 2006). Our result provides additional evidence of the organelle-to-nucleus horizontal functional gene transfer in species of P. trichocarpa, where the majority of cp coding sequences are transferred to the nucleus genome, and most of them remain as coding regions in the nucleus genome (Supplementary Tables S5, S6). However, further investigation is needed to determine whether these plastid-coding sequences in nuclear genome could function properly. Because chloroplast-to-nucleus transferred functional genes must acquire promoter and terminator sequences (Eckardt, 2006), it is necessary to identify the regulatory motifs of these chloroplast-to-nucleus transferred genes. Our discovery that large cp genome segments are transferred into the nuclear genome, implies that the transfer of the coding genes might be accompanied by the transfer of their promoter and terminator sequences, which could provide materials for those cp genes to properly function in the nuclear genome. However, further experiments (e.g., western bolt and genetic screen) and analyses are needed to confirm the function of the genes transferred from the chloroplast to the nucleus. Furthermore, the identification of the whole cp genome transfer event in P. trichocarpa is based on single genome analysis. It remains unclear whether this event is usual or even fixed in P. trichocarpa. Therefore, sequencing the genomes of more individuals is needed to further verify the whole cp genome transfer event in P. trichocarpa.