Reconstructing the history of a WD40 beta-propeller tandem repeat using a phylogenetically informed algorithm

Preliminaries

A tandem repeat is a sequence composed of the repetition of approximate copies of a motif, where the first or the last repetition may be truncated. For example, the tandem repeat:

CC GTTAC GTTACC GTAGC​

contains approximate repetitions, separated by spaces, of the motif GTTACC.

When the motif and the copies have the same length, we call this length a period. With insertions and deletions in the copies, it is always possible to transform a tandem repeat such that all the complete copies have the same length. For example, the above sequence can be rewritten, by inserting two gaps, as: CC GTTA-C GTTACC GT-AGC. Once all copies have the same length, the trivial multiple alignment of the successive copies is called a self-alignment.

The problem of self-aligning a tandem repeat sequence while optimizing any interesting score—sum-of-pairwise-distances, star and tree alignments, etc.—is NP-hard (Elias, 2006), and is at the heart of every reconstruction strategy (Rivals, 2004). However, in the case of tandem repeats that code for protein domains whose structure is resolved, such as in Fig. 1, structure-guided alignments can be used to infer the corresponding self-alignments (Chaudhuri, Soding & Lupas, 2008).

The history of a tandem repeat is a description of how it went from one copy of the motif in the ancestor to the many copies in the observed sequence. This can be formalized using various models of evolution. All these models include mutations, insertions and deletions of nucleotides that affect a single copy of the motif, and duplications that replace a segment of k copies of the motif by two consecutive segments both identical to the initial segment (Benson & Dong, 1999). More elaborate models will consider operations that transform a segment into two or more identical segments (Gascuel, Bertrand & Elemento, 2005), deletion of copies (Sammeth & Stoye, 2006), or even inversions and speciations (Tremblay Savard, Bertrand & El-Mabrouk, 2011). The main hurdle in the reconstruction problem is the fact that more than a billion years of evolution can significantly blur the similarity of the sequences, thus rendering quantitative comparison tools meaningless.

Algorithmic approaches vary according to the evolutionary model, and more fundamentally according to assumptions on the behavior of duplications, linked to the concept of parsing. Given a period p, a parsing of a tandem repeat is the determination of the boundaries of the copies of a repeated motif of length p (Matroud et al., 2012). For example, here are two different parsings of the example sequence:

CC | GTAA-C | GTTACC | GT-AGC
CCGTA | A-CGTT | ACCGT- | ACC​

In the fixed boundaries model, the boundaries of all duplications are boundaries of the same parsing; in the dynamic boundaries model, parsing may vary from one duplication to the next. With the fixed boundaries hypothesis the reconstruction problem can be reformulated as the construction of a duplication tree, opening up the vast repertoire of phylogenetic tools. In contrast, the more general dynamic boundaries hypothesis implies that exhaustive searches must be repeated at each step (Benson & Dong, 1999), and the evolutionary history cannot be represented by a tree (Jaitly et al., 2002).

Architecture of the GNB1 beta-propeller

The guanine nucleotide-binding protein (GNB1) is a protein characterized by a tandem repeat of seven copies of the WD40 protein motif (PFAM00400). These fold in a three-dimensional cylindrical structure, whose top view, as seen in Fig. 1, looks like a propeller with seven blades. Each blade is a beta sheet consisting of four strands, denoted A, B, C, and D, from the centre of the propeller to its circumference. However, the strands appear in the tandem repeat in the following order:

D1A2B2C2D2A3B3C3D3A4B4C4D4A5B5C5D5A6B6C6D6A7B7C7D7A1B1C1​

Thus, there are at least two competing parsings of the tandem repeat: one that preserves the integrity of the 3D blades, with motif ABCD, and one that maximizes the number of complete copies, with motif DABC. In view of this particularity, we reject the fixed boundaries hypotheses in favor of dynamic boundaries, leaving open the possibility of eventually detecting fixed boundaries.

The WD40 protein motif, which actually codes for the components of a blade, is quite common (Stirnimann et al., 2010), and is responsible for a variety of propellers ranging from 4 to 8 blades, whose ancestry can be traced back to both recent rounds of duplications, and very ancient ones (Chaudhuri, Soding & Lupas, 2008). The GNB1 propeller falls in the latter category, and has been found in all sequenced eukaryotes species.

A GNB1 ortholog tandem repeats dataset

