Adel Ait-Hamlat

Diego Javier Zea

Antoine Labeeuw

Lélia Polit

Hugues Richard
email: hugues.richard@upmc.fr

Elodie Laine
email: elodie.laine@upmc.fr

Transcripts' evolutionary history and structural dynamics give mechanistic insights into the functional diversity of the JNK family

Keywords: Alternative splicing, Molecular modeling, Evolution, Transcript phylogeny, Kinase

Alternative splicing and alternative initiation/termination transcription sites, have the potential to greatly expand the proteome in eukaryotes by producing several transcript isoforms from the same gene. Although these mechanisms are well described at the genomic level, little is known about their contribution to protein evolution and their impact at the protein structure level. Here, we address both issues by reconstructing the evolutionary history of transcripts and by modeling the tertiary structures of the corresponding protein isoforms. We reconstruct phylogenetic forests relating 60 protein coding transcripts from the c-Jun N-terminal kinase (JNK) family observed in 7 species. We identify two alternative splicing events of ancient origin and show that they induce subtle changes on the protein's structural dynamics. We highlight a previously uncharacterized transcript whose predicted structure seems stable in solution. We further demonstrate that orphan transcripts, for which no phylogeny could be reconstructed, display peculiar sequence and structural properties. Our approach is implemented in PhyloSofS 1 Evolution and structural impact of alternative splicing for JNKs.

2

(Phylogenies of Splicing Isoforms Structures), a fully automated computational tool freely available at https://github.com/PhyloSofS-Team/PhyloSofS.

Introduction

Alternative splicing (AS) of pre-mRNA transcripts and alternative transcription initiation/termination are essential eukaryotic regulatory processes. They can impact the regulation of gene expression, for instance by introducing changes in the three prime untranslated region [START_REF] Sandberg | Proliferating cells express mrnas with shortened 3' untranslated regions and fewer microrna target sites[END_REF] . Or they can directly modify the content of the coding sequence (CDS) [START_REF] Graveley | Alternative splicing: increasing diversity in the proteomic world[END_REF] , leading to different protein isoforms. Virtually all multi-exons genes in vertebrates are subject to AS [START_REF] Wang | Alternative isoform regulation in human tissue transcriptomes[END_REF] and about 25% of the AS events (ASEs) common to human and mouse are conserved in other vertebrates 3;49;51 . This suggests an important role for AS in expanding the protein repertoire through evolution. AS has also gained interest for medicinal purpose, as the ratio of alternatively spliced isoforms is imbalanced in several cancers 43;70 .

The extent to which the ASEs detected at the gene level actually result in functional protein isoforms in the cell remains largely unknown. Transcriptomics and proteomics studies suggested that most highly expressed human genes have only one single dominant isoform 22;25 , but the detection rate of these experiments is very difficult to assess [START_REF] Kim | A draft map of the human proteome[END_REF] and likely suffer from strong experimental detection bias [START_REF] Wang | Detection of proteome diversity resulted from alternative splicing is limited by trypsin cleavage specificity[END_REF] . A recent analysis of ribosome profiling data suggested that a major fraction of splice variants is translated, with direct implications on specific cellular functions [START_REF] Weatheritt | The ribosomeengaged landscape of alternative splicing[END_REF] . Moreover, a large scale assessment of isoforms present in the cell revealed that the majority of isoform pairs share less than 50% of their interactions [START_REF] Yang | Widespread expansion of protein interaction capabilities by alternative splicing[END_REF] . From a structural perspective, very few alternatively spliced isoforms have been characterized and are available in the Protein Data Bank (PDB) 7;27 . It was shown that the boundaries of single constitutive exons or of co-occurring exon pairs tend to overlap those of compact structural units, called protein units [START_REF] Gelly | Selective constraint on human pre-mRNA splicing by protein structural properties[END_REF] . Moreover, tissuespecific alternatively spliced exons are enriched in disordered regions containing binding motifs 13 . It was also suggested that splicing events may induce major fold changes 9;10 , and a few cases of isoforms displaying domain atrophy while retaining some activity have been reported [START_REF] Prakash | Domain atrophy creates rare cases of functional partial protein domains[END_REF] .

The elusiveness of the significance of AS for protein function and fold diversification though evolution has stimulated the development of knowledge bases, such as APPRIS [START_REF] Rodriguez | APPRIS: annotation of principal and alternative splice isoforms[END_REF] and Exon Ontology [START_REF] Tranchevent | Identification of protein features encoded by alternative exons using Exon Ontology[END_REF] . They provide functional and sequence-based information at the level of the transcript or the exon. A method reconstructing transcripts' phylogenies was also proposed and proved useful for enhancing transcriptome reconstruction from ESTs and investigating proteins functional features (domains, sites) 16;17 .

In this work, we combine sequence-and structure-based information to shed led on the evolution of AS. We have developed PhyloSofS (Phylogenies of Splicing Isoforms Structures), an automated tool that infers plausible evolutionary scenarios explaining an ensemble of protein coding transcripts observed in a set of species and predicts the tertiary structures of the protein isoforms. We show how PhyloSofS can be used to identify and date ASEs, and also shed light on the molecular mechanisms underlying their functional outcome in the c-Jun N-terminal kinase (JNK) family. This choice was motivated by the fact that JNKs are among the few tens of families for which alternative transcripts performing different functional tasks have been experimentally identified and characterized 8;36;64 . Moreover, they play essential regulatory roles by targeting specific transcription factors (c-Jun, ATF2...) in response to cellular stimuli. The deregulation of their activity is associated with various diseases (cancer, inflammatory diseases, neuronal disorder...) which makes them important therapeutic targets [START_REF] Manning | Targeting JNK for therapeutic benefit: from junk to gold?[END_REF] . About ten JNK splicing isoforms have been documented in the literature [START_REF] Kyriakis | Mammalian MAPK signal transduction pathways activated by stress and inflammation: a 10-year update[END_REF] . They were shown to perform different context-specific tasks 12;30;32;66 and to have different affinities for their substrates 11;67 . By reconstructing the phylogeny of JNK transcripts across seven species, we identify two ASEs of ancient origin. We further identify key residues that may be responsible for the selective recognition of JNK substrates by different isoforms and characterize the behaviour of these isoforms in solution by biomolecular simulations. One of the ASEs involves a 80-residue deletion and has never been documented before. We find that its predicted structure is stable in solution. Both ASEs are supported by sequencing evidence from transcriptomics studies.

Our work allows to put together, for the first time, two types of information, one coming from the reconstructed phylogeny of transcripts and the other from the structural modeling of the produced isoforms, and this to shed light on the molecular mechanisms underlying the evolution of protein function. It goes beyond simple conservation analysis, by dating the appearance of ASEs in evolution, and beyond general structural considerations regarding AS, by characterizing in details the isoforms' shapes and motions. We find that the effect of functional ASEs on the structural dynamics of the isoforms may be subtle and require such a detailed investigation. Our results also open the way to the identification and characterization of new isoforms that may be targeted in the future for medicinal purpose.

Computational method

PhyloSofS can be applied to single genes or to gene families. Given a gene tree and the observed protein coding transcripts at the leaves (Fig. 1a, on the left), it reconstructs a phylogenetic forest embedded in the gene tree (Fig. 1a, on the right) representing plausible evolutionary scenarios explaining the transcripts. Each tree in the forest (in orange, green or purple) represents the phylogeny of one transcript. In other words, each root indicates the appearance of a new transcript and its corresponding ASE(s).

Transcript losses are possible (triangles in Fig. 1a), and the exon usage of a transcript can change along the branches upon the inclusion or exclusion of one or several exons (which we refer to as "mutations"). The underlying evolutionary model is comprised of two levels, following [START_REF] Christinat | Inferring transcript phylogenies[END_REF] . At the level of the gene, exons can be absent, constitutive, or alternative (i.e. involved in at least one ASE), whereas at the level of the transcript, exons are either present or absent. The cost associated to the mutation of an exon naturally depends on its impact on the status of the exon at the gene level. For instance when the gain of an exon at the gene level shifts its status from absent to constitutive, the mutation will not be penalized (see Methods and Table II).

PhyloSofS algorithm seeks to determine the scenario with the smallest number of evolutionary events, following the maximum parsimony principle. It is inspired from that reported in [START_REF] Christinat | Inferring transcript phylogenies[END_REF] . Our main contribution was to develop heuristics in order to treat complex cases in a computationally tractable way. Specifically, we have implemented a multistart iterative strategy combined with a systematic local exploration around the best current solution to efficiently search the space of phylogenetic forests (see Methods and Supplementary Fig. S1 ). Moreover, we have designed a branch-and-bound algorithm adapted to the problem of assigning transcripts between parent and child nodes (see Methods and Supplementary Text S1 ). The reconstructed forests are provided with a userfriendly visualization (Fig. 1b-c). In addition to phylogenetic reconstruction, PhyloSofS predicts the 3D structures of the protein isoforms. The predictions are performed based on comparative modeling using the HH-suite [START_REF] Hildebrand | Fast and accurate automatic structure prediction with HHpred[END_REF] . Furthermore, PhyloSofS annotates the generated models with sequence (exon boundaries) and structure (secondary structure, solvent accessibility, model quality) information. For example, it is very easy to visualize the location of each exon on the modeled structure. Here, we present the application of PhyloSofS to the c-Jun N-terminal kinase (JNK) family across 7 species (H. sapiens, M. musculus, X. tropicalis, T. rubripes, D. rerio, D. melanogaster and C. elegans). This case represents a high degree of complexity with 60 observed transcripts assembled from a total of 19 different exons. Most of these transcripts comprise more than 10 exons and the number of transcripts per gene per species varies from 1 to 8 (Fig. 1b-c). The forest is comprised of 7 trees, 19 deaths (triangles) and 14 orphan transcripts (in grey). Mutation events are indicated on branches by the symbol + orfollowed by the number of the exon being included or excluded (e.g. +11 ). The cost of the phylogeny is 69 (with C B = 3, C D = 0 and σ = 2). On the top right corner are displayed the exons present in each current species (in black). On the bottom left corner are displayed the exon compositions of the human isoforms for which a phylogeny could be reconstructed.

+ 1 0 ' -1 3 -12 +1 4 -1 0-13 JNK2 JNK3 +1 4-1 2 -1 3 -10 α β δ Zebra-a Fugu-a JNK1 0 1' 1

Results

Transcripts' phylogeny for the JNK family.

The observed transcripts were collected from the Ensembl 75 database (see Methods). Phy-loSofS algorithm was run for 10 6 iterations on the JNK family gene tree and we retained the most parsimonious evolutionary scenario (cost = 69, see Methods for a detailed description of the parameters). The reconstructed forest is comprised of 7 transcript trees (Fig. 1b, each tree is colored differently). Each transcript is described as a collection of exons, numbered from 0 to 14 (Fig. 1b, top right and bottom left corners, and see Methods for more details on the numbering). We could reconstruct a phylogeny for 46 out of the 60 observed transcripts. The 14 "orphan" transcripts (leaves in grey) are not conserved across the studied species, and thus likely result in non-functional protein products. Mutations occurring along the branches of the trees are labelled (Fig. 1b, see +/-symbols followed by the number of the included/excluded exon). In total, JNK transcripts' phylogeny comprises 11 mutations.

The sequences of the JNK genes are highly conserved through evolution (Table I).

While Drosophila melanogaster and nematode are the most distant species to human, their unique JNK genes share as much as 78% and 56% sequence identity, respectively, with human JNK1 (Table I). The sequence identities with human JNK2 and JNK3 are slightly lower (Table I, in grey). This suggests that the most recent common ancestor of the 7 studied species contained one copy of an ancestral JNK1 gene. Under this assumption, we propose an evolutionary scenario to reconcile the JNK family gene tree (Supplementary Fig. S2a) and the species tree (Supplementary Fig. S2b). In this scenario, early duplication events led to the creation of JNK2 and JNK3 in the ancestor common to mammals, amphibians and fishes (Supplementary Fig. S2b). JNK1 was then further duplicated in fishes while JNK2 was lost in Xenopus tropicalis. One can see that there is no conflict between the gene and species trees under these hypotheses (Supplementary Fig. S2a-b). genes (Fig. 1b). The transcripts of the same color belong to the same tree and share the same exon composition, even if they come from different gene loci and hence have different amino acid sequences. For instance, the transcript structure including exons 6, 8

and 12 and excluding exons 0, 1', 7 and 13 (in yellow) is shared by 3 human transcripts present in JNK1, JNK2 and JNK3 (Supplementary Fig. S2c). Note that this may not be the case in general, for any protein family: the leaves of a tree may have different exon compositions if mutations occur along the branches.

Two pairs of exons, namely 6-7 and 12-13, are mutually exclusive (Supplementary Fig. S2c). The associated ASEs can be dated early in the phylogeny (Fig. 1b), before the gene duplication (Supplementary Fig. S2b). Neither Drosophila melanogaster nor nematode contain any of exons 12-13. Hence, it is equivalent to consider that exon 12

or exon 13 appeared first (compare Fig. 1b and Supplementary Fig. S3 ). By contrast, exon 7 is clearly predicted as appearing before exon 6 (Fig. 1b, compare purple tree with yellow and orange trees). Noticeably, The two transcripts expressing exons 6 (in orange and yellow) are consistently absent from Zebrafish JNK1b and Xenopus tropicalis parsimonious explanation would be that the gene annotation in the Ensembl database is incomplete. We searched for direct experimental evidence of the expression of the transcripts in these two organisms using transcriptome sequencing data from hundreds of RNA-Seq libraries (see Methods). In Xenopus tropicalis, the analysis of exon-exon junctions revealed the expression of transcripts containing a 72bp-long exon with a translated sequence very similar to that of exon 6 in other species. The sequencing support for this exon is strong as it is present in more than two third of the Xenopus tropicalis RNA-seq libraries we studied. Additionally, a transcript containing this exon is predicted in Refseq (Refseq ID: XM 012966153.2) and the corresponding genomic region is strongly conserved. In D. rerio, the analysis of exon-exon junctions in JNK1b also identified one new 72bp-long exon with a translated sequence very similar to exon 6 in other species.

Hence, there is significant evidence of the expression of transcripts containing exon 6 in both X. tropicalis JNK1 and D. rerio JNK1b, although they are not annotated in Ensembl. This observation gives support to our choice of not penalizing death (C D = 0) when we reconstruct the transcripts' phylogeny as a way to account for the incomplete transcript data.

Among the three transcripts appearing after the gene duplication events (Fig. 1b, in pink, green, and red), one transcript features a large deletion encompassing exons 6, 7

and 8 (JNK1 sub-forest, internal node A11, in pink). Its exon composition is perfectly conserved along the phylogeny (no mutation). We looked for additional RNA-seq support for this transcript and found evidence in 3 Human RNA-seq libraries (out of 166) for reads aligning to the exon-exon junction between exon 7 and 8'. There was no evidence in the RNA-Seq mouse libraries. The two other transcripts are created at the root of the JNK3 sub-forest (ancestor node 10, in green and red). They are characterized by the presence of exons 0 and 1', not found in the other paralogues. Another characteristic feature can be observed for the JNK3 gene, namely exon 7 is completely absent from the associated sub-forest. The genomic sequence of exon 7 is present at the JNK3 locus in all species, but it diverged far more in this gene compared to JNK1 and JNK2.

Mapping of the gene 1D structure onto the protein 3D structure.

About eighty structures of human JNKs are available in the PDB (Supplementary Table S1 ). This abundance of structural data can be explained by the fact that JNKs are important therapeutic targets and they were crystallized with different inhibitors. The three paralogues share the same fold, which is highly conserved among protein kinases. The structures are highly redundant, with an average root mean square deviation (RMSD) of 1.96 ± 0.71 Å, computed over more than 80% of the protein residues. In order to visualize the correspondence between the gene structure and the protein secondary and tertiary structures, the exons were mapped onto a high-resolution PDB structure (3ELJ 15 ) of human JNK1 (Fig. 2, each exon is colored differently). One can observe that the organization of the protein 3D structure is preserved by the 1D structure of the gene. Most of the secondary structures (10 over 12 α-helices and 7 over 9 β-strands) are completely included in single exons. Moreover, each one of the regions important for the structural stability and/or function of protein kinases (Fig. 2, labelled in black) is included in one single exon (see also Supplementary Table S2 ). So are the N-terminal hairpin and the MAPK insert (labelled in grey), two structural motifs specific to the mitogen-activated protein kinase (MAPK) type, to which the JNKs belong. By contrast, binding sites for cofactors and substrates (green circles, see also Supplementary Table S2 ) are comprised of residues belonging to different exons. This is expected as binding sites are comprised of segments that can be very far from each other along the protein sequence. Of note, the block formed by exons 1 to 5, comprising the N-terminal lobe and the A(ctivation)-loop (Fig. 2, from blue to white), is constitutively present in all transcripts belonging to the colored trees on Fig. 1b. The correspondence was also analyzed for the JNK protein from Drosophila melanogaster, whose 3D structure is very similar to that of human JNK1 (Supplementary Fig. S4, RMSD of 0.68 Å).The JNK gene from the Drosophila melanogaster genome comprises much fewer exons than the human gene. The match between the borders of these exons and the borders of the secondary structures and known important regions is even better in that species. Considering the high degree of conservation of JNK sequences, one may hypothesize that a good match also exists in all studied species. Our observation is in agreement with a previous study establishing a relationship between exon boundaries and structurally consistent protein regions [START_REF] Gelly | Selective constraint on human pre-mRNA splicing by protein structural properties[END_REF] .

Properties of the orphan transcripts.

We investigated whether the orphan transcripts, for which no phylogeny could be reconstructed (Fig. 1b, grey leaves), displayed peculiar sequence and structural properties compared to the "parented" transcripts (Fig. 1b, colored leaves). Our assumption is that an orphan transcript is less likely to have functional importance. First, the orphan transcripts are significantly smaller than the parented ones (Fig. 3a). While the minimum length for parented transcripts is 308 residues, with an average of 406 ± 40 residues (Fig. 3a, in white), the orphan transcripts can be as small as 124 residues, with an average of 280 ± 88 residues (Fig. 3a, in grey). Second, regarding secondary structure content, both types of transcripts contain about 40% of residues predicted in α-helices or β-sheets (Fig. 3b). Third, the 3D models generated by PhyloSofS molecular modeling routine for the orphan transcript isoforms are of poorer quality than those for the transcripts belonging to a phylogeny (Fig. 3c-d). The quality of the models was assessed by computing Procheck 42 G-factor and Modeller 47 normalized DOPE score (Fig. 3c-d). A model resembling experimental structures deposited in the PDB should have a G-factor greater than -0.5 (the higher the better) and a normalized DOPE score lower than -1 (the lower the better). The distributions obtained for the parented isoforms are clearly shifted toward better values and are more narrow than those for the orphan transcripts. Finally, the proportion of protein residues being exposed to the solvent (relative accessible surface area rsa > 25%) is significantly higher for the orphan isoforms (Fig. 3e), as is the proportion of hydrophobic residues being exposed to the solvent (Fig. 3f). Overall, these observations suggest that simple sequence and structure descriptors enable to distinguish the orphan transcripts from the ones within a phylogeny and that the formers display properties likely reflecting structural instability (large truncations, poorer quality, larger and more hydrophobic surfaces).

Subtle changes in the protein's internal dynamics linked to substrate differential affinity.

The two mutually exclusive exons 6 and 7 are particularly important for JNK cellular functions, as they confer substrate specificity. The inclusion or exclusion of one or the other results in different substrate-binding affinities 11;67 . From a sequence perspective, the two exons are homologous, highly conserved through evolution, and differ only by a few positions (Supplementary Fig. S5 ). From a structural perspective, they both fold into an α-helix, known as the F-helix, followed by a loop (Fig. 2, in light pink). The F-helix was shown to play a central role in the structural stability and catalytic activity of protein kinases 38;53 . It serves as an anchor for two clusters of hydrophobic residues, namely the catalytic and regulatory spines (see illustration on the PKA kinase on Supplementary Fig. complex on Supplementary Fig. S7a). In the following, we will use these known structural features as proxies for the stability and catalytic competence of the studied isoforms.

The available JNK crystallographic structures and the 3D models generated by Phy-loSofS do not display any significant structural change upon exchanging exons 6 and 7. The catalytic and regulatory spines, together with their anchors in the F-helix, are present in both types of isoforms (Supplementary Fig. S6b-c). The N-terminal aspartate (D207) of the F-helix, which serves as an anchor for the spines, is 100% conserved in both exons 6 and 7 in the 7 studied species (Supplementary Fig. S5, indicated by an arrow).

The two other anchor points are also present, namely I214 and L/M218 (Supplementary Fig. S5, indicated by arrows). Moreover, the characteristic H-bond pattern with the HRD motif and the associated strained backbone conformation are also observed in both types of isoforms (Supplementary Fig. S7b-c). Consequently, both exons 6 and 7, and thus the isoforms containing them, possess the structural features known to be important for kinase catalytic activity and/or regulation.

To further investigate the potential impact of the inclusion/exclusion of exon 6 or 7 on the dynamical behavior of the protein, we performed all-atom molecular dynamics (MD) simulations of the human isoforms colored in orange and purple on Figure 1b. We shall refer to these isoforms as JNK1α (with exon 6 ) and JNK1β (with exon 7 ), in agreement with the nomenclature found in the literature [START_REF] Waetzig | Context-specific inhibition of JNKs: overcoming the dilemma of protection and damage[END_REF] . JNK1α and JNK1β were simulated in explicit solvent for 250 ns (5 replicates of 50 ns, see Methods). The backbone atomic fluctuation profiles of the two isoforms are very similar (Fig. 4a, orange and purple curves), except for the A-loop which is significantly more flexible in JNK1α: the region from residue 176 to 188 displays averaged Cα fluctuations of 1.55 ± 0.28 Å in JNK1α and of 0.98 ± 0.16 Å in JNK1β (Fig. 4a). We should stress that this loop displays the highest deviations among the JNK structures available in the PDB and often comprises unresolved residues. The two exons, 6 and 7, have similar backbone flexibility. In the F-helix, the anchor residues for the spines, D207, I214 and M218 adopt stable and very similar conformations (Fig. 4b). Moreover, the HRD backbone strain and the associated H-bond pattern are maintained along the simulations of both systems (Supplementary Fig. S8a-b). Consequently, the observations realized on the static 3D models hold true when simulating their dynamical behavior: the 6 /7 variation does not induce any drastic change on the protein's overall shape and behaviour.

Nevertheless, we observe differences in the side-chain flexibilities of a few residues lying in the loop following the F-helix between the two isoforms (Fig. 4b). On the one hand, in exon 6 (in orange), the polar and positively charged residues H221, K222 and R228 are exposed to the solvent and display large amplitude side-chain motions. These amino acids are 100% conserved in exon 6 across all species (Supplementary Fig. S5 ).

On the other hand, in exon 7 (Fig. 4b, in purple), G221, G222 and T228 have small side chains with much reduced motions. While G221 is conserved across all species, position 222 is variable and position 228 features G, T or S (Supplementary Fig. S5 ). This region of the protein is involved in the binding of substrates (see Fig. 2, F-site). Moreover, in both isoforms, we predicted residues 223-230 as directly interacting with cellular partners (see Methods). Consequently, one may hypothesize that the differences highlighted here may be crucial for substrate molecular recognition specificity. The positive charges, high fluctuations, high solvent accessibility and high conservation of residues H221, K222 and R228 in JNK1α support a determinant role for these residues in selectively recognizing specific substrates.

Structural dynamics of a newly identified isoform.

Our reconstruction of the JNK transcripts' phylogeny highlighted a JNK1 isoform (Figure 1b, in pink) that has not been documented in the literature so far. It is expressed in human, mouse and fugu fish (Figure 1b), suggesting that it could play a functional role in the cell. To investigate this hypothesis, we analyzed the 3D structure and dynamical behavior of this isoform in human. We refer to it as JNK1δ.

JNK1δ displays a large deletion (of about 80 residues), lacking exons 6, 7 and 8. It does not contain the F-helix, shown to be crucial for kinases structural stability [START_REF] Kornev | Surface comparison of active and inactive protein kinases identifies a conserved activation mechanism[END_REF] , nor the MAPK insert, involved in the binding of the phosphatase MKP7 [START_REF] Liu | A conserved motif in JNK/p38-specific MAPK phosphatases as a determinant for JNK1 recognition and inactivation[END_REF] (Fig. 2). The 3D model generated by PhyloSofS superimposes well to those of JNK1α and JNK1β, with a RMSD lower than 0.5 Å on 245 residues. This is somewhat expected as we use homology modeling. Nevertheless, cases were reported in the literature where homology modeling detected big changes in protein structures induced by exon skipping [START_REF] Nicolas | Becker muscular dystrophy severity is linked to the structure of dystrophin[END_REF] . In the model of JNK1δ, the F-helix present in JNK1α and JNK1β (residues 207 to 220) is replaced by a loop (residues 282 to 288) corresponding to exon 8' (Fig. 4c, indicated by the two stars). The sequence of this loop (exon 8') does not share any significant identity with the F-helix (N-terminal parts of exons 6 and 7 ), except for the N-terminal residue which is an aspartate, namely D282 (D207 in JNK1α and JNK1β). This replacement results in the regulatory spine being intact in JNK1δ (Supplementary Fig. S6d, in red). Moreover, the HRD motif's strained backbone conformation and the associated H-bond pattern, which are stabilized by the aspartate, are maintained (Supplementary Fig. S7d). By contrast, the catalytic spine lacks its two anchors (Supplementary Fig. S6d, in yellow).

JNK1δ was simulated in explicit solvent for 250 ns (5 replicates of 50 ns). The isoform displays stable secondary structures (Supplementary Fig. S9, at the bottom) and atomic fluctuations comparable to those of JNK1α and JNK1β (Fig. 4a, pink curve to be compared with the purple and orange curves). The Cα atomic fluctuations averaged over the loop replacing the F-helix are of 0.88 ± 0.18 Å. This is higher than the values computed for the F-helix in JNK1α and JNK1β (0.57 ± 0.10 Å and 0.53 ± 0.09 Å), but it still indicates a limited flexibility. Moreover, the N-terminal aspartate D282 establishes stable H-bonds with the HRD motif along all but one of the replicates (Supplementary Fig. S8a, on the right) and the HRD motif's backbone remains in a strained conformation (Supplementary Fig. S8b, on the right), as was observed for JNK1α and JNK1β. Consequently, JNK1δ seems stable in solution, and, as observed on the static 3D model, the absence of the F-helix in this isoform is partially compensated by the presence of D282, which is sufficient to maintain H-bonds with the HRD motif and a resulting backbone strain of the motif, important for kinase structural stability.

The main difference between JNK1δ and the two other isoforms lies in the amplitude of the motions of the A-loop. In JNK1δ, the C-terminal part of the A-loop can detach from the rest of the protein along the simulations (Fig. 4c). The amplitude of the angle computed between the most retracted conformation (in grey) and the most extended one (in black) is 107 • . By contrast, in JNK1α and JNK1β, the A-loop always stays close to the rest of the protein, with amplitude angles of 18 • and 19 • , respectively. The A-loop contains two residues, T183 and Y185 (Fig. 4c, highlighted in sticks), whose phosphorylation is required for JNK activation. We hypothesize that the large amplitude motion in JNK1δ might favor their accessibility and, in turn, the activation of the protein. Unresolved residues in the 3D models.

In the 3D models generated by PhyloSofS, the N-terminal exons 0 and 1' and the Cterminal exons 12 and 13 are systematically missing. This is due to the lack of structural templates for these regions. Using a threading approach instead of PhyloSofS homology modeling routine (see Methods) did not enable to improve their reconstruction. In fact, the models generated by the threading algorithm are very similar to those generated by PhyloSofS.

At the C-terminus, exons 12 and 13 are completely predicted as intrinsically disordered (Supplementary Fig. S12a and Supplementary Fig. S12b, blue curve). At the N-terminus, exons 0 and 1' contain two segments of about 10 residues predicted as disordered protein-binding regions (Supplementary Fig. S12b, orange curve), i.e regions unable to form enough favorable intra-chain interactions to fold on their own and likely stabilized upon interaction with a globular protein partner [START_REF] Meszaros | Prediction of protein binding regions in disordered proteins[END_REF] . These exons are present in only two JNK3 transcript isoforms (Fig. 1b, colored in dark red and green).

Alternative phylogenies and robustness to parameter changes

PhyloSofS phylogenetic reconstruction's algorithm may find several solutions with equivalent costs, depending on the input data and the set of parameters. The forest described above (Fig. 1b, and Supplementary Fig. S3 with branch swapping), comprising 7 trees, 19 deaths and 14 orphans, was visited 1 219 times over 10 6 iterations of the program. An alternative phylogeny was visited 310 times, that comprises the same number of trees and orphans, but 2 more deaths (Supplementary Fig. S10 ). The difference between the two forests lies among the fugu JNK1 transcripts, where one transcript belongs to the orange tree (Supplementary Fig. S10 ) instead of the yellow one (Fig. 1b). The two trees differ by the inclusion or exclusion of exon 12 or 13, and the re-assigned transcript lacks both exons. Consequently, the new branching results in the loss of exon 13 between the internal nodes A11 and A18 (Supplementary Fig. S10 ), instead of the loss of exon 12 between A24 and fugu JNK1 (Fig. 1b). Another forest with the same cost comprising 8 trees, 23 deaths and 13 orphans was visited 190 times (Supplementary Fig. S11 ). The additional tree is created in the internal node A10 and links two observed JNK3 transcripts: one from the mouse that was previously orphan (Fig. 1b) and one from zebrafish that previously belonged to the green tree. These two transcripts are very similar to the green and dark red transcripts, in terms of exon composition and of structural properties. The only difference is that they lack exons 12 and 13. Consequently, this new branching avoids the loss of exon 12 between A16 and zebrafish JNK3. Overall the differences between the three solutions are minor and these ambiguities do not impact our interpretation of the results.

To further assess the robustness of our results, we ran PhyloSofS algorithm with different parameters and analyzed the output phylogenies. The three main parameters of the algorithm are the costs C B , C D and σ, associated to the creation (birth), the loss (death)

and the mutation of a transcript, respectively. The forests described above were obtained by setting the death cost to zero (C D = 0) and the ratio between the birth and mutation costs below 2 (C B /σ = 3/2 = 1.5). The choice of not penalizing death was motivated by the fact that the transcriptome data and annotations we are working with may be incomplete. Indeed, the different genomes available in Ensembl are not annotated with the same accuracy. The choice of tolerating few mutations within each tree was motivated by the fact that several pairs of transcripts from the same species differ by only two exons (Fig. 1b, compare exon compositions in the bottom left corner). Varying slightly the birth-to-mutation ratio while maintaining the death cost to zero did not impact the results (Fig. 5 and Supplementary Table S3, compare combinations 5.0.3 and 5.0.4 with 3.0.2). A bigger parameter change, resulting in a birth-to-mutation ratio larger than 2, had a minor impact on the results (Fig. 5, combination 5.0.2). The numbers of orphans and deaths were slightly modified, but the main information contained in the phylogenies remained the same (Supplementary Table S3 ). Penalizing death reduced the number of deaths, as expected, and increased the number of orphans (Fig. S3 ), except when the number of trees was constrained by a high birth cost (combination 5.3.2, see also Fig. 5). These observations illustrate well how the presence of under-annotated species in the input data may artificially scatter very similar or identical transcripts in different trees. We also investigated the influence of gene-level changes on the phylogenetic reconstruction. To do this, we changed the priority rule for the determination of the exon states at the gene level and we changed the cost associated to gene-level induced transcript mutations from zero to σ. The former did not have any impact on the results (Supplementary Table S3, combination 3.0.2 a ). The latter resulted in more trees and more deaths, which can be explained by the fact that all mutations being penalized it may become more advantageous to create a new transcript than to pair up two transcripts. However, the dating of the main ASEs did not change (Supplementary Table S3, combination 3.0.2 b ). We should stress that in all our simulations, exon 7 is present at the root of the forest, while exon 6 appears afterwards (Supplementary Table S3 ). Altogether, these results validate our choice of parameters and show that our JNK phylogeny is robust to small parameter changes.

been subject to much debate 48;56 . The present work contributes to elaborating strategies to answer it, by crossing sequence analysis and phylogenetic inference with molecular modeling. We report the first joint analysis of the evolution of alternative splicing across several species and of its structural impact on the produced isoforms. The analysis was performed on the JNK family, which represents a high interest for medicinal research and for which a number of human isoforms have been described and biochemically characterized.

Firstly, our results allowed dating an ASE consisting of two mutually exclusive homologous exons (6 and 7 ) in the ancestor common to mammals, amphibians and fishes.

We find that the most ancient of these two exons is exon 7. By characterizing in details the structural dynamics of two human isoforms, JNK1α and JNK1β, bearing one or the other exon, we could emphasize subtle changes associated to this ASE and identify residues that may be responsible for the selectivity of the JNK isoforms toward their substrates. Secondly, we highlighted an isoform that was not previously described in the literature, namely JNK1δ. Despite displaying a large deletion (about 80 residues), it is conserved across several species and short MD simulations suggest that it is stable in solution. According to the APPRIS database v20 57 , there are 4 peptides matching this isoform in publicly available proteomics data. By comparison, the other human JNK1

isoforms with a phylogeny have between 5 and 7 matching peptides, while the orphan transcripts identified by our analysis have between zero and 2 matching peptides, suggesting that JNK1δ is indeed translated and stable in solution. Hence, considering that the catalytic site is intact in JNK1δ, we propose that this isoform might be catalytically competent and that the large amplitude motion of the A-loop observed in the simulations might facilitate the activation of the protein by exposing a couple of tyrosine and threonine residues that are targeted by MAPK kinases. The validation of this hypothesis would require further calculations and experiments that fall beyond the scope of this study. Already, this interesting result suggests that our approach could be used to identify and characterize new isoforms, that may play a role in the cell and thus serve as therapeutic targets. Thirdly, we found characteristics specific to the JNK3 isoforms, namely the ab-sence of exon 7 and the presence of twos exons (0 and 1') containing regions predicted to be disordered and involved in interactions. These observations suggest specific competences or functions for this gene. Studies investigating the gain/loss of alternative splice forms associated to gene duplication at large scale 1;58 have highlighted a wide diversity of cases and have suggested that it depends on the specific cellular context of each gene.

Although we did not have a sufficient sample resolution to confirm it with RNA-Seq data, JNK3 is reported to be specifically expressed in the heart brain and testes [START_REF] Waetzig | Context-specific inhibition of JNKs: overcoming the dilemma of protection and damage[END_REF] Our approach enables to go beyond a description of transcript variability across species and/or across genes. Indeed, by reconstructing phylogenies, we do not only cluster transcripts but we also add a temporal dimension to the analysis. Previous methods reported in the literature were only applied on simple cases [START_REF] Christinat | Inferring transcript phylogenies[END_REF] and/or largely simplified the evolution model to increase computational efficiency [START_REF] Christinat | A transcript perspective on evolution[END_REF] , such that they could not be used for ancestral transcripts' reconstruction. PhyloSofS algorithm makes the reconstruction of transcripts' phylogenies feasible for any gene family. For the JNK family, the execution of 1 million iterations took about two weeks on a single CPU. This case represents a high level of complexity as most of the transcripts contain more than 10 exons (the average number of exons per gene being estimated at 8.8 in the human genome [START_REF] Sakharkar | Distributions of exons and introns in the human genome[END_REF] ) and up to 8 transcripts are observed within each species (it is estimated that about 4 distinct-coding transcripts per gene are expressed in human 57 ). To reduce the computing time, the user can easily parallelize the multi-start iterative search on multiple cores and s/he has the possibility to give as input a previously computed value for the lower bound (to increase the efficiency of the cut). We should stress that the problem of pairing transcripts across homologous and paralogous genes between different species, addressed here, is much more complex than that of inferring the transcripts' phylogeny of each gene separately. Indeed, in the former case, the problem size is bigger, one needs to reconcile the gene tree with the species tree, and the sequences are more divergent.

Our phylogenies may be impacted by two main sources of error coming from the input data. Specifically, under-annotation of transcripts can lead to missing distant evolutionary relationships. To deal with this issue, we set the cost associated to transcript death to zero. This enables to construct trees that can relate transcripts possibly very far from each other in the phylogeny (i.e. expressed in very distant species, because some species in between are under-annotated). This parameter may be tuned by the user depending on the quality and reliability of the input data. A second source of error comes from annotated transcripts that are not translated or not functional at the protein level. However, we do not expect that these transcripts will significantly pollute the phylogenetic reconstruction. Indeed, they are likely not conserved across species and thus will be attributed the status of orphans in the phylogenetic reconstruction. Moreover, we have emphasized an independent source of evidence coming from their structural characterization which can help us flag them. The reliability of the transcript expression data clearly constitutes a present limitation of the method. However, as experimental evidence accumulate and precise quantitative data become available, computational methods such as PhyloSofS will become instrumental in assessing the contribution of AS in protein evolution.

Although PhyloSofS was applied here to study the evolution of transcripts in different species, it has broad applicability and can be used to study transcript diversity and conservation among diverse biological entities. The entities could be at the scale of (i) one individual/species (tissue/cell differentiation), (ii) different species (matching cell types), (iii) population of individuals affected or not by a multifactorial disorder. In the first case, the tree given as input should describe checkpoints during cell differentiation and PhyloSofS will provide insights on the ASEs occurring along this process. In the second case, PhyloSofS can be applied to study one particular tissue across several species in a straightforward manner (explicitly dealing with the dimension of different tissues requires further development). In the third case, the tree given as input may be constructed based on genome comparison, a biological trait or disease symptoms. PhyloSofS can be used to evaluate the pertinence of such criteria to relate the patients, with regards to the likelihood (parsimony) of the associated transcripts scenarios. This case is particularly relevant in the context of medical research.

Materials and Methods

PhyloSofS workflow

PhyloSofS can be applied to single genes or to gene families. The input is a binary tree (called a gene tree) describing the phylogeny of the gene(s) of interest for a set of species (Fig. 1a, on the left), and the ensemble of transcripts observed in these species (symbols at the leaves). Only transcripts annotated as coding for a protein are considered.

PhyloSofS comprises two main steps:

• It reconstructs a forest of phylogenetic trees describing plausible evolutionary scenarios that can explain the observed transcripts (Fig. 1a, on the right). The forest is embedded in the input gene tree and is reconstructed by using the maximum parsimony principle. The root of a tree corresponds to the creation of a new transcript, each leaf stands for an observed transcript and a dead end (indicated by a triangle on Fig. 1a, on the right) indicates a transcript loss. Transcripts can mutate along the branches of the trees.

• It predicts the three-dimensional structures of the protein isoforms corresponding to the observed transcripts by using homology modeling. The 3D models are then annotated with quality measures and with exon labels.

PhyloSofS comes with helper functions for the visualization of the output transcripts' phylogeny(ies) and of the isoforms' molecular models. The program is implemented in Python 3 and freely available at GitHub under MIT license: https://github.com/PhyloSofS-Team/PhyloSofS.

Step a. Transcripts' phylogenies reconstruction

For simplicity and without loss of generality, we describe here the case of one gene of interest studied across several species. The gene is represented by an ensemble E of n e exons. The identification and alignment of the n e homologous exons between the different transcripts must be performed prior to the application of the method (see below for details on data preprocessing for the JNK family). The n s transcripts of species s are described by a binary table T s of n e × n s elements, where T s i,j = 1 if exon i is included in transcript j (Supplementary Fig. S1a, see colored squares), 0 if it is excluded (white squares).

We model transcripts evolution as a two-level process, at the gene and transcript levels, as described by Christinat and Moret [START_REF] Christinat | Inferring transcript phylogenies[END_REF] . At the level of the gene, each exon can be either absent, alternative or constitutive. This status is inferred from the occurrence of the exon in the transcripts. Hence, for a given species s, a vector g s of length n e encodes the state of each exon by the values {0, 1, 2} for absent, alternative and constitutive, respectively (Supplementary Fig. S1b, white, black/white and black squares). At the leaves (current species), the components of g s are calculated as:

g s i = ns j=1 T s i,j + 1 - ns j=1 (1 -T s i,j ) (1)
As in [START_REF] Christinat | Inferring transcript phylogenies[END_REF] , the g s vectors for internal nodes (ancestral species) are determined by using Sankoff's algorithm 62 . Dollo's parsimony principle is also respected, such that an exon cannot be created twice 2 . If different exon states have equal cost, we follow the priority rule 1 > 0 > 2.

Three evolutionary events are considered, namely creation, death and mutation of a transcript with costs C B , C D and σ, respectively. The mutation cost σ is accounted for only when the associated evolutionary change occur at the level of the transcript (Table II). This reflects the fact that changes at the level of the gene affects the expression of exons in the transcripts but changes at the level of the transcripts do not affect the gene structure. For instance, if an exon is absent in a parent and becomes present in the child, then this change of status at the transcript level will be penalized by σ only if the exon could be absent in the child, i.e. its status at the gene level is "alternative". If the "constitutive" exon is the child, then the mutation is not penalized (Table II, compare the cells (0,0)→(1,1) and (0,0)→(1,2).

Each internal node of the gene tree, representing an ancestral species, is expanded in several subnodes, representing the transcripts of the gene in this ancestral species (Supplementary Fig. S1c). There exist three types of subnodes: binary (two transcript and C death (S) are the total costs of creation and loss of transcripts, expressed as

C birth (S) = C B × |S| (2)
C death (S) = C D × nodes N n l (N ) + n r (N ), (3)
where |S| is the number of trees in the forest.

Given a forest structure, a transcripts' phylogeny determines the pairings of transcripts at each internal node (Supplementary Fig. S1d). The cost of the transcripts' phylogeny ϕ complying with the forest structure S is calculated as:

C ϕ = C S + A tree of ϕ Γ(A) (4)
where Γ(A) is computed for each tree A of ϕ by evaluating the changes of exon states along the branches of ϕ:

Γ(A) = t k i →t l j branch of A Γ(t k i → t l j ) (5)
where t k i is the parent transcript, i th subnode of node k, t l j is the child transcript, j th subnode of node l and Γ(t k i → t l j ) = e∈E σ((T k e,i ; g k e ), (T l e,j ; g l e )), with g y e ∈ {0, 1, 2} the state of exon e at the level of the gene at node y and T y e,x ∈ {0, 1} the state of exon e at the level of the x th transcript of node y. The evolution costs σ are given in Table II.

PhyloSofS algorithm seeks to determine the scenario with the smallest number of evolutionary events, i.e. the transcripts' phylogeny with the minimum cost (Supplementary Fig. S1c-d). It proceeds as follows:

Initialization:

C min ← ∞
Choose the forest structure S 0 that maximizes the n b values Iteration:

for i = 0 to t max -1 do if C S i < C min then
Find the most parsimonious phylogeny ϕ i given structure

S i if C ϕ i < C min then C min ← C ϕ i end if end if
Choose forest structure S i+1 by setting n b , n l and n r at every internal node

end for

To efficiently search the space of all possible forest structures (Supplementary Fig. S1c), PhyloSofS relies on a multi-start iterative procedure. Random jumps in the search space are performed until a suitable forest structure S i (with C S i < C min ) is found. The cost C S i of the forest structure S i serves as a lower bound for the cost C ϕ i of the phylogeny

Step b. Isoforms structures prediction

The molecular modeling routine implemented in PhyloSofS relies on homology modeling.

It takes as input an ensemble of multi-fasta files (one per species) containing the sequences of the splicing isoforms. For each isoform, it proceeds as follows:

1. search for homologous sequences whose 3D structures are available in the PDB (templates) and align them to the query sequence;

2. select the n (5 by default, adjustable by the user) best templates;

3. build the 3D model of the query;

4. remove the N-and C-terminal residues unresolved in the model (no structural template);

5. annotate the model with sequence and structure information.

Step 1 makes extensive use of the HH-suite [START_REF] Hildebrand | Fast and accurate automatic structure prediction with HHpred[END_REF] and can be decomposed in:

Retrieval and pre-processing of JNK annotated transcriptome data

The peptide sequences of all splice variants from the JNK family observed in human, mouse, Xenopus tropicalis, zebrafish, fugu, Drosophila melanogaster and nematode were retrieved from Ensembl 75 release 84 (March 2016) along with the phylogenetic gene tree.

Only the transcripts containing an open reading frame and not annotated as undergoing nonsense mediated decay or lacking 3' or 5' truncation were retained.The isoforms sharing the same amino acid sequence were merged. The homologous exons between the different genes in the different species were identified by aligning the sequences with MAFFT [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: improvements in performance and usability[END_REF] , and projecting the alignment on the human annotation. They do not necessarily represent exons definition based on the genomic sequence and this can be explained by two reasons.

First, the gene structure may be different from one species to another. For instance, the third and fourth exons of human JNK1 genes are completely covered by a single exon in the Drosophila melanogaster JNK gene (Supplementary Fig. S4 ). In cases like this, we keep the highest level of resolution and define two exons (e.g. numbered 3 and 4 ).

Secondly, it may happen that a transcript contains only a part of an exon in a given species translated in another frame. In that case, we define two exons sharing the same number but distinguished by the prime symbol (e.g. exons 8 and 8'). In total, 64 transcripts comprised of 38 exons were given as input to PhyloSofS.

PhyloSofS' parameter setting

To set the parameters, two criteria were taken into consideration. First, the different genomes available in Ensembl are not annotated with the same accuracy and the transcriptome data and annotations may be incomplete. This may challenge the reconstruction of transcripts' phylogenies across species. To cope with this issue, we chose not to penalize transcript death (C D =0). Second, the JNK genes are highly conserved across the seven studied species (Table I), indicating that this family has not diverged much through evolution. Consequently, we set the transcript mutation and birth costs to σ = 2

and

C B = 3 (C B < σ × 2)
. This implies that few mutations will be tolerated along a phylogeny. Prior to the phylogenetic reconstruction, PhyloSofS removed 19 exons that appeared in only one transcript (default option), reducing the number of transcripts to 60. This pruning enables to limit the noise contained in the input data and to more efficiently reconstruct phylogenies. PhyloSofS algorithm was then run for 10 6 iterations.

The 3D models of all observed isoforms were generated by PhyloSofS molecular modeling routine by setting the number of retained best templates to 5 (default parameter) for every isoform.

Analysis of JNK tertiary structures.

The list of experimental structures deposited in the PDB for the human JNKs was retrieved from UniProt 5 . The structures were aligned with PyMOL [START_REF] Delano | The PyMOL Molecular Graphics System[END_REF] and the RMSD between each pair was computed. Residues comprising the catalytic site were defined from the complex between human JNK3 and adenosine mono-phosphate (PDB code: 4KKE, resolution: 2.2 Å), as those located less than 6 Å away from the ligand. Residues comprising the D-site and the F-site were defined from the complexes between human JNK1 and the scaffolding protein JIP-1 (PDB code: 1UKH, resolution: 2.35 Å28 ) and the catalytic domain of MKP7 (PDB code: 4YR8, resolution: 2.4 Å44 ), respectively. They were detected as displaying a change in relative solvent accessibility >1 Å2 upon binding.

The I-TASSER webserver 59;73;76 was used to try and model the regions for which no structural templates could be found. DISOPRED [START_REF] Ward | The DISOPRED server for the prediction of protein disorder[END_REF] and IUPred 21 were used to predict intrinsic disorder. JET2 40 was used to predict binding sites at the surface of the isoforms.

Molecular dynamics simulations of human isoforms.

The 3D coordinates of the human JNK1 isoforms JNK1α (369 res., containing exon 6 ), JNK1β (369 res., containing exon 7 ) and JNK1δ (304 res., containing neither exon 6 nor exon 7 ) were predicted by PhyloSofS pipeline. The 3 systems were prepared with the LEAP module of AMBER 12 [START_REF] Case | Amber 12. University of California[END_REF] , using the ff12SB forcefield parameter set: (i) hydrogen atoms were added, (ii) the protein was hydrated with a cuboid box of explicit TIP3P

water molecules with a buffering distance up to 10 Å, (iii) Na + and Cl -counter-ions were added to neutralize the protein.

The systems were minimized, thermalized and equilibrated using the SANDER module of AMBER 12. The following minimization procedure was applied: (i) 10,000 steps of minimization of the water molecules keeping protein atoms fixed, (ii) 10,000 steps of minimization keeping only protein backbone fixed to allow protein side chains to relax, (iii) 10,000 steps of minimization without any constraint on the system. Heating of the system to the target temperature of 310 K was performed at constant volume using the Berendsen thermostat 6 and while restraining the solute C α atoms with a force constant of 10 kcal/mol/ Å2 . Thereafter, the system was equilibrated for 100 ps at constant volume (NVT) and for further 100 ps using a Langevin piston (NPT) [START_REF] Loncharich | Langevin dynamics of peptides: The frictional dependence of isomerization rates of n-acetylalanyl-N'-methylamide[END_REF] to maintain the pressure.

Finally the restraints were removed and the system was equilibrated for a final 100 ps run.

Each system was simulated during 250 ns (5 replicates of 50 ns, starting from different initial velocities) in the NPT ensemble using the PMEMD module of AMBER 12. The temperature was kept at 310 K and pressure at 1 bar using the Langevin piston coupling algorithm. The SHAKE algorithm was used to freeze bonds involving hydrogen atoms, allowing for an integration time step of 2.0 fs. The Particle Mesh Ewald (PME) method [START_REF] Darden | Particle mesh ewald: An nlog(n) method for ewald sums in large systems[END_REF] was employed to treat long-range electrostatics. The coordinates of the system were written every ps.

Standard analyses of the MD trajectories were performed with the ptraj module of AMBER 12. The calculation of the root mean square deviation (RMSD) over all atoms indicated that it took between 5 and 20 ns for the systems to relax. Consequently, the last 30 ns of each replicate were retained for further analysis, totaling 150 000 snapshots for each system. The fluctuations of the C-α atoms were recorded along each replicate.

For each residue or each system, we report the value averaged over the 5 replicates and the standard deviation (see Fig. 4a). The secondary structures were assigned by DSSP algorithm over the whole conformational ensembles. For each residue, the most frequent secondary structure type was retained (see Fig. 4a and Supplementary Fig. S9 ). If no secondary structure was present in more than 50% of the MD conformations, then the residue was assigned to a loop. The amplitude of the motion of the A-loop compared to the rest of the protein was estimated by computing the angle between the geometric center of residues 189-192, residue 205 and either residue 211 in the isoforms JNK1α and JNK1β or residue 209 in the isoform JNK1δ. Only C-α atoms were considered.

RNA-Seq Data integration

To obtain additional support for transcript isoform expression, we queried the Bgee database (v.14) 4 for a list of all RNA-Seq experiments related to the selected species.

Using SRA tools, we downloaded raw sequences from H. sapiens (224 samples), M. musculus (155 samples), Xenopus tropicalis (69 samples) and D. rerio (67 samples) and then aligned the reads using STAR v.2.5.3a 20 with default parameters. T. rubripes is not annotated in Bgee and was not integrated in this part of the analysis. Reads overlapping exon-exon boundaries (e.g. splice-junction reads) next to alternative splicing events provide direct evidence for the expression of specific transcripts isoforms. Combined with sample annotation, they could also inform on tissue specific isoform expression. We thus considered all reads included within one of the JNK genes and monitored the alignment of splice junctions between different exons as support for the transcripts isoforms: exons 5-6 and 6-8 for JN1α, exons 5-7 and 7-8 for JNK1β, and exons 5-8' for JNK1δ.

Figure 1 :

 1 Figure1: Transcripts' phylogenies reconstructed by PhyloSofS. (a) On the left, example of a phylogenetic gene tree where 8 transcripts (represented by geometrical symbols) are observed in 4 current species (leaves of the tree, colored in different grey tones). These data are given as input to PhyloSofS. In the middle, the problem addressed by PhyloSofS is that of a partial assignment: how to pair transcripts so as to maximize their similarity? On the right, example of a solution determined by PhyloSofS. The transcripts' phylogeny is a forest comprised of 3 trees (colored differently). The nodes of the input gene tree are subdivided into subnodes corresponding to observed (current) or reconstructed (ancestral) transcripts. The root of a tree stands for the creation of a new transcript and is associated to a cost C B . Triangles indicate transcript deaths and are associated to a cost C D . Mutation events occur along branches and are associated to a cost σ. The grey node corresponds to an orphan transcript for which no phylogeny could be reconstructed. (b) Transcripts' phylogeny reconstructed by PhyloSofS for the JNK family. The forest is comprised of 7 trees, 19 deaths (triangles) and 14 orphan transcripts (in grey). Mutation events are indicated on branches by the symbol + orfollowed by the number of the exon being included or excluded (e.g. +11 ). The cost of the phylogeny is 69 (with C B = 3, C D = 0 and σ = 2). On the top right corner are displayed the exons present in each current species (in black). On the bottom left corner are displayed the exon compositions of the human isoforms for which a phylogeny could be reconstructed.

Figure 2 :

 2 Figure 2: Exons mapped onto the tertiary structure of human JNK1. The protein (residues 7 to 364) is represented as a cartoon and the different exons are colored from blue through white to red. The residues in yellow are at the junction of 2 exons.It should be noted that exons 8 and 8' used in PhyloSofS actually correspond to only one genomic exon (see Methods). The regions labelled in black are common to kinases and were reported in the literature (see[START_REF] Huse | The conformational plasticity of protein kinases[END_REF] ) for playing important roles in their structural stability and/or function. The regions labelled in grey are specific to MAP kinases. The green circles indicate the catalytic site and binding sites for JNK cellular partners28;44 . The structure was solved by X-ray crystallography at 1.80 Å resolution (PDB code: 3ELJ[START_REF] Chamberlain | Optimization of 4,6bis-anilino-1H-pyrrolo[2,3-d]pyrimidine IGF-1R tyrosine kinase inhibitors towards JNK selectivity[END_REF] ).

Figure 3 :

 3 Figure 3: Structural features of the transcript isoforms. Distributions are reported for the parented transcripts (in light gray) and the orphan transcripts (in dark grey) in the transcripts' phylogeny (see Fig. 1b). (a) Length of the transcript (in residues). (b) Predicted secondary structure content (in percentages of residues). (c) Overall G-factor computed by Procheck 42 . (d) Normalized DOPE score computed by Modeller 47 . (e) Fraction of protein residues being exposed to the solvent (rsa > 0.25). (f) Fraction of hydrophobic protein residues being exposed to the solvent (rsa > 0.25).

5 Figure 4 :

 54 Figure 4: Dynamical behavior of the human JNK1 isoforms in solution. (a) The secondary structures for JNK1α (with exon 6 ) are depicted on top (the profiles for the 2 other isoforms are very similar, see Supplementary Fig. S9 ). The atomic fluctuations (computed on the Cα) averaged over 5 50-ns MD replicates are reported for JNK1α in orange, JNK1β in purple and JNK1δ in pink. The envelopes around the curves indicate the standard deviation. (b) Representative MD conformations obtained by clustering based on position 228 (RMSD cutoff of 1.5 Å). There are 8 conformations for JNK1α (in orange) and only 1 for JNK1β (in purple). (c) Superimposed pair of MD conformations illustrating the amplitude of the A-loop motion in JNK1δ (see Materials ad Methods for details on the calculation of the angle). Exons 5, 8' and 9 are indicated by colors and labels. For clarity, 8' is also indicated by two stars on the structure.

  children), left (one transcript child in the node's left child) and right (one transcript child in the node's right child). Left and right subnodes imply that a transcript death occurred along the branch. A forest structure S is fixed by setting n b , n l and n r the respective numbers of binary, left and right subnodes for every internal node of the gene tree. The cost associated to structure S is calculated as C S = C birth (S) + C death (S), where C birth (S)

  (a) search for homologous sequences and building of a multiple sequence alignment (MSA), by using HHblits 55 , (b) addition of secondary structure predictions, obtained by PSIPRED 34 , to the MSA, (c) generation of a profile hidden markov model (HMM) from the MSA, (d) search of a database of profile HMMs for homologous proteins, using HHsearch 63 . Step 3 is performed by Modeller 47 with default options. Step 5 consists in: (a) inserting the numbers of the exons in the β-factor column of the PDB file of the 3D model, (b) computing the proportion of residues predicted in well-defined secondary structures by PSIPRED 34 , (c) assessing the quality of the model with Procheck 42 and with the normalized DOPE score from Modeller, (d) determining the by-residue solvent accessible surface areas with Naccess 31 and computing the proportions of surface residues and of hydrophobic surface residues.

Table I : Percentages of sequence identity of JNK genes to human.

 I

		JNK1		JNK2 JNK3
	Mouse	99		97	100
	Xenopus tropicalis	89		-	98
	Fugu	79		81	96
		82	(a)
	Zebrafish	87	(a)	85	93
		87	(b)
	Drosophila melanogaster	78		73	77
	Nematode	56		54	56
	Each gene of each species was aligned to its orthologous gene in human. Human and mouse
	genomes contain 3 paralogues: JNK1, JNK2 and JNK3. Xenopus tropicalis contains only
	JNK1 and JNK3. The fishes contain 4 paralogues: JNK1, JNK1a, JNK2 and JNK3 in fugu,
	JNK1a, JNK1b, JNK2 and JNK3 in zebrafish. Drosophila melanogaster and nematode contain
	only one gene, whose sequence identities with human JNK1, JNK2 and JNK3 are displayed in
	black, grey and grey, respectively. In addition to the values reported in the table, here are
	some sequence identities computed between paralogues: (i) 83% between human JNK1 and
	JNK2, and between human JNK1 and JNK3; (ii) 86% between fugu JNK1 and JNK1a; (iii)
	92% between zebrafish JNK1a and JNK1b.

  The numbers of trees, deaths and orphans are given in y-axis for seven combinations of parameters. The parameter values are indicated on the x-axis, in the following order: C B .C D .σ. The distributions are computed over 1756, 4204, 2880, 576, 47, 326 and 1100 solutions of equal costs, found during 10 6 iterations of the algorithm.

		12		25
		11		20		20
	#(trees)	8 9 10	#(deaths)	10 15	#(orphans)	16 18
		7		5		14
		6				12
		302 503 504 502 532 332 222 3.0.2 5.0.3 5.0.4 5.0.2 5.3.2 3.3.2 2.2.2		302 503 504 502 532 332 222 3.0.2 5.0.3 5.0.4 5.0.2 5.3.2 3.3.2 2.2.2		302 503 504 502 532 332 222 3.0.2 5.0.3 5.0.4 5.0.2 5.3.2 3.3.2 2.2.2
	Figure 5: Statistical analysis of the phylogenies obtained with different parame-
	ters.

5, combinations 5.3.2, 3.3.2 and 2.2.2), while the resulting scenarios were less parsimonious (Supplementary Table

S3

). For example, the JNK1δ transcript was created twice, at the internal nodes exons 6/7 and 12/13 also appeared multiple times in independent sub-forests (Supplementary Table

Table II :

 II Exon states and associated costs σ.Each cell is associated to the evolution of the state of the exon e from a parent transcript to a child transcript. The first and second terms in parenthesis correspond to the status of the exon at the transcript and gene levels, respectively. Only evolutionary changes taking place at the transcript level, without being directly induced by a gene-level change, are penalized (σ). Zeros highlighted in bold indicate transcript-level changes being a direct consequence of a gene-level change. This table of costs was taken from[START_REF] Christinat | Inferring transcript phylogenies[END_REF] .

	child / parent (0,0) (0,1) (1,1) (1,2)
	(0,0)	0	0	0	0
	(0,1)	0	0	σ	σ
	(1,1)	σ	σ	0	0
	(1,2)	0	σ	0	0

Acknowledgments A grant of the French national research agency (MASSIV project, ANR-17-CE12-0009) provided a salary to D.J.Z. and A.L.. We thank Y. Christinat for providing information on the algorithm he developed for the reconstruction of transcript phylogenies.

Complexity and comparison with other methods

The size of the search space for the transcripts' phylogeny reconstruction grows exponentially with the number of observed transcripts (leaves). To explore that space, the heuristic algorithm implemented in PhyloSofS relies on a multi-start iterative procedure and on the computation of a lower bound to early filter out unlikely scenarios (see Methods).

As the algorithm finds better and better solutions, the filtering procedure becomes more and more efficient. Transcripts assignment (Fig. 1a, in the middle) is performed only when the computed lower bound passes the filter. At each internal node, k (among n l ) transcripts from the left child must be paired with k (among n r ) transcripts from the right child. We solve this problem with a branch-and-bound algorithm whose complexity is of O(n 3 ) for n l = n r = n (see details in Supplementary Text S1 ). Hence, it requires about (s-1)n 3 operations to pair up transcripts from the leaves to the root, with s the number of current species (s-1 being the number of ancestral species). Another heuristic algorithm, making use of a neighbor-joining operation, has been proposed in the literature [START_REF] Christinat | A transcript perspective on evolution[END_REF] . This algorithm requires to look at up to (s-1)n 2

possible pairs of transcripts at the level just below the root, and recursively applies a neighbor-joining procedure to each one of these pairs down to the leaves. One advantage of this algorithm is that it constrains the space of considered phylogenies and thus may in practice be more efficient than PhyloSofS. However, it does not consider transcript mutations and hence cannot reconstruct ancestral transcripts (see Fig. 1A in [START_REF] Christinat | A transcript perspective on evolution[END_REF] ). Moreover, the information about transcript ancestry is virtually lost as trees can be merged by the neighbor-joining operation (see Fig. 1B in [START_REF] Christinat | A transcript perspective on evolution[END_REF] ). As a consequence, this algorithm is not suitable for inferring evolutionary scenarios explaining observed transcripts nor for dating AS events.

Discussion

To what extent the transcript diversity generated by AS translates at the protein level ϕ i . Forest structures that are too costly are simply discarded, without calculating the corresponding phylogenies. As the algorithm finds better and better solutions, the cut becomes more and more efficient. The phylogeny ϕ i is reconstructed by using dynamic programming. Sankoff's algorithm is applied bottom up to compute the minimum pairing costs between transcripts (Supplementary Fig. S1d, each transcript is represented by a matrix of costs). At each internal node, the pairings are determined by using a specific version of the branch-and-bound algorithm [START_REF] Land | An automatic method of solving discrete programming problems[END_REF] (see Supplementary Text S1 ). If the reconstructed phylogeny is more parsimonious than those previously visited (C ϕ i < C min ), then the minimum cost C min is updated. There may be more than one phylogeny with minimum cost that comply with a given structure S i . The next forest structure S j will be randomly chosen among the immediate neighbors of S i (Supplementary Fig. S1d). Two structures are immediate neighbors if each one of them can be obtained by an elementary operation applied to only one node of the other one (Supplementary Fig. S13 ). If the phylogeny ϕ j is such that C ϕ j < C min , then the next forest structure will be chosen among the neighbors of S j , which serves as a new "base" for the search. Otherwise, the algorithm continues to sample the neighborhood of S i . This step-by-step search is applied until no better solution can be found. At this point, a new random jump is performed.

The total number of iterations t max is given as input by the user (1 by default).

We should stress that PhyloSofS algorithm is designed to deal with much more complex cases than those reported in [START_REF] Christinat | Inferring transcript phylogenies[END_REF] in a computationally tractable way. Hence, it differs from the algorithm reported in 16 in several respects. First, our multi-start iterative strategy relies on random jumps in the forest structure space combined with systematic local exploration around the best current solution, while Christinat and Moret [START_REF] Christinat | Inferring transcript phylogenies[END_REF] proposed an exhaustive generation and evaluation of forest structures. Secondly, we have designed a branch-and-bound algorithm specifically adapted to the problem of determining the best phylogeny complying with a given forest structure (see Supplementary Text S1 ). Both aspects contribute to PhyloSofS efficiency in reconstructing transcripts' phylogenies.

PhyloSofS generates PDF files displaying the computed transcripts' phylogenies using a Python driver to the Graphviz 23 DOT format.