Ancestral Reconstruction and Investigations of Genomic Recombination on some Pentapetalae Chloroplasts

3.1 Method I: Naked eye investigation

As stated above, this method is not an algorithm that automatically builds the ancestors of the provided genomes, but it is a method applied manually, as follows. Some python codes were produced in order to graphically represent each triplet constituted by two sister species and their closest cousin as three parallel lines, as described in Figure 2 (the three lines at the top of the figure). On each line are located numerous equidistant vertices, one per coding sequence in the associated genome, and all the sequences having the same gene name according to Dogma are linked by a red edge.

Figure 2:

Simulation of ancestral reconstruction process between two genomes.

ADD means that the considered sequence has been added to the list of genes in the ancestral genome.

The ancestral genome of each couple of sister species has then been manually deduced by a consensus approach. Each part shared in common in its genomes is put in the ancestor (see, for instance, the cases of YCF68 and YCF68_1: they are both in Genome-1 and Genome-2, so they were put on the ancestor). In case of a difference, both brother genomes are compared locally with their closest cousin. If the latter agrees with one of the two brothers, the agreement (sequence of genes) is put on the ancestor and this part of the lines is considered as resolved. For instance, PSBI is in Genome-1 but not in Genome-2. As it is in Genome-3, the parcimonious hypothesis is that it was in the ancestor of the 3 species, and that Genome-2 lost it.

If the closest cousin cannot help us to resolve the situation, because the cousin presents locally a third pattern different from the two brothers, then one or more new close cousins are considered recursively. It is the solution that reduces the largest number of rearrangement operations that is finally chosen, leading to the local ancestor gene list (again a parsimonious approach).

Figure 3:

Names of the ancestors in the Apiales and Asterales phylogenetic tree, provided with supports and branch lengths: the numbers shown at each branch are bootstrap scores computed by RAxML [4], while branch lengths are indicated (cf. Br Lengths).

A letter has also been associated to each internal node (the ancestors are in red), while the number of genes per genome is indicated under the label ”Genome length”.

By doing so and verifying our results three times, a trustworthy ancestral list of genes at each internal node was obtained. Our next objective was to recover these ancestors automatically.

3.2 Method II: Ancestor Prediction based on Gene Contents

This method encompasses the following general steps:

1. Step 1: Preliminary stage. In this step, all internal nodes from leaf nodes to the root one are named following the alphabetical order. Each letter in an internal node represents an ancestor genome. This latter can be the ancestor of two leaves, of an internal ancestor and a leaf node, or of two internal ancestors. The result of this step can be seen in Figure 3. In this tree, a bottom-up procedure was applied to predict the ancestor genome at each internal node.

2. Step 2: Genome Selection.Figure 3 presents the most supported topology for two Campanulids subgroups. The algorithm starts by automatically selecting the two closest sister species according to the Needleman-Wunsch distance applied to lists of genes. The other species are then ordered according to their distance in the tree (number of nodes between it and one of the two sister species, and Needleman-Wunsch distance to solve ex-aequo cases), defining what can be called an ordered list of cousins.

3. Step 3: Genes Investigation. In the easiest situation, all gene couples completely match between the two sister species (which thus have the same length). In this case, the frequency of occurrences depends on the considered family, the ancestor is easily deduced as being the same as its children. If there is at least one problematic situation between the selected genomes (that is, if there is at least one deleted, duplicated, or inserted gene in one genome), as for example between E. senticosus and B. hainla, then a deeper investigation is initiated using one or more cousin genome(s). For instance, in this example, M. selavayi and K. septemlobus will be considered first as cousin genomes to take the final decision in the treatment of such problematic situations.

   In this case, all genes are iterated in the considered two brother genomes U₁ and U₂, if gene g_i in U₁ matches properly in name, position, and orientation with gi′ in U₂, then it is added in the ancestor genome γ at position i. Otherwise, consider the gene gi′′ at the same location in the first cousin genome: if g_i or gi′ is equal to gi′′ then add the most frequent gene to the ancestor genome γ in position i, else this gene is considered as an insertion.

   Figure 4 offers a simulation example of the considered procedure. In this figure, suppose that A, B, C, D, and E in the leaves are genes, and the objective is to predict the ancestor α₁. Note that genes A, B, and D match in positions. Concerning the problematic C gene between these two genomes, a cousin will be needed to determine whether it is present in the α₁ ancestor genome or not. One or both genomes in α₂ subtree are considered to be cousin(s) to treat the problem of gene C. The two cousin genomes have one copy of gene C in their gene lists. According to our voting system, gene C will be in α₁ ancestor and a delete operation is recorded (that is, AB_D). An insert state is also marked in α₂ subtree, where gene E did not appear in either cousin genomes of α₁ tree, nor in its brother. Such a deletion is illustrated in Figure 5a and b.

   The simple conflict resolution presented above has been refined by considering the Gestalt pattern matching method [17] based on dynamic programming like Needleman-Wunsch, as it is implemented in the SequenceMatcher method of the difflib Python library.

4. Step 4: Ancestor update. After applying the previous step to all the genes of the sister species, their ancestor is then reconstructed. The subtree of the two brother genomes is then replaced by the list of genes of their ancestor.

5. Step 5: Loop-back. Repeat from Step 2 until the final root ancestor is constructed.

Figure 4:

Simulation of gene investigation step between two genomes

So all matching genes are directly assigned to the ancestor. For non-matching genes, the process consists of the selection of a third genome, among the cousins according to the provided tree. The selected cousin is the closest one to the two considered genomes, according to the chosen distance. It is then compared to the two sister species for each non-matching gene, if the cousin agrees with one sister, then the considered gene is added to the ancestor.

Figure 5:

Graphical presentation of genes alignment between two genomes, useful in the naked eye investigation.

Red curves indicate duplicated genes within a genome, while green lines indicate homology (same gene) between two distinct genomes. (A) At this particular location of the genome, a gene loss can be seen in B. hainla (or one gene gain in K. septemlobus), namely the YCF69. (B) As the same situation occurs between B. hainla and M. delavayi, we can conclude to a gene loss (as gaining twice the same gene at the same location is less likely). (C) At this location, the gene order is not modified between the genomes of the two considered sister species, namely K. septemlobus and M. delavayi.

4.1 The Apiales order

Let us first consider the Apiales order for which both the phylogenetic relationship and the gene lists of leaves are known at the beginning of the study. We then apply the manual and the automatic approaches to infer ancestral states at each internal node of the tree. The results were convergent and lead to the following conclusions regarding the evolution of gene content among the tree. We focused first on the evolution of duplications. Table 2 (supplementary material) contains all duplicated genes among the Apiales order. For each duplication, the number of copies is specified as well. Let us now enter into details regarding the leaves of the phylogenetic tree shown in Figure 3.

1. Sister species E. senticosus and B. hainla have been considered first, with K. septemlobus playing the role of the cousin. After manual and automatic comparisons, we found that the gene YCF1 is present twice in E. senticosus, while it is in three copies in B. hainla. As the cousin has only two sequences of YCF1, we suggest that the latter is present twice in the ancestor, one gene has been inserted in B. hainla. Similarly, YCF68 is in 4 copies in B. hainla and in 6 copies in the sister species. As the cousin presents 6 copies as well, it can be deduced that the common ancestor of E. senticosus and B. hainla contains 6 copies of this gene. In other words, two copies of YCF68 have been removed in B. hainla. All the other genes are similar in both names and locations, and thus the ancestral genome (I) can be deduced.

2. The brother genomes A. undulata and P. ginseng have exactly the same ordered list of genes, which is thus assigned to their last common ancestor (E).

3. Similarly, all couples of sister species A. undulata and P. ginseng, M. delavayi and K.septemlobus, S.delavayi and M. delavayi, and finally K. septemlobus and E. sentucosus match perfectly when considering each couple of brother genomes. In other words, they have not deviated from their respective last common ancestors, which presents the same sequence as their children species. Selected genomes are aligned graphically as shown in Figure 5. We then identify, by using naked eyes investigation and human thinking, the most parsimonious scenario applied on a deduced ancestor, which can lead to these two children using the lowest number of edit operations (such as inserted and deleted genes). Figure 5c shows this matching process applied on M. Delavayi and K. septemlobus, which have a core genome of 169 genes (only the 23 first genes are depicted). Note that, in this example, M. Delavayi (L) and K. septemlobus (M) match completely, so the ancestor H is very easy to obtain (H = L ∩ M has 169 genes).

4. Let us finally compare A. cerefolium and D. carota. The YCF68 gene exists in 2 copies in A. cerefolium while it is missing in D. carota, as shown in Figure 6a. The cousin, for its part, also contains the gene YCF68 (in 6 copies), and so our algorithm concludes the presence of this gene (in 2 copies) in the ancestor of A. cerefolium and D. carota. Additionally, D. carota contains 4 copies of ORF56, while this gene is only represented twice in its sister. As the cousin genome has 4 representatives of ORF56, we can reasonably deduce that this is the case too in the ancestor of these two sister species: two copies of the gene ORF56 have been deleted from the genome A. cerefolium. Such decisions are depicted in Figure 6b, which shows a specific region of the ancestor genome (C). This region has been generated by our algorithm, which has been applied on A. cerefolium and D. carota, and it has been manually validated.

Figure 6:

Green lines indicate homologuous genes while red vertices are for paraloguous ones. C line in figure (B) is the reconstructed ancestor of these two chloroplasts.

(A) Example of gene correspondences in brother genomes D. carota and A. cerefolium. For instance YCF68 is found in position 111 in A. cerefolium while it is missing from D. carota. Additionally, ORF56 is in 2 copies, positions 114 and 115, in D. carota, while this gene is only represented once at position 115 in its sister. (B) Comparison between two brothers. The result is the ancestor genome (C). We can reasonably deduce two copies of ORF56 have been deleted from genome A. cerefolium. The ancestor, for this part, contains also the gene YCF68, which has been deleted from the genome D. carota.

The process detailed above continues with the obtained ancestors and is repeated until reaching the root of the tree: the Last Universal Common Ancestor (LUCA) of Apiales. By operating this reconstruction stage, we found that chloroplasts of this order have not faced so much deletion or insertion in their genomes. Indeed, in most cases, the disparity comes from the variation in numbers of gene copies. The obtained results are summarized in Figure 7.

Figure 7:

Insertion and deletion events found during ancestor reconstruction on Apiales order.

Letters in red refer to ancestor genomes (their lengths are provided too).

4.2 The Asterales order

The Asterales order, which is close to the Apiales one, was then examined. Table 3 in supplementary material contains what has been deduced from our experiments on this order. As can be seen, Asterales genomes have undergone many more changes compared to the Apiales ones. This difference between both orders lead to a larger variation in the lengths of Asterales genomes. For the sake of illustration, let us consider for instance the chloroplast of H. annuus. It only contains 161 coding sequences while its sister species, namely P. argentatum, has 183 genes. The matching process previously described has led in this case to an ancestor of size 162. More precisely, 23 genes have been inserted and two other ones have been deleted in P. argentatum, while only one gene has been removed in H. annuus, as described in Figure 8.

In this second order and in most cases, genes are comparable in both names and locations, having the same positions if we do not consider duplications. Indeed, almost all differences in this set of chloroplastic genomes come from a variation in the number of copies. Let us now investigate a third order, to compare it with Apiales (only a few variations of genomes) and Asterales (large variety in duplications).

4.3 The Fabids order

Figure 8:

Summary of the complete ancestral genomes reconstruction of the Asterales order.

Insertion and deletion events are provided, with names and length of each internal node.

It was easy to deal with Apiales order, while Asterales improved the complexity of the ancestral reconstruction, due to duplications. However, in both cases the proposed algorithm was able to recover results that have been inferred manually (naked eye investigation). Let us now consider a larger and more complicated order, namely the Fabids, to evaluate the performances of our proposal when facing a complex collection of genomes.

Indeed, the main problem with this new order is that it contains large scale inversions in some branches, while it was not the case with both previously studied orders. In this case, a single inversion detection algorithm has been able to signal helpful information regarding such regions, like the beginning and the end (insertion or deletion) of reversals. However, the most difficult case where insertions or deletions are inside the inversion zone is difficult to handle.

Figure 9:

A phylogenetic tree in the reconstruction of the Fabids ancestor and the unambiguous reconstruction accuracies of our algorithms on this tree.

Alphabetic characters represent the ancestors. L stands for the length of each node (number of genes), [D, I] describes the number of deletions and insertion, while inversions are also indicated.

In this situation, our proposition consisted in selecting one of the two brothers to operate as a reference. The best cousin was then searched within the same clade, and the status of each gene in this region (matching, or need insertion or deletion) was then compared. Most of the considered reversal regions match at the genes names level, but with reverse positions. Figure 9 presents our finding on Fabids order, with information about the length of each node. Various rearrangement information were also provided such as the number of insertion, deletion, and inversion.

Figure 10:

The variation in comparison results of ancestral genomes nodes on Asterales order with MLGO.

4.4 Comparison with MLGO

For the sake of comparison, we have examined the ancestral genome contents of both Apiales and Asterales species with MLGO tool, which stands for Maximum Likelihood for Gene Order Analysis^[1]. This latter is, to the best of our knowledge, the first web tool for phylogeny and ancestral genomes reconstruction compatible with genome rearrangements [12].

On the one hand, the ancestors in the Apiales order, provided either by our approach or with MLGO, are very similar in terms of gene contents. However our method outperforms MLGO when investigating the specific location of genes and their number of duplications. On the other hand, the results are very different for some nodes in the Asterales case, as summarized in Figure 10. For instance, when gene YCF1 in internal node d has 2 copies in both its two children and their closest cousin has 2 occurrences of this gene too. In this case, our algorithm proposes to set the number of YCF1 in d to 2, while MLGO produced only one copy.

Figure 11:

The variation in ancestral genomes nodes which were achieved by comparing our method results with MLGO tool on Fabids order.

Similar consequences can be outlined in the most difficult case, namely the Fabids order, as highligthed by Figure 11, each time, our algorithm outperforms the MLGO results by producing what is the most likely ancestral state (numbers of genes and their positions) in each situation. For instance, considering the ancestral node M, it was found that gene INFA is missing in the first child while it is present in the second one. It was also found in both its closest cousins, with one copy at each time. The most reasonable scenario is to consider that the ancestral node under consideration also has a single copy of INFA. This result is produced by our algorithm, while MLGO considers that M must not have INFA in its genome. Other divergent results can be reported, as in the ordinary case of node E: gene ATPF is present once in each of the two childrens. So our algorithm considers that it is present once in E, while with MLGO, this node must contain two copies of ATPF. Other nodes are problematic in the MLGO case, for example, {a, b, c, d} as can be seen in Figure 11. Each time, our algorithm produces results in agreement with the one that have been deduced manually, while in some cases MLGO has yielded surprising results. Indeed, these results were produced following a random extraction among the differences in the order of genes in the current species: on a large random extraction of our large set of reconstructed ancestral genes, leading to 18 different evolutionary scenarios according to our method and that of MLGO, we systematically discovered the right ancestral situation, while MLGO systematically misled itself.

Figure 12:

Green lines indicate homologuous genes while red vertices are for paraloguous ones. C line in figure (B) is the reconstructed ancestor of these two chloroplasts.

(A) Example of comparison with MLGO on Apiales order. We show similarity in gene contents between our results, ancestor node (C), and ancestor node (A1) from MLGO. (B) Apiales order tree produced by MLGO.