We selected a set of GNB1 orthologs using the Ensembl (Flicek, Amode & Barrell, 2014) database, and BLAST (Altschul et al., 1990) in the case of Arabidopsis. Table 1 gives the accession numbers and the positions of the WD40 tandem repeat in each sequence. Figure 2 shows the current consensus of the phylogeny of the six species (Dunn et al., 2008).

Table 1:

Dataset.

Sequences used to construct the dataset.

Code	Species	Accession	Positions
HUM	Homo sapiens	NP_002065	48–340
DAN	Danio rerio	NP_997774	48–340
DRO	Drosophila melanogaster	NP_996462	48–340
CEL	Caenorhabditis elegans	AAK55963	48–340
NEU	Neurospora crassa	XP_956704	63–357
ARA	Arabidopsis thaliana	NP_001078491	53–364

DOI: 10.7717/peerj-cs.6/table-1

We used the structure-guided alignment of GNB1 given in Chaudhuri, Soding & Lupas (2008) to identify insertions and deletions of the human sequence within the WD40 motif. We then aligned the human sequence with A. thaliana and N. crassa sequences in order to transfer the insertions and deletions from the human, and to identify further deletions in A. thaliana and N. crassa. The alignment is shown in Fig. 3. The three remaining species inherited the same insertions and deletions as the human.

The resulting amino acid sequences have all the same length. These sequences are ortholog tandem repeats with 6.8 copies of the WD40 motif, as shown in the self alignment of the human sequence in Fig. 4. The underlying genomic sequences have a length of 864 bp; in the rest of this paper, they will be denoted by the names HUM, DAN, DRO, CEL, NEU, ARA, as in Table 1.

For comparison purposes, we computed the initial alignment with the standard tools Clustal Omega (Sievers et al., 2011) and Muscle (Edgar, 2004). Of the eleven or so mutations proposed by Chaudhuri et al., Clustal finds 7 of them, and Muscle only 2, we thus used the Clustal alignment. No classical aligner can predict, as Chaudhuri et al. do with structural knowledge, the two deletions that occur in all sequences. These deletions are necessary information for the computation of self-alignments.

Pinpointing the most recent duplication

Here we describe how to detect the most recent duplication in a tandem repeat, first by using species separately, then by processing the six species in parallel. The first model follows closely the Benson & Dong model (1999). The second is an enhanced model that uses phylogenetic trees of ortholog domains to distinguish mutations that occurred before all speciation events, from mutations that occurred after speciation.

Using a single species

After a duplication of a segment, the two copies are identical. If each copy accumulates mutations independently and uniformly, which is equivalent to a molecular clock hypothesis among the repeated segments, the most recent duplication should be indicated by the pair of segments that have the least differences. This leads to the following definition of a cost function, whose value is proportional to the number of mutations that separates two segments, and called the contraction cost (Benson & Dong, 1999). Let S be a nucleotide sequence, ℓ be a period length, and the edit distance given by δ(n₁, n₂) = 0 if nucleotides n₁ and n₂ are equal, otherwise δ(n₁, n₂) = 1. We define the contraction cost of two consecutive segments of length ℓ at position pos of the sequence S as: (1)${\displaystyle C\left(pos,\ell \right)=\frac{\sum _{i=pos}^{pos+\ell -1}\delta \left(S_i,S_{i+\ell }\right)}{\ell }.}$ The normalization of the cost will allow the comparison of costs between different values of ℓ.

Given a motif of length p and a tandem repeat sequence S, the position of the most recent duplication of one or more copies of the motif should be a position in S that minimizes C(pos, ℓ), for all values of pos and ℓ ∈ {k∗p | k ≥ 1} for which C(pos, ℓ) is defined.

For the three GNB1 ortholog sequences HUM, NEU and ARA, we tested all possible periods and obtained minimal values for p = 126. Figure 5 shows the values of C(pos, 126) for each species. The prediction that no meaningful minimal values would emerge is confirmed both by the random nature of the HUM curve, and by the fact that, although these three sequences share the same duplication history, they do not agree on its most recent development.

Using multiple species

In this section, we use the same global approach as the preceding section, but we estimate the cost of a contraction based on the reconstruction of the ancestral tandem repeat of n extant tandem repeats who share the same duplication history, and whose phylogeny is known. We replace the distance δ(n₁, n₂) between nucleotides—whose value is 0 or 1—by a cost that ranges in the interval [0..1]. The following small example illustrates the method.

Consider the following three ortholog sequences, each with a tandem repeat of two repetitions of a motif of length 3.

     1  2  3     1′ 2′ 3′
Sp1  T  A  G  |  T  T  T
Sp2  C  A  A  |  A  T  T
Sp3  T  C  C  |  A  T  T

At positions 1 and 1′, both species Sp2 and Sp3 have a mutation compared to Sp1. If the phylogeny is (Sp1, (Sp2, Sp3)), then all parsimonious reconstructions of the ancestral nucleotide give T for both positions 1 and 1′ as in Tree 1 of Fig. 6. In the case of positions 2 and 2′, all parsimonious reconstructions, such as Trees 2a and 2b of Fig. 6, give A for position 2, and T for position 2′.

For positions 3 and 3′, some reconstructions have a mutation that precedes speciation and some do not. Figure 6 gives four possible reconstructions: Trees 3a, 3b and 3c have one mutation before speciation and Tree 3d does not.

Formally, given n ortholog tandem repeats, from n species whose phylogeny is given by tree $\mathcal{T}$, construct the tree $\mathcal{D}$ whose root is a duplication node, with two copies of tree $\mathcal{T}$ labeling its two children. The leaves of the left subtree are labeled by the nucleotides at position pos for each species, and the leaves of the right subtree are labeled by the nucleotides at position pos + ℓ.

Consider the set $\mathcal{L}\left(pos\right)$ of all labelings of all nodes of tree $\mathcal{D}$ such that the total number of mutations on the edges is minimal. Let m(pos) be the number of labelings in $\mathcal{L}\left(pos\right)$ with one mutation between the root and one of its children. Then the contraction cost for position pos is defined by: $m\left(pos\right)/\left|\mathcal{L}\left(pos\right)\right|$. Note that, in a parsimonious labeling, there is at most one mutation between the root and its children.

The contraction cost of two consecutive segments of length ℓ at position pos of the sequence S is defined by: (2)${\displaystyle C\left(pos,\ell \right)=\frac{\sum _{i=pos}^{pos+\ell -1}m\left(i\right)/\left|\mathcal{L}\left(i\right)\right|}{\ell }}$

We tested this new cost function with the three ortholog tandem repeats ARA, HUM and NEU, in parallel, to obtain the bottom curve of Fig. 5: the positions where the minimal, or near minimal costs are attained are in a small interval, and the contraction costs are nearly halved with respect to the former cost function.

Since the contraction cost for a position, $m\left(pos\right)/\left|\mathcal{L}\left(pos\right)\right|$, is a number between 0 and 1, the total number of mutations as computed by C(pos, ℓ)∗ℓ is not necessarily an integer. When we want to evaluate the number of mutations necessary to perform a contraction, we use the formula: ${\displaystyle M\left(pos,\ell \right)=\lceil C\left(pos,\ell \right)\ast \ell \rceil .}$

Contractions and iterations

The labels at the root of different labelings in $\mathcal{L}\left(pos\right)$ are the putative ancestors of the nucleotides at position pos and pos + ℓ, this set of ancestors is denoted by A(pos).

To perform the contraction of two segments of length ℓ on a tandem repeat at position pos, we replace the two segments by a single segment. This segment is the sequence of sets A(i) for pos ≤ i < pos + ℓ. For example, the contraction of the tandem repeat in the three sequences TAGTTT, CAAATT, and TCCATT used in Fig. 6 would be:

{T}{A,T}{A,C,G,T}

After a contraction is performed, it is possible to iterate the procedure by choosing the next more recent duplication. In this case, the leaves of the tree $\mathcal{D}$ will be labeled by sets of nucleotides instead of single nucleotides.

Clearly, when using multiple species, iterated contractions may be performed only on ancient duplications, that is those which preceded speciation. In the case of GNB1 orthologs, strong evidence supports the fact that all duplications of the WD40 motif preceded all speciations, at least up to the common ancestor of plants and animals. A compelling example of this was shown by Schaper and colleagues in the recent large scale study (Schaper, Gascuel & Anisimova, 2014). All approaches are based on the same principle that, when analyzed by phylogenetic tools, the motifs regroup by their position in the tandem sequence rather than by the species they belong to. We provide, in Fig. 7, yet another evidence of this.

The contraction procedure may be iterated as often as the sequence is long enough to be contracted. However, it is clear that in some circumstances further contractions are not recommended, as the experiment with the HUM sequence showed in Fig. 5. Before performing a contraction, one should make sure that the interval of positions with minimal costs is small, and that these costs are well below the random costs of performing a contraction anywhere else. The next section contains some further examples of these situations.

Algorithms and complexity

The complete procedure is composed of the following steps, taking as input a motif of known length p and a set of n ortholog tandem repeats of the motif:

1. Create a dataset of n tandem repeats of the same length, using insertions and deletions within copies of the motif. If a structure-guided alignment exists, use it.

2. Find the values of (pos, k∗p), k ≥ 1 that minimize, or nearly minimize, the contraction cost C(pos, ℓ), as defined by Eq. (1) or Eq. (2).

3. Decide whether there is sufficient evidence to confirm the position of a further contraction. If not, stop.

4. Perform the contraction, and go to Step 2.

Only Steps 2 and 4 can be solved by algorithms, and, fortunately, efficient algorithms are available for these tasks. The computation of contraction costs as defined by Eq. (1) is rather trivial. The computation of C(pos, ℓ) as defined by Eq. (2), which calls for the enumeration of all parsimonious labeling of a tree is more complex. Fortunately, there exists an algorithm developed by Sankoff (1975) whose data structure can be used to compute both $m\left(pos\right)/\left|\mathcal{L}\left(pos\right)\right|$ and the set of A(pos) of possible ancestors.

Define the cost of a labeling of a tree to be the number of edges that have different labels at their extremities. At each vertex v of the tree $\mathcal{D}$, we compute, for each symbol b in {A, C, T, G, - }, the minimal cost M(v, b) of a labeling of the subtree rooted at v, with base b labeling vertex v; and the number N(v, b) of such minimal cost labelings. A classical dynamic programming algorithm (Sankoff, 1975) is used to fill out the structure. The values of these parameters at the root R of $\mathcal{D}$, and at its left and right children L_R and R_R, are sufficient to compute A(pos) and the contraction cost at position pos. namely:

1. A(pos) is the set of bases b for which M(R, b) is minimal.

2. $\left|\mathcal{L}\left(pos\right)\right|=\sum _{b\in A\left(pos\right)}N\left(R,b\right)$.

3. $m\left(pos\right)=\left|\mathcal{L}\left(pos\right)\right|-\sum _{b\in E}N\left(L_R,b\right)\cdot N\left(R_R,b\right),\text{ where }E=\left\{b|b\in A\left(pos\right)\text{ and }M\left(R,b\right)=M\left(L_R,b\right)+M\left(R_R,b\right)\right\}$.

These are all polynomial algorithms, but the extensive computations involved in making tens of thousand of experiments prompted us to implement these with various speed-up techniques (Lavoie-Mongrain, 2014).

The 4-bladed ancestor

Of the 27,000 experiments, 6 optimal scenarios predicted a total of 83 mutations to perform three contractions, and 66 near optimal scenarios predicted a total of 84 mutations. Table 3 shows the same statistics as Table 2 for these 72 scenarios. The complete results for these experiments are in Table S1.

We analyzed all 72 reconstructions of the 4-bladed ancestor corresponding to the 72 best results. These reconstructions are remarkably similar, which is not completely surprising given the narrow ranges of possible first, second and third contraction, as can be seen in Table 3, and Fig. S1.

For each triplet of positions, we performed the three contractions, and then translated the contracted part using the extended one-letter amino acid code (Cornish-Bowden, 1985) that assigns the letter B to the sets of amino acids {D, N}, and the letter Z to the set {E, Q}, see Fig. 9.

All reconstructions left undisturbed both the first two repeats—positions 1–251—and the last 39 positions (see Fig. S1). All reconstructions agreed on the following translated middle part sequence (see Fig. S2), where dots indicate unresolved amino acids:

GHTG.V....F.PBD.....G..D.TC..WD...GZ.V.....H.N......FS..G​

and 69 sequences (96%) agreed with the pattern:

GHTG.V....F.PBD..L..G..D.TC..WD...GZ.V.....H.N......FS..G​

Searching these patterns in various protein databases yields no result. However, searching for the longest (sub)pattern that has at least one hit, gives very interesting results. This pattern is:

L..G..D.TC..WD...GZ.V.....H.N​

The pattern was not found in the six sequences that were chosen for the current experiment. However, there were 32 hits in the TrEMBL protein database (Bairoch & Apweiler, 2000) using ScanProsite (de Castro et al., 2006). Of these hits, two were for the slime molds Dictyostelium discoideum and Dictyostelium purpureum, whose divergence is presumed to date back to 400 million years (Parikh et al., 2010). The rest of the hits are all Basidiomycota, which is a fungus phylum: these hits include mostly rots and smuts, with some yeasts, and two edible mushrooms, the Shitake and the Champignon de Paris.