dataesr/scientific-paragraphs-categorization huggingface dataset:
g. , regulatory elements,  Andolfatto, 2005) ; (ii) we implicitly modeled the effects of selection against migrant genes by including heterogeneous effective migration rates across the genome. Data set 1: "exome capture" We used the dataset already published in  Fraïsse et al. (2016)  and available on http://www. scbi. uma. es/mytilus/index. php. Briefly, a set of 890 EST contigs was used as a reference for a pre-capture multiplex DNA enrichment in samples of eight individuals from two geographical populations in each species (M. edulis: North Sea and Bay of Biscay; M. galloprovincialis: Brittany and Mediterranean Sea, Table  1 ). In addition, we used a sample of four individuals of M. trossulus to serve as an outgroup (Table  1 ). Each DNA library was sequenced twice to increase the per-base coverage (Miseq or GA2X followed by HiSeq2000).
After trimming and quality-filtering, reads of each individual were aligned against the same EST reference sequences using the BWA program (bwa-mem,  Li & Durbin, 2009). Because of the relative divergence between the two species (∼2%, Table  2 ), we adjusted the default parameter of BWA to allow less stringent mapping (minimum seed length k = 10 [default: 19], clipping penalty L = 3 [5], mismatch penalty B = 2  [4] , and gap open penalty O = 3 [6]). Full methods are described in  Fraïsse et al. (2016). We used a maximum-likelihood method, implemented in the program read2snps  (Tsagkogeorga, Cahais & Galtier, 2012; Gayral et al. , 2013) , to call genotypes directly from read numbers at each position. The method computes the probability of each possible genotype after estimating the sequencing error rate.
To limit bias in the site frequency estimation  (Han, Sinsheimer & Novembre, 2013) , a minimum of 10Â coverage, that is, Data set 2: "rna-seq" Data set 2 is made of sequenced transcriptomes (RNA-seq) previously generated in a wide meta-analysis study comparing levels of polymorphism across 76 animal species  (Romiguier et al. , 2014). It includes the transcriptomes of four individuals sampled in the same populations as described above and one individual of M. trossulus (Table  1 ). Briefly, for each individual, cDNA libraries were prepared with total RNA extracted from whole body and sequenced on HiSeq2000. Illumina reads (100 bp, paired-end) were mapped with the BWA program on de novo transcriptomes, independently assembled for each species with a combination of the programs Abyss and Cap3, following the strategy B and D in  Cahais et al.
Data analysis We restricted our analysis to loci assembled in all individuals, and longer than 300 bp after filtering positions containing missing data, or more than two segregating alleles when the three species were aligned. Only synonymous positions were used. The total number of loci and SNPs retained in each dataset is given in Table  2. Technique: rna sequencing ("rna-seq") vs. exome enrichment sequencing ("exome capture"); n: number of individuals analyzed in each species; n locus : total number of locus n SNP : total number of polymorphic sites. The following statistics were calculated for each locus. Their average (in black) and standard deviation (in bold) across all loci are given. S: number of polymorphic sites; Sf: number of fixed differences; Sx: number of exclusive polymorphic sites; Ss: number of shared polymorphic sites; π: number of pairwise differences  (Tajima, 1983) ; u w : Watterson's u  (Watterson, 1975) ; D: Tajima's D  (Tajima, 1989a (Tajima, , 1989b)) ; FST = 1-π S /π T : level of species differentiation, where π S is the average pairwise nucleotide diversity within species and π T is the total pairwise nucleotide diversity of the pooled sample across species; div: total inter-specific divergence; netdiv: net molecular divergence measured at synonymous positions.
Site frequency spectrum We first computed the jSFS for each dataset (Fig. 2  and Fig. S1 ). Each derived allele, oriented by treating the outgroup sequence as a fixed ancestral allele, was assigned to one cell of the jSFS, depending on its frequency in each of the two populations. From the full jSFS, different classes of polymorphism were extracted and used as summary statistics  (Tellier et al. , 2011). Specifically, we used the four Wakeley-Hey classes (jsfs = 4 in Fig. 2 ): fixed differences, Sf; private polymorphisms for each species Sx 1 and Sx 2 ; and shared polymorphisms, Ss  (Wakeley & Hey, 1997). We also considered a summary in which the Sf and Sx classes were split depending on whether the derived allele was fixed or absent in the other species (jsfs = 7 in Fig. 2 ;  Ramos-Onsins et al.
Model choice To choose the best supported model, we followed the methods previously described in  Roux et al. (2013 Roux et al. ( , 2014)). Briefly, posterior probabilities for each of the eleven speciation scenario were estimated with a neural network using the R package abc  (Csilléry, François & Blum, 2012). It implements a non-linear multivariate regression by considering the model itself as an additional parameter to be inferred. The 0. 1% replicate simulations nearest to the observed values of the summary statistics were selected. Moreover, to evaluate the relative probability of the heterogeneous model of migration rates across loci, we compared the alternative models ("homo" vs. "hetero") within each scenario including gene exchange. We approximated Bayes factors (BFs) as the posterior probability of the best supported model divided by that of the model with the highest posterior probability from the remaining candidates.
Model checking We checked the ability of our ABC framework to correctly recover the true model by a "leave-one-out cross validation" from our simulations. We randomly extracted 100 simulated datasets from the million simulations performed for each model. For each of the 100 of datasets simulated under a given model (i. e. , pseudo-observed datasets), we applied the model choice procedure described above to compute the posterior probability of all competitive models. The accuracy rate for model M was calculated as the proportion, among pseudo-observed data inferred to correspond to model M, of those actually generated under model M. The ambiguity rate was computed as the proportion of pseudo-observed data generated under model M whose best model was not strongly supported, that is, its posterior probability was below an arbitrary threshold (P min ) set to be one-third above the expected value given the total number of models compared (1/11 for 11 models, 1/6 for six models, or 1/2 for two models).
Patterns of polymorphism We obtained polymorphism data for two mussels species, M. edulis and M. galloprovincialis, in samples of increasing size (n = 2, n = 4, and n = 8) and for the two datasets ("exome capture" vs. "rna-seq"). The jSFS of each data set, orientated with one outgroup (M. trossulus), is shown in Fig. S1 , while Table  2  gives summary statistics of genetic polymorphism. The data produced by the two sequencing techniques differed in the total number of SNPs (n snp = 3,993 (n = 2); n snp = 5,092 (n = 4); n snp = 5,000 (n = 8) in the "exome capture" data; n snp = 17,275 (n = 2); n snp = 17,902 (n = 4) in the "rna-seq" data). The substantially lower number of SNPs in the "exome capture" data reflects the lower number of loci retained for the analysis due to a reduced and more heterogeneous sequencing depth compared to the "rna-seq" data (while applying identical coverage filters to call SNPs).
g. , BF = P AM. hetero /P AM. homo = 1. 25 with n = 2 and jsfs = 4) compared to the "exome capture" data (BF = P AM. homo /P AM. hetero = 1. 26). This is consistent with the higher heterogeneity of the jSFS in the "rna-seq" data, involving a higher proportion of fixed differences and shared polymorphic sites (Table  2  and Fig. S1 ). Secondly, we evaluated the effect of the number of individuals sampled. As with the number of SNPs, it is clear that sample size had little effect. The best supported scenario remained consistent across the different sampling size (n = 2, 4, or 8 individuals; Table  3A ). For example, when considering 23 classes in the jSFS (jsfs = 23), the heterogeneous PSC scenario was the best supported model whatever the sampling size in both datasets.
Effects on model selection of the number of classes in the jSFS We next investigated the effects of binning the jSFS: (i) on model choice in the mussel datasets (Fig. 2 ; Table  3 ) and (ii) on the ability of the method to discriminate between different speciation scenarios based on simulated datasets (Table  S3 ). Given the limited effects of the sampling strategy, the ABC performance results are presented for a simulated dataset resembling the "exome capture" data with n = 2 individuals only. This allowed us consider the full jSFS when using 23 classes of polymorphism.
Using simulated datasets, we then assessed the performance of the method in identifying the correct model when homogeneous and heterogeneous models were compared (Table  S3A ). The correct model was always recovered (i. e. , an accuracy rate of 1) for the different binnings in all speciation scenarios including gene flow; however, the ambiguity rate did strongly decrease when more information from the jSFS was included. With only four classes (jsfs = 4), none of the replicates showed a posterior probability higher than 0. 83 (the threshold set for the two-model comparison), which corresponds to an ambiguity rate of 1. Similarly, all ambiguity rates were above 0. 97 with seven classes (jsfs = 7). In contrast, when considering the full jSFS (jsfs = 23), we could correctly recover with a strong support the simulated speciation models (e.
(2017)  similarly found that the statistical power of their ABC model selection increases with the number of classes until reaching a plateau of error rates (but above which computation efforts continues to increase). Inferences based on the 23-bin jSFS consistently support a model of PSC with genomic heterogeneity in gene flow (PSC. hetero) in both the "exome capture" and "rna-seq" mussel datasets, and simulations showed that this model has a high accuracy rate, and one of the lowest ambiguity rates (Table  S3 ). Although its relative posterior probability was moderately higher than that of the other models with migration (SC and IM, Table  S1 ), we showed that the method has some power to distinguish PSC from these models (Fig. S3A ). These results are in agreement with previous ABC-based studies revealing clearly a secondary contact history between the two mussel species  (Fraïsse et al.
, 2014; Roux et al. , 2014) , although they relied on eight nuclear loci and on a different set of summary statistics (those that we called "mscalc" here). The periodic connectivity models (PSC and PAM) were not included in these previous studies because of the lack of power to test for intermittent gene exchange since secondary contact. In the present study, we directly compared the use of the jSFS vs. "mscalc" statistics by performing additional ABC inferences on the "mscalc" statistics presented in Table  2. Results were similar to those obtained with jsfs = 23 (Table  S1 ), that is, the best supported models included one or two SCs, except for two datasets ("exome capture" data with n = 2 and n = 8 individuals) for which SI was chosen with a posterior probability of 37% and 34%, respectively.
CONCLUSION In this work, we have shown that two high-throughput sequencing datasets ("exome capture" and "rna-seq"), imply the same history of divergence in mussels, regardless of the number of individuals or SNPs sampled, but conditional on the inclusion of informative classes in the jSFS. Thus, genome-wide data coupled with flexible inference methods allow us to test for more complex scenarios of divergence, by providing a comprehensive picture of the gene histories across the genome. Here, we incorporate in our ABC framework heterogeneity in migration rates among loci, to account for the semi-permeability of the barrier to gene flow between recently diverged species  (Barton & Bengtsson, 1986). This variation in the rate of effective migration results in variation of genetic divergence along the genome.
3 (Thermo Scientific) as interface. Two different workflows were created to handle MSA or ETD data separately. For ETD spectra interpretation, the workflow included a nonfragment filter node that removed a 4-Da window around the selected precursor ion, a 2-Da window around chargereduced precursors and a 2-Da window around known neutral losses up to 120 Da in mass listed in  (39). Database searches were performed with the Mascot server v2. 2. 07 while specifying the following parameters: database TAIR10 (release 2010/12/14, 35386 sequences); enzymatic specificity: tryptic with two allowed missed cleavages; fixed modification of cysteine residues (Methylthio(C)); possible phosphorylation of S, T, and Y residues; 5 ppm tolerance on precursor masses and 0. 6 Da tolerance on fragment ions. Fragment types considered were those specified in the configuration "ESI-trap" in the MSA workflow, and "ETD-trap" in the ETD workflow.
The possible oxidation of methionine residues was ignored, because we verified that its addition did not allow increasing the number of unique phos-phopeptides identified at a fixed false discovery rate (FDR)  (37). The data was then processed using FragMixer  (37) , an in-house developed computational tool that collates the information provided by the database search engine Mascot, from two different spectra with respect to amino acid sequence and location of phosphosite using two simple filtering rules relying on the peptide scores and Mascot Delta scores (MD-scores). It is publically available and can be downloaded from http://proteomics. fr/FragMixer. FragMixer was used to reach an estimated FDR of 1% for all analyses and a false localization rate (FLR) of phosphosites below 5%. We specified MD-score thresholds associated to a FLR below 5% according to  (38).
Relative phosphopeptide abundance was obtained using the program MassChroQ  (40)  which integrates in all the compared runs the areas of the extracted ion chromatograms (XICs) for every phosphopeptide identified in at least one given run. The XIC areas reconstituted by MassChroQ were then further handled as described in the paragraph "Statistical evaluation of relative phosphopeptide quantification. " LC-MS/MS analyses of cytoplasmic samples with the Q-Exactive instrument were interpreted using MaxQuant version 1. 4. 1. 2. Search parameters were as follows: database TAIR10; enzymatic specificity: tryptic with two allowed missed cleavages; possible oxidation of Met residues and acetylation of protein N terminus; 5 ppm tolerance on precursor masses and 20 ppm tolerance on fragment ions. FDR was set to 1% for both peptide-spectrum matches and proteins.
Three biological repeats of cytoplasmic protein extracts were prepared for mpk3, mpk4, and mpk6 mutant plants and four repeats were obtained from WT plants. Within each biological repeat, one sample was treated for 15 min with water and another sample for the same duration with 1 M flg22 (final concentration). Five LC-MS/MS injections were acquired on each phosphopeptide sample using an updated exclusion list. These iterative injections permitted improving quantification, by limiting the risks of missing values because of MassChroQ failing to realign properly the runs to be compared. The lists of phosphopeptides identified and quantified in these cytoplasmic proteins were finally processed as follows. In a first step, we performed WT/mapk pairwise comparisons of the phosphoproteomics data to identify phosphopeptides being completely nondetected in at least one condition.
The digestion was stopped by the addition of 1% formic acid and the peptides were recovered by incubating the gel pieces in ACN. The recovered peptide solution was desalted using C18 ZipTip® (Millipore, Burlington, MA, Cat. No. ZTC18S096) and analyzed by LC-MS/MS. Briefly, peptide samples were separated on a C18 column (Acclaim PepMap C18, 25 cm length ϫ 75 m I. D. ϫ 3 m particle size, 100 Å porosity, Dionex) connected to an LTQ-Orbitrap Velos or a Q-Exactive HF instrument. The LC gradient ramped from 5% solvent B (water/ACN/formic acid, 20/80/0. 1, v/v/v) to 45% solvent B over 45 min, then to 90% solvent B for 10 min. The MS instrument acquired fragmentation spectra on the top 10 peptides using CID fragmentation in the LTQ-Orbitrap or HCD in the Q-Exactive instrument. RAW data files obtained were converted to MGF files using Proteome Discoverer interface (version 1.
4). Database searches were performed with the Mascot server v2. 4 specifying the following parameters: database TAIR10 (release 2010/12/14, 35386 sequences); enzymatic specificity: trypsin permitting two allowed missed cleavages; fixed modification of cysteine residues (Methylthio(C)); possible variable modifications of phosphorylation on S, T, and Y residues; 5 ppm tolerance on precursor masses and 0. 5 Da tolerance on fragment ions. The results were filtered based on Mascot scores and MD-scores. Bimolecular Fluorescence Complementation (BiFC)-To obtain the expression vectors, coding sequences of candidate genes and MAPKs (kindly provided by J. Colcombet) were cloned in fusion with the N-and C-terminal parts of YFP, either as N-or C-terminal fusions, under the control of the cauliflower mosaic virus 35S (CaMV-35S) promoter in the pBIFC1,2,3 and 4 vectors.
Overview of Phosphoproteomes of WT and mapk Mu- tants-To identify phosphoproteins that exhibit changes in phosphorylation levels in any of the mpk3, mpk4 and mpk6 mutants, we performed a phosphoproteomic analysis of at least three independent biological repeats of WT plants and of mpk3, mpk4, and mpk6 mutants after mock or 15 min treatment with flg22 (Fig. 1-A ). Phosphopeptides were enriched from cytoplasmic protein extracts using IMAC and identified by LC-MS/MS. 1,279 different phosphopeptides differing in amino acid sequence and/or number of phosphorylation sites were identified, that corresponded to 731 unique proteins (supplemental Table  S1 , Fig. 1A ). Compared with the list of phosphopeptide sequences present in the database Phos-PhAt (21,646 phosphopeptides), we identified 193 novel phosphopeptides (supplemental Table  S1 ).
When considering the sequences with precise localization of phosphorylation, the total of 1919 phosphosites were identified within 79. 2% of singly, 16. 9% of doubly and 3. 9% of Ն triply phosphorylated peptides (Fig. 1B ). They were additionally distributed into 80. 4% Ser, 16. 0% Thr and 3. 6% Tyr residues (Fig. 1C ). To extract phosphorylation motifs from the dataset, the phosphopeptides identified with precise phosphosites were submitted to motif-x (49) (Fig. 1D , supplemental Table  S1 ). The combined low-stringency (S/T)*P and high-stringency Px(S/T)*P MAPK motifs accounted for the largest peptide population, followed by the casein kinase motif SxxS* and the SnRK/CDPK "simple motif 1" (K/R)xx(S/T)*  (50) , which fits to the involvement of CDPKs in flg22 and biotic stress signaling  (51, 52).
However, the motifs S*DxE, S*ExE, S*xDD, S*xxE, and S*xE correspond to unknown acidophilic kinases. To identify the differentially phosphorylated peptides (supplemental Table  S2 ) that depend on any one of the three MAPKs in the absence of MAMP challenge (Fig. 1E ) or in response to MAMP treatment (Fig. 1F ), we performed WT/ mapk pairwise comparisons for peptides that are absent in any one condition (supplemental Table  S3 ) as well as analyzed the entire dataset by a two-way ANOVA statistical test for relative changes in phosphorylation levels (supplemental Table  S4 ). A total of 152 different phosphorylated sequences (corresponding to 138 proteins) were identified in this way, of which 70 harbored an (S/T)*P site. The differentially phosphorylated peptides were grouped into two categories: those that depend on the genetic background of the three MAPKs under nonimmunity and those that depend on flg22 (immunity-dependent).
Interaction of Immune MAPKs with Their Putative Substrates-The assessment of the quantitative phosphoproteomic data suggested complex kinase-substrate relation- ships among MPK3, MPK4, and MPK6 and several tens of proteins as substrates. To test our conclusions, we selected seven proteins for further assays. These proteins were detected by (S/T)*P-containing phosphopeptides that belonged to the different categories: absent in at least one mapk mutant (Calcium-dependent lipid-binding family protein (CaLB) and Adenine nucleotide alpha hydrolases-like superfamily protein, also called Universal Stress Protein A (USPA)); induced by flg22 (UBX domain-containing protein (PUX9)); mapk genotype-dependent (Target of Myb protein 1 (TOM1); 60S ribosomal protein L13 (BBC1) and UDP-glucose dehydrogenase 2 (UGD2)); and finally, both flg22-induced and mapk genotypedependent (Aluminum-induced protein with YGL and LRDR motifs (AYL1)).
Among 35 putative MPK3/6 substrates, nine are found in our data set and seven harbor the same phosphosites (Table  S8 ). Likewise, Lassowskat et al. reported a phosphoproteomic study on adult A. thaliana plants expressing the constitutively active PcMKK5 DD which activates MPK3 and MPK6  (33). The authors identified 538 probable phosphoproteins with altered abundance downstream of MPK3/6 activation, of which 30 are found in our dataset (Table  S8 ). Common putative MAPK substrates were thus identified in these studies, although plant age, growth conditions and the technical approaches differed. Besides, we identified AT2G20760 (clathrin light chain) which was shown to be phosphorylated by MPK3 and MPK4 on protein microarrays  (71) , and AT1G11360 (adenine nucleotide alpha hydrolases-like superfamily protein) and AT1G30450 (CCC1) which were predicted to be MPK3 or MPK6 substrates based on the experimentally determined preferential phosphomotif of these kinases  (72).
Similarly, some of the acidophilic motifs detected in the genotypedependent group could be the targets of casein kinases (CKs)  (78). CK2 protein kinases are involved in multiple biological processes and in particular in immunity  (79, 80). Our study might thus reveal new putative CDPK and CK targets that work in conjunction with MAPKs. Our global phosphoproteomic analysis of WT plants and mapk mutants suggests the existence of a high intricacy between the immune-regulated MAPKs and their substrates. A selection of putative MAPK substrates was analyzed to assess the ability of immune MAPKs to both interact and phosphorylate these candidates. BiFC and pull-down data revealed interaction of these putative substrates with none up to the three MAPKs. BBC1 did not interact with any of the three MAPKs in these assays, and several candidate proteins were shown to interact with the MAPKs only in BiFC or pulldown assays.
Acknowledgments-We would like to thank  Veronique Legros  and  Huoming Zhang  for technical assistance in MS, and  Christophe Bruley ,  Veronique Dupierris  and  Emmanuelle Lastrucci  (EDyP) for informatics support. We would like to thank  Carlos Xavier Pita  and  Heno Hwang  from  KAUST for scientific Illustrations.
* This work was supported by the  Agence Nationale de la Recherche  ( ANR-2010-JCJC-1608  to D. P. ) and by the  King Abdullah University of Science and Technology  ( BAS/1/1062-01-0 ). The IPS2 benefits from the support of the  LabEx Saclay Plant Sciences-SPS  ( ANR-10-LABX-0040-SPS ). □ S This article contains supplemental material. Center for Desert Agriculture, 4700 King Abdullah University of Science and Technology (KAUST) ,.