dataesr/scientific-paragraphs-categorization huggingface dataset:
Photosynthesis Functional annotation with eggNOG revealed that many genes related to Energy production and conversion (group C, Figure  4C ) were differentially methylated in several groups (K2, K6, K7 and K9 especially). Behind this wide-scale annotation, functions related to photosynthesis were found. Our dataset indeed suggests a strong control of the photosynthesis process by DNA methylation. GO term enrichment (Figures  S4  and S5 ) clearly revealed a significant over-representation of genes related to photosynthesis in comparisons involving roots and flowers. We found a good correlation between hypomethylation and increased expression for these genes (Figure  4B  and Table  S2 ). It is therefore likely that the establishment of photosynthesis in green tissues may be triggered by a strong hypomethylation of photosynthesis-related genes.
Defense Response It is well-established that salicylic acid treatments or challenges with phytopathogens trigger strong modification of the DNA methylation state  [61, 62]. In our dataset, the GO term "defense response" was found to be significantly enriched in several comparisons (Figure  S4 ). For instance, a large reprogramming in the plant defense system characterized the Yleafcat sample (Figure  4D ). These leaves were collected from plants one week after having been fed by caterpillars. In this organ, an important increase of MIA has previously been found in contrast to control leaves, suggesting a systemic activation of the plant defense system. DEMs genes in this Yleafcat sample particularly corresponded to transcription factors and kinases (Figure  4D  and Table  S2 ). Interestingly, we found an ortholog (CRO_T110712) of the A.
The resulting DNA library was purified and bisulfite converted. A qPCR assay was used to determine the optimal number of PCR amplification cycles (between 10 and 15 cycles) required to obtain a high diversity library with minimal duplicated reads prior to final library amplification. The sequencing was performed with paired ends (2 × 150 bp) on an Illumina HiSeq4000 platform. Raw data were stored in FASTQ files (SRA record is under the reference PRJNA635601, see Table  2 ).
Epigenomics Data Generation and Analysis The bioinformatics pipeline used in this study is adapted from the ENCODE pipeline (https: //www. encodeproject. org/wgbs/) and installed in a Galaxy framework (http://bioinfo. univ-perp. fr/, Perpignan, France). First, quality control and cleaning of the raw data were carried out by only considering nucleotides with a quality score over 26, and reads which had more than 95% of their nucleotides over this quality threshold. The second step was the alignment of the WGBS reads on the latest available genome sequence (https://doi. org/10. 5061/dryad. 08vv50n) [3]  by using BISMARK (version 0. 16. 3,  [77]  and Bowtie 2 tools (version 2. 1. 0; [78] ). Due to sequenced genome quality (2090 scaffolds in version 2, N residues in intergenic sequences), BSMAP version 2.
74 [79,80]  was used to measure DNA methylation levels at each cytosine site on a reduced representation consisting on gene models (34,363) including 5 Kbp flanking regions. All data were processed in R  [81]. Multidimensional scaling was performed using only sites covered by at least seven reads were kept as recommended in the EdgeR differential methylation detection procedure  [82]. Differential methylation in pairwise sample comparisons was detected following the DSS procedure  [83]. It models the biological variation using neighboring C sites as pseudo-replicates and automatically adjusts library coverage differences. DMC were kept when the differential methylation p-value was below 0. 05 (gene-wise Wald test on two-group comparisons). DMR were detected when at least 3 DMC were present in the same region of at least 50 bp.
Transcriptome Analysis Gene expression was obtained from previously published RNA-seq data as well as new samples (flowers, basal, and distal parts of young leaves, deposited as project accession number PRJEB38854, see Table  2 ). RNA were obtained from plants as described above using a Trizol-based extraction  [14]. Libraries prepared with the Illumina TruSeq RNA Library Prep Kit v2 were sequenced by Eurofins Genomics using the Illumina HiSeq2000/2500 technology, in the paired-end mode (2 × 100 pb). Sequencing runs yield a total of 25,352,619, 20,959,861, and 17,313,293 reads for samples ERR4235477, ERR4235478, and ERR4235479, respectively. Reads were first trimmed with FastP v0. 20  [84]  using default parameters and next pseudo-aligned on transcripts predicted from the second C. roseus genome version  [3]  and quantification was performed with Salmon v0.
Statistical Procedures All analyses were made within the R programming language  [81]. Tables were all used in the data. table format  [87]. All graphs were made with the ggplot2 package  [88] , together with the ggdendro  [89]  and ggrepel  [90]  extensions. The PCA in Figure  2  was constructed with the FactoMineR package  [91]. Correlations in Figure  3C ,D were calculated with the Pearson coefficient and tested with Pearson's product moment correlation coefficient. Heatmaps were plotted with the pheatmap package  [92]. Hierarchical trees were constructed with the Ward agglomeration criterion on Euclidean distances. Groups of DEM in Figure  4  were detected using a dynamic tree cutting algorithm  [56]  with a minimum cluster size of 30. Functional analysis of genes was performed by transferring annotations from Uniprot homologs (blastx,  [93] ) or using the eggNOG mapper  [94].
Functional term (GO or Uniprot keywords) enrichment was tested using a hypergeometric distribution Figures  S4  and S5 ). P-values were adjusted using the Bonferroni correction. The gene co-expression network (Figure  6 ) was constructed using a Spearman's rho-based distance matrix and the igraph package  [95]. Significant correlations were retained if rho > 0. 95. When applicable, differences between treatments were calculated using the Student's t-test or the non-parametric Wilcoxon rank sum test.
Acknowledgments: We also thank  Jean-Louis Rouet  and  Laurent Catherine  for access and support to the CCSC computing resources (Cascimodot Federation, CNRS, Orléans). We acknowledge the  GenOuest bioinformatics core facility  (https://www. genouest. org)and the  IHPE Galaxy Perpignan  (https://bioinfo. univ-perp. fr/) for providing the computing infrastructure. Hao Wu  is acknowledged for his help with the DSS package,  Pr. Elisabeth Huguet  and  Karine Musset  for the help with the caterpillar experiments and  Florent Lafontaine  for the technical help with DNA extraction. We thank  Anne-Laure Le Gac ,  Régis Fichot ,  Christoph Grunau , and the collaborative GDR3E network (https://rtp-3e. wixsite. com/rt3e)for helpful discussions concerning methylome analysis.
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  37  and likely suffer from strong experimental detection bias  69. A recent analysis of ribosome profiling data suggested that a major fraction of splice variants is translated, with direct implications on specific cellular functions  72. 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  74. From a structural perspective, very few alternatively spliced isoforms have been characterized and are available in the Protein Data Bank (PDB)  7;27.
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.
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  16.
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  29. 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.
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.
Although this can correspond to a real loss of transcripts in those species, a more 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.
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.
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.
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  17. This algorithm requires to look at up to (s-1)n 2 + (s-1)(n-1) 2 + (s-1)(n-1)(n-2) 2 +. +1 O(s 3 n 3 ) 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.
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.
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  60  ) 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.
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.
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 par- • It predicts the three-dimensional structures of the protein isoforms corresponding to the observed transcripts by using homology modeling.
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  35  , 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.
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  19  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.