dataesr/scientific-paragraphs-categorization huggingface dataset:
5. This first application of the eco-evolutionary multi-species model Ev-OSMOSE to the North Sea opens the field of eco-evolutionary studies to marine ecosystems models. This study underlines the parameterization feasibility in spite of the high data quality requirement to parameterize the phenotypic and genotypic variances of life-history traits. Ev-OSMOSE-NS is the first configuration to account for genotypic and phenotypic variances of several interacting species and succeeds to improve the simulated variances of life-history traits. It is an important step toward more realism notably in representing length-at-age distribution and the maturation process. Ev-OSMOSE-NS is also, to our knowledge, the first multi-species model applied to a marine ecosystem that accounts for mendelian inheritance of traits from parents to their offspring for all the species of a food web simultaneously, thus allowing to account for the micro-evolution of exploited species in response to selective pressures such as fishing and climate change together with their coevolution due to trophic interactions.
Acknowledgements This work has been partially funded by the  BiodivErsA and Belmont Forum  project  SOMBEE   (BiodivScen programme ,  ANR  contract n• ANR-18-EBI4-0003-01 ). Alaia Morell  was supported by a PhD grant from  Ifremer  and  Région Hauts-de-France. The  POLCOMS-ERSEM  projections were produced with funding from the  European Union   Horizon 2020 research and innovation programme  under grant agreement No  678193  ( CERES ,  Climate Change and European Aquatic   Resources). The authors acknowledge the  Pôle de Calcul et de Données Marines  (PCDM, http://www. ifremer. fr/pcdm) for providing DATARMOR storage, data access, computational resources, visualization and support services. Yunne-Jai Shin  acknowledges funding support from the  European Union 's  Horizon 2020 research and innovation program  under grant agreement No  869300  ( FutureMARES ) and the  Pew marine fellows program.
Metagenomic Sequence Processing Quality control of sequences was performed by first removing adapters and filtering raw sequences with bbtools v37. 78 (bbduk with parameters fastawrap ¼ 300 k ¼ 23) and the software Skewer v0. 2. 2  32  with the parameters of -n -l 100 -q 25. The bbmap command, part of the bbtools package (https://sourceforge. net/projects/bbmap/), was then used to remove any human reads potentially contaminating the samples (parameters minratio ¼ 0. 9 maxindel ¼ 1 bwr ¼ 0. 16 bw ¼ 12 fast minhits ¼ 2 qtrim ¼ r trimq ¼ 10 untrim idtag printunmappedcount kfilter ¼ 25 maxsites ¼ 1). The fastqc software (v0. 11. 7) and MultiQC sofware (v1. 5. dev0 88e1e4f) was used to generate summaries of final QC sequence counts. Reads passing quality control were then taxonomically profiled and their relative abundances estimated with Kracken2  33  and Bracken (v2.
Estimation of transcript abundance was calculated, and the count values were normalized to the upper quartile of the fragments per kilobase of transcript per million fragments mapped reads using Cufflinks (cuffdiff). The processed data were then analyzed on software R for preparations of linear gene expression comparison, PCA, and heatmap figures. Heatmap figure represent the gene presenting a 2-fold increase and a P value <. 05 between CDD and Psy diet. Webgestalt was used for the gene set enrichment analysis. A gene list for Webgestalt analysis was generated using the output from the cuffdiff program using the same set of gene than for the heatmap.
Statistical Analyses Statistical significances of results were analyzed using GrapahPad 9 software. Normality test was performed using the Shapiro-Wilk test. Experiments in which results showed normal distribution were tested by unpaired Student t test or analysis of variance. Experiments yielding data not following a normal distribution were analyzed by the Mann-Whitney test or Kruskal-Wallis test. When specified experiments ware performed multiple times, each showing a similar pattern of results, which mean that P values were similar between experiments (± 1 log 10 ). Alpha and beta diversity were calculated using QIIME2. 36 Bracken derived abundances were imported into QIIME2 and rarefied to 82502 sequences per sample. Shannon alpha diversity (with Kruskal-Wallis tests between diet groups) and Bray-Curtis dissimilarity were calculated using the diversity coremetrics plugin.
The Shannon alpha diversity values were plotted using the R package ggplot2. Differences between the microbial taxonomic abundances of experimental groups were tested with LEfSe. 37 The bracken OTU table was collapsed to the family level of taxonomy and rarefied to 82502 sequences per sample using QIIME2. The resulting table was then imported into the galaxy tool for LEfSe (https://huttenhower. sph. harvard. edu/galaxy/) for the 2 diets, CCD and psyllium, being compared. The default parameters were used for LEfSe, except for the threshold on the logarithmic LDA score, which was increased from 2. 0 to 4. 0. Significance is expressed as *P <. 05, **P <. 01, ***P <. 001, and ****P <. 0001.
Similarly, coalescent species delimitation was consistent among similarity threshold values. However, we identified a strong effect of the bCT on the branch lengths of concatenation and species trees, with higher bCTs yielding trees with shorter branches, which might be a pitfall for downstream inferences of evolutionary rates. Our results suggest that the choice of assembly parameters for RADseq data in the context of shallow phylogenomics might be less challenging than previously thought. Finally, 3 we propose a pipeline for empirical optimization of the iCT and bCT, implemented in optiRADCT, a series of scripts readily usable for future RADseq studies.
Finally, enrichment PCR was performed to amplify the library using forward and reverse RAD primers. Sequencing was conducted on an Illumina Next-Seq at Glasgow Polyomics to generate 75 bp paired-end reads. After quality check using FastQC (www. bioinformatics. babraham. ac. uk/projects/fastqc/), raw reads were demultiplexed with stacks v2  (Rochette et al. 2019 ) and trimmed to 60 bp. Raw reads were submitted to the NCBI Short Read Archive under Bioproject ###### (to be added upon manuscript acceptance)  (open boxes = "split" scenarios, solid boxes = "lump" scenarios, cf. methods for more details).
Hybrid-Enrichment analyses As a comparative basis for evaluating the quality of phylogenetic analyses performed with RADseq data, nuclear loci were also sampled through a Hybrid-Enrichment protocol, using a modified version of the FrogCap sequence capture probe set  (Hutter et al. 2022). In brief, our goal was to create a set of universal amphibian markers that could be targeted using a single set of probes. We started by using the FrogCap markers  (Hutter et al. 2022;  https://github. com/chutter/FrogCap-Sequence-Capture) that were successfully captured broadly across Anura in the data from  Hutter et al. (2022). Next, these markers were matched against the Salamander Ambystoma mexicanum genome  (Keinath et al. 2018 ) and retained if they matched with ≥65% similarity. The total number of markers was 7,720 (Ultra Conserved Elements [UCEs]: 2,122; exons: 5,598).
We further aimed to add additional universal markers that could potentially be compared and included in phylogenetic analyses with other groups of organisms. We selected the USCO (Universal Single Copy Orthologs) set of markers, which are orthologous genes found across deep evolutionary scales (i. e. , Kingdom, Phylum, Class). To identify USCO markers present across Amphibia, we selected genomes of Ambystoma mexicanum and the frogs Xenopus tropicalis, Rana catesbeiana, Nanorana parkeri, Rhinella marina, and Oophaga pumilio  (Sun et al. 2015; Hammond et al. 2017; Edwards et al. 2018; Keinath et al. 2018; Rogers et al. 2018). We used the program BUSCO v3. 0. 2  (Seppey et al. 2019)  to locate and identify the Metazoa and Tetrapoda USCO markers in these genomes. We included all metazoan genes found in all the genomes, which totalled 677 out of 978 genes.
The genomic libraries for the samples were prepared by Arbor BioSciences library preparation service. Library preparation is detailed in the Supplementary Material S1. Illumina sequence data were de-multiplexed using the Illumina software bcl2fastq. The bioinformatics pipeline for filtering adapter contamination, removing biological contamination, assembling contigs, and preparing sequence data for alignment is available at https://github. com/chutter/FrogCap-Sequence-Capture, and described in more detail in Supplementary Material S1. Next, the retained set of markers was aligned on a marker-by-marker basis using MAFFT local pair alignment (settings: max iterations=1000; ep=0. 123; op=3; --adjustdirection;  Katoh & Standley 2013). We screened each alignment for samples ≥40% divergent from consensus sequences, which were almost always incorrectly assigned to contigs.
RADseq data assembly and parameter optimization For each genus, RADseq de novo assemblies were generated from the demultiplexed reads using ipyrad v0. 9. 50  (Eaton & Overcast 2020) , a software that processes RADseq data through seven sequential analytical steps (referred to later in this paragraph, see https://ipyrad. readthedocs. io/en/latest/7-outline. html for a detailed overview). In this pipeline, the CT is determined by a single parameter, the clustering threshold (#14 in the parameter file) which determines the percentage of similarity needed between two sequences to be clustered (i. e. , considering sequences of 100 bp, a clustering threshold of 0. 90 would cluster sequences with ≤10 substitutions). This parameter will be used at step 3 for intra-sample clustering of the reads (intrasample Clustering Threshold, hereafter referred to as iCT) and at step 6 for between-samples clustering of the consensus sequences (between-samples Clustering Threshold, hereafter bCT).
g. ,  Eaton & Ree 2013; Paetzold et al. 2019). All the other parameters were left as default. The s7 summary files were used to collect loci-related metrics: number of loci in the final assembly, proportion of loci shared by ≥50% and ≥80% of the samples, total number of variable and parsimony-informative sites, and percentage of missing data in the concatenation and SNPs matrices. Custom scripts in R (R core team 2020) were used to visualize the impact of the iCT on these variables (all scripts are available at github. com/rancilhac/optiRADCT). Particularly, for each variable, a rate of change Δ was calculated as ∆ሺܺሻ ்ୀே ൌ ܺ ்ୀே െ ܺ ்ୀேିଵ. Finally, the concatenated loci from each assembly were used to calculate ML phylogenetic trees using RAxML v8 under a GTR+Γ substitution model. The rapid bootstrapping algorithm was run with 100 replicates.
The resulting trees were compared to the respective HE reference ML tree using two metrics: the Robinson-Foulds distance (RFD;  Robinson & Foulds 1981 ) that measures topological difference only, and Kuhner & Felsenstein's Branch Lengths Score (BLS;  Kuhner & Felsenstein 1994) , an extension of the previous that accounts for differences in both topology and branch lengths. It should be noted that four samples were missing from the HE datasets (three Discoglossus and 1 Triturus), and one from the RADseq datasets (Appendix 1). These samples were removed to compare the trees. Distances were computed using custom R scripts based on ape v5. 6-2  (Paradis et al. 2004) , which was also used to assess changes in phylogenetic uncertainty (i. e. , bootstrap support of the branches [BS]) and branch lengths, depending on the iCT.
Impact of the bCT on phylogenetic and Multi-Species Coalescent analyses For the 92 assemblies, the concatenated loci were used to calculate ML phylogenetic trees with RAxML, using the same settings as before. The quartet-based method tetrad (https://github. com/ eaton-lab/tetrad) was also used to infer species trees from the unlinked SNP matrices outputted by ipyrad (i. e. , the best-covered SNP is selected from each locus or taken as random in case of equality). The impact of the bCT on the topologies, resolution and branch lengths of the inferred trees was investigated with the same approach as the iCT. As tetrad does not compute branch lengths, only topologies and supports were compared in this case. For the Discoglossus and Rana data sets, we conducted further MSC-based inference of species tree and species limits using the SNAPP/BFD* framework  (Bryant et al.
2012; Leaché et al. 2014)  implemented in BEAST 2  (Bouckaert et al. 2019). Because this approach is very sensitive to missing data  (Schmidt-Lebuhn et al. 2017) , unlinked SNP matrices were filtered to keep only sites covering all samples. Prior to that, individuals with missing data at >85% of all SNPs were removed to mitigate the data loss. In the Discoglossus data set, this resulted in the exclusion of D. montalentii as an effect of loci dropout but including it would have drastically reduced the SNP matrices, making it difficult to determine whether uncertainty was due to an improper bCT or the low amount of data. First, SNAPP was run seven times for each genus, using SNP matrices assembled with the "optimal" bCT as identified in the previous steps (0. 91 in both cases, cf. Results), and six alternative bCTs representing various levels of over-merging (0.
RADseq data: effect of the intra-sample Clustering Threshold For this first step, 23 assemblies were generated for each data set. Loci sequences and SNP matrices outputted by ipyrad, as well as spreadsheets providing detailed metrics for each assembly, have been uploaded to https://zenodo. org/record/7829243 (DOI: 10. 5281/zenodo. 7829243). In all cases, the number of constructed clusters increased, and their depth decreased, with increasing iCT (Fig. S5 ). Both variables showed a break at iCT=0. 95 in all four data sets, above which the variation generated by each increment of the iCT was substantially higher. Finally, the number of clusters rejected by the heterozygosity filter followed a plateau for iCTs between 0. 80 and 0. 91 before decreasing drastically and stabilizing near zero for iCT > 0.
1%) and concatenation (34. 5% -71. 8%) matrices was relatively stable for iCTs<0. 90, above which it increased drastically (Fig. 2 ). All data sets were largely concordant in supporting two distinct breaks in the rate of change of these variables, for values of iCT=0. 90 and 0. 95, although the trend was somewhat weaker for Triturus (Fig. 2 ). The first of these breaks was not very strong, but iCT=0. 95 marked a threshold above which each increment of the iCT generated a strong drop in the number of shared loci and variable and parsimony informative sites, as well as a strong increase in the proportion of missing data. Despite some dissimilarities, the concatenation-based phylogenetic analyses yielded overall consistent topologies across iCT values. All inferred trees in Newick parenthetical format are available at https://zenodo.
RADseq data: impact of the between-sample Clustering Threshold In the four genera, clusters-and loci-related assembly metrics were relatively stable for 0. 80<iCT< 0. 90 (Fig. 2 , Fig. S5 ). Especially, the change of three clusters-related metrics -the number of clusters, depth of the clusters and number of clusters of high heterozygosity -showed a break at iCT=0. 95 in all four genera. Loci-related metrics, on the other hand, were less straightforward to interpret, particularly because Triturus showed less variation than the three other taxa. While some signs of under-merging were noticeable from iCT=0. 90 (e. g. , increase of missing data in Discoglossus and Rana), we find that the most striking decrease of shared loci and variable sites, and increase of missing data, occur for iCT>0. 95 in all four taxa.
We interpret this pattern as indicating that iCT=0. 95 represent a threshold above which under-merging of heterozygous loci becomes pervasive (cf. Discussion for more details). Hence, we considered this value as the "optimal" iCT for the next analytical steps. The clusters constructed with iCT=0. 95 were assembled into loci with 23 bCTs ranging from 0. 50 to 0. 99. The effect of the bCT on the informativeness and completeness of the assemblies was overall similar to that of the iCT (Fig. 4 ; assemblies and metrics are available at https://zenodo. org/record/7829243). Trends in loci-related metrics identified two bCT values inducing noticeable breaks. In the one hand, the proportion of loci shared by ≥50% of the samples and proportions of missing data in concatenation and SNP matrices remained relatively constant for 0.
9% and an increase of 6. 8% in these two metrics, respectively. As for the iCT, concatenation-based phylogenetic analyses were resilient to changes in the bCT, and all assemblies yielded the same topologies as at the previous step with high support for interspecific relationships (all inferred trees are available in Newick parenthetical format at https://zenodo. org/record/7829243). In Discoglossus, the BS of two inter-specific branches decreased to minima of 76 and 87, respectively, for bCT≥0. 95, whereas bCTs≤0. 60 failed to support the monophyly of D. g. galganoi (BS≤59). In Rana, the BS of the branch separating R. parvipalmata and its closest relative decreased to a minimum of 87 for bCTs≥0. 98. Otherwise, all bCT values recovered species monophyly as well as inter-specific relationships with maximum support in all four genera.
Discussion In this study, we explore the effects of varying intra-and between-samples clustering thresholds for de novo assembly of multispecies RADseq data sets on downstream analyses widely used in systematics studies. Although our results are applicable to studies using RADseq markers at any scale, they perhaps bear relevance for studies with limitations on number of individuals sequenced per species, such as is typical for phylogenomic, species delimitation, systematic and taxonomic endeavours.
These changes reflect the "under-merging phenomenon"  (McCartney-Melstad et al. 2019; Hühn et al. 2022) , where an increase of the iCT will ultimately lead to splitting alleles from heterozygous loci into independent clusters of lesser depth, while lower thresholds can lead to the clustering of non-orthologous sequences, creating fewer clusters of higher depth  (Paris et al. 2017 ). This will in turn affect the properties of the assembled loci because a higher iCT will increase the number of clusters discarded due to low read coverage, and artificially increase the number of invariant loci. In this context, proper optimization of the iCT appears crucial for de novo assembly of multispecies RADseq data sets. Hence, it is of interest to note that a noticeable break in the rate of change of most clusters-and loci-related metrics could be identified for iCT>0.
95, although trajectories varied depending on the metrics and taxa considered. For example, the number of clusters above minimum depth did not show a clear break when increasing the iCT, likely because the number of clusters being lost to the depth filter for high iCTs was compensated by the overall increase of the number of clusters. Regarding the loci-related metrics, the increase of missing data in Discoglossus and Rana at lower iCTs compared to the two newt genera could be due to their deeper crown age and the higher number of samples in these data sets, which could accentuate the effect of allelic dropout. The attenuated variation in most loci-related metrics in the Triturus data set remains an unexplained idiosyncrasy; given that no such pattern was observed in the second newt genus in our analysis (Lissotriton), the unusually large genome size of these urodeles cannot be invoked as an explanation for this pattern.
In the absence of clear evidence for "optimal" CT values, we conservatively propose to use the highest CT below the "under-merging threshold". Doing so, we could identify iCT and bCT values that yielded good quality data sets (i. e. , supporting well resolved trees concordant with reference topologies, in both concatenation and species tree analyses, and mitigating branch contraction), following a twostep approach. We thus propose the following decision framework to help optimize these parameters in future RADseq-based studies that include samples of different species-level lineages (summarized in Tab. 1): 1) Optimization of the iCT by identifying an "under-merging threshold". As a first step, the data set is assembled with a range of iCT (and a standard bCT). We employed a wide range of iCT values, but it can be narrowed as intra-specific diversity is expected to remain relatively low (e.
60) or high (≥0. 96) bCTs, although the topology itself remained constant. These relationships were also consistent with those recovered from the reference HE data sets, apart from the position of R. parvipalmata relative to R. temporaria (which could be due to influences of incomplete lineage sorting, as highlighted by concatenation-vs. species-tree differences in t he RADseq data). On the other hand, intra-specific relationships were more variable, but no clear relationship with either the iCT or the bCT could be established. Thus, this variation likely reflects genuine phylogenetic uncertainty at the "grey zone " between intra-and inter-specific levels, where gene tree heterogeneity is expected to be higher due to gene flow and ancestral allele sharing among populations. In contrast to the topology, we found that the branch lengths of the trees were strongly affected by the bCT, and to a lower extent by the iCT, with both parameters having distinct effects.
On the one hand, using greater iCTs induced a gradual, near-linear increase of the branch lengths, without any clear break useful for iCT optimization. On the other hand, increasing the bCT resulted in shorter branches, and the amplitude of the changes was much larger than that observed with the iCT. Branch length estimates are expected to vary according to several parameters, including the size and informativeness of the data set  (Schwartz & Mueller 2010)  and the proportion of missing data  (Leaché et al. 2015b ). Especially, in RADseq studies, using a higher Clustering Threshold will likely bias the assemblies toward slowevolving loci, resulting in under-estimated branch lengths. Concordantly, we observe that branch length follows a near-linear relationship with the number of PIS and proportion of missing data for bCT>0.
These results are consistent with the observed reduction of the concatenation trees' branch lengths and might be an effect of the smaller size of the SNP matrices used, directly resulting from the loss of variability due to under-merging of loci. Furthermore, Bayes Factor species delimitation analyses consistently supported the same model ranking, regardless of the SNP matrix used. In both genera, "split" scenarios were slightly favoured compared to the current taxonomy, while "lump" scenarios always produced BF<1, a result consistent with recent studies showing that MSC species delimitation detects population structure rather than species limits (e. g. ,  Sukumaran & Knowles 2017). However, the stable results recovered for both phylogenetic and MSC species delimitation analyses suggest that candidate species inferred from RADseq data based on either approach would be negligibly affected by suboptimal clustering threshold specifications.
Conclusions Our results are encouraging for future RADseq-based phylogenomic studies. Indeed, other genome-wide data collection approaches, such as sequence capture, have often been preferred to RADseq for multiple reasons, including de novo assembly pitfalls  (Leaché et al. 2015a; McKain et al. 2018; Hühn et al. 2022). We show here that RADseq and sequence capture data can yield equally accurate results, and that the clustering threshold is not a primary factor influencing the phylogenetic relationships inferred from the former. While unrealistic values, particularly those inducing strong under-merging, had an impact on the reconstructed relationships, moderate deviations from the "optimal" values produced well-supported topologies consistent with independent lines of evidence. Furthermore, in congruence with a previous study on plants  (Clugston et al.