dataesr/scientific-paragraphs-categorization huggingface dataset:
S14 ). Except for the root, the points of the infinite-sites plot from the combined dataset formed a straight line, indicating that the relatively high uncertainty in time estimates was mostly due to uncertainties in fossil calibrations  (Rannala and Yang 2007). Setting the root maximum bound to 38 Ma had the strongest impact on the divergence time estimates and resulted in more ancient estimates, with larger differences observed near the root. When we used constraints on minimum bounds, posterior estimates tended to be more recent and credibility intervals larger. Despite the fact that using four calibrations resulted in a slight truncation of the prior density on the age of the Calonectris -Ardenna node, mean posterior time estimates were nearly identical and variances were higher when only using three calibrations (Fig.
S18 ). Nonetheless, the observed values at the extremes tended towards the expected values when increasing the BS threshold for collapsing branches in the gene trees. When branches with BS < 50 were collapsed, the TICR test was no longer significant, suggesting that the excess of outlier quartets was caused by including noise in form of inaccurate gene tree branches. Using a program like BUCKy  (Larget et al. 2010)  to calculate the CF from gene trees can provide an advantage, since it also considers uncertainty in gene tree estimation. It is noteworthy that in all PE-ddRAD datasets (but not in UCE datasets), P. boydi and P. lherminieri were found more frequently in 4-taxon sets that did not fit the tree model with ILS compared to any other taxon. The TICR test for Ardenna was not significant and all 4-taxon sets fitted the tree model with ILS.
Consistent with the TICR test results, we detected a lack of sharp improvement in the slope heuristic in both Puffinus and Ardenna datasets. Nonetheless, we detected a general sharper decrease up to h=2 in Puffinus, although inferred reticulation events had γ < 5%. The only introgression event with γ > 5% was from P. lherminieri to P. boydi. For candidate networks assuming introgression between P. lherminieri and P. boydi, log pseudo-likelihood values obtained from PE-ddRAD data were either lower than or similar to those obtained in optimised networks with h=1, showing support for introgression between these taxa (Supplementary Table  S6 ). The directionality of introgression could not be determined reliably because optimisation of candidate networks showed very similar log pseudolikelihood values for pairs of candidate networks with reversed direction of introgression.
DISCUSSION Our results demonstrate the power of integrating UCE and RAD markers for resolving the phylogenetic relationships of a group of pelagic seabirds characterised by rapid diversification events that have confounded previous phylogenetic studies. To our knowledge, our study is the first to compare and integrate UCE and paired-end ddRAD datasets in a phylogenomic context using comparably phased sequence alignments for both datasets. Here, we propose a strategy to optimise RAD-Seq data for phylogenetic analyses, we consider aspects of our methodological approach that may be of help to future studies, we discuss case-by-case how the integrative use of PE-ddRAD-Seq and UCE data with phylogenetic and introgression analyses allows identification of the causes of phylogenetic discordance, and we discuss the systematic implications of our phylogenetic results.
RAD-Seq Dataset Optimisation for Phylogenetic Analyses There are two main factors that affect the number of orthologous loci recovered in RAD-seq datasets: 1) divergence times between lineages and 2) filtering and assembly parameters applied to orthology inference. Considerable research attention has focussed on exploring how the number of orthologous loci decreases with increasing divergence times  (Rubin et al. 2012; Cariou et al. 2013 ) and how filters based on taxon coverage affect dataset size and the ability to resolve phylogenetic relationships  (Wagner et al. 2013; Leaché et al. 2015; Díaz-Arce et al. 2016; Tripp et al. 2017). However, phylogenomic studies have generally neglected the optimisation of assembly parameters in order to minimise the inclusion of paralogous or repetitive loci (but see  Hosegood et al.
2020) , a procedure that is common practice in population genomics analyses  (Paris et al. 2017; Rochette and Catchen 2017). To fill this gap, we used a two-step optimisation process to assess the impact of http://mc. manuscriptcentral. com/systbiol common issues in RAD-Seq, such as sequencing error and paralog content, on phylogenetic reconstruction. However, our analyses revealed a minor effect of assembly parameters on phylogenetic reconstruction, compared to taxon coverage. Our study adds to previous evidence showing a major effect of taxon coverage on phylogenetic reconstruction when using RAD-Seq datasets  (Díaz-Arce et al. 2016; Tripp et al. 2017 ). We recommend that phylogenetic studies using RAD-Seq data should explore several taxon coverage filters in order to maximise the phylogenetic informativeness of their datasets.
Considerations on Methodological Approaches for Phylogenetic Inference As expected, differences in methodological approaches can severely affect the phylogenetic inference. First, in cases of extreme levels of ILS, concatenation analyses may result in an average topology that differs from the true species tree  (Mendes and Hahn 2018). On the other hand, low phylogenetic information per locus in PE-ddRAD and UCE datasets might result in poorly resolved gene trees, and in such cases, concatenation methods can be more accurate than summary coalescent approaches  (Mirarab et al. 2016; Springer and Gatesy 2016). In our case, concatenation and summary coalescent approaches produced largely congruent topologies that expressed a lower degree of incongruence than analyses using different datasets and different levels of missing data.
This was due to the inability of accurate phasing towards the extremes of UCEs because of lower read coverages (we required a minimum of 5 reads per haplotype to include a position). This approach is more conservative than using the full contig sequence and therefore results in a reduction of the number of PIS. However, it delivers a higher reliability in base calling and a true representation of polymorphism that can be particularly useful at recent timescales. Researchers working with UCE data face a tradeoff between longer alignments with higher amounts of PIS and more reliable base calling allowing the inclusion of polymorphism data. Third, partitioning strategy can affect the outcome of phylogenetic analyses. For example, maximum likelihood, when used to analyse an unpartitioned concatenated alignment from different loci, can converge to a tree other than the species tree as the number of loci increases  (Roch and Steel 2015).
Our concatenation analyses using the partitioned UCE 75 dataset recovered exactly the same topology compared to unpartitioned analyses, with only slight changes in branch support, showing that the topology was relatively insensitive to the partitioning scheme used. This was likely due to a smaller effect of evolutionary rate heterogeneity among loci at shallow phylogenetic scales like the one we were working than at deep timescales. Fourth, increasing the number of individuals strongly improves species tree estimation in ASTRAL-III when branch lengths are extremely short  (Rabiee et al. 2019). We found evidence for this pattern in our UCE ASTRAL-III trees, where we found incongruences affecting nodes where individuals had been removed to avoid the inclusion of fragmentary sequences. We also observed an improvement in precision when analyses were performed after contracting very low support branches (BS < 10) in the gene trees, as previously suggested by  (Zhang et al.
Divergence dating with UCE and PE-ddRAD Collins and Hrbek (2018) showed consistent discrepancies in divergence time estimates using different phylogenomic markers under strict and relaxed clock models. The discrepancies were explained by temporal differences in phylogenetic informativeness (PI). However, their analyses showed that UCE and ddRAD datasets showed similar temporal patterns of PI and resulted in similar divergence time estimates. Our analyses using UCE and PE-ddRAD empirical datasets also showed very similar estimates. However, UCE estimates tended to be slightly older and have wider credibility intervals than PE-ddRAD estimates probably due to UCE datasets having lower PI at the timescales covered in this study. Theoretically, for infinitely long alignments, an infinite-sites plot, which measures the uncertainty in the divergence time posterior (width of the credibility interval vs.
posterior mean of node ages) should converge onto a straight line  (Rannala and Yang 2007). The slope of this line represents the amount of uncertainty in time estimates per 1 million years (Myr) of divergence solely due to uncertainties in the fossil calibrations. The points of the infinitesites plot from the UCE and the PE-ddRAD datasets did not form a straight line, suggesting that uncertainties in time estimates were due both to limited data as well as uncertainties in the fossil calibrations. On the other hand, the points from the TENT dataset formed a straight line (Supplementary Fig. S14 ), indicating that the relatively high uncertainty in time estimates was mostly due to uncertainties in fossil calibrations and that increasing the amount of data would only marginally reduce the credibility intervals  (Rannala and Yang 2007).
Integrative Approach using UCE and PE-ddRAD to Disentangle Phylogenetic Discordance The utility of RAD-seq and target capture approaches for phylogenetic estimation across different timescales has been widely demonstrated  (Faircloth et al. 2012; Cruaud et al. 2014; Smith et al. 2014; McCluskey and Postlethwait 2015). However, only a handful of studies have explored the utility of both approaches in phylogenetic studies  (Leaché et al. 2015; Harvey et al. 2016; Manthey et al. 2016; Collins and Hrbek 2018). Using a higher number of taxa and loci than these previous studies, we show the advantages of integrating PE-ddRAD-Seq and UCE data to infer the phylogenetic relationships of a challenging group, the shearwaters, across a range of timescales using concatenation and coalescent approaches.
Despite finding only minor data-type effects, datasets from different markers and levels of missing data tended to better resolve short internodes at different timescales (Fig. 4  and Supplementary Fig. S5 ). For instance, the UCE 95 dataset, which contained the lowest number of PIS, showed the highest support across methods (with the exception of SVDQuartets) for the sister relationship between Ardenna and Calonectris (18. 5 Ma), but showed the poorest performance at recovering the monophyly of recently diverged taxa (Fig. 4 ). This shows that the phylogenetic signal of this dataset is stronger near the root and weaker at shallow timescales. Conversely, the dataset with the highest number of PIS (PE-ddRAD 75) recovered alternative topologies (with low support) for the relationships between the three genera, but recovered the monophyly of all the species (and subspecies).
Finally, phylogenetic analyses using the TENT 75 matrix performed well at both deep and shallow timescales, showing the advantages of combining different data types in a single analysis. Our study is consistent with previous findings that combining data types leads to higher resolution on short internodes  (Jarvis et al. 2014) , and it is noteworthy that, due to high rate heterogeneity between PE-ddRAD loci and UCEs, the observed data-type effects could also reflect poor model fit  (Reddy et al. 2017). Despite the strong support in most of our phylogeny, we found phylogenomic conflict associated with short internodes in three areas of the tree  (Degnan and Rosenberg 2006). Short internode lengths are usually associated with topological conflict. This conflict arises because short speciation intervals 1) accumulate few substitutions, resulting in a low number of informative sites and, 2) they increase the probability of finding different gene histories due to ILS  (Alda et al.
2019). The first area of conflict was the short internode near the root separating Ardenna and Calonectris from Puffinus (triangle in Fig. 3 ). When using the PE-ddRAD dataset with the highest number of PIS, a sister relationship between Ardenna and Puffinus received the highest quartet support in ASTRAL-III analyses. The remaining datasets provided the highest quartet support for a sister relationship between Ardenna and Calonectris and the support for this arrangement decreased when increasing the number of PIS in the dataset (Table  1  and 2 ). These results show how high levels of ILS resulted in different phylogenetic signals across datasets and could indicate that our most variable dataset is affected by homoplasy due to saturation at this timescale. In addition, the increased support for the Ardenna + Puffinus clade compared to the alternative topology in PE-ddRAD datasets could also indicate molecular convergence between Ardenna and Puffinus, which show similarities in morphology and diving behaviour and were historically placed in the same genus  (Penhallurick and Wink 2004).
The split of the three North Atlantic lineages of Puffinus (1: P. lherminieri, P. baroli and P. boydi, 2: P. puffinus, and 3: P. mauretanicus and P. yelkouan) (Fig. 3a ) represents the second conflict associated with a short internode. In this case, our analyses largely supported a sister relationship between P. puffinus and the Mediterranean clade. Nonetheless, only the PE-ddRAD 75 and the TENT 75 datasets clearly rejected the polytomy test (Table  2 ). When two speciation events occur at the same time (hard polytomy) and gene tree heterogeneity is mostly caused by ILS, the multispecies coalescent model (MSC) predicts equal frequencies for each of the three topologically-informative unrooted quartet topologies defined around the short internode  (Degnan and Rosenberg 2009). Our data show that these frequencies are nearly equal for most analyses (Table  2 ).
In addition, our introgression analyses confirmed that no introgression had occurred between P. puffinus and the P. boydi, P. baroli and P. lherminieri lineage. These observations suggest a nearly simultaneous divergence of the three lineages that resulted in high levels of ILS. The last short internode-associated phylogenetic conflict was the relationship between the two all-dark coloured Ardenna species (A. tenuirostris and A. grisea). Concatenation analyses showed strong support for the sister relationship between A. grisea and the A. gravis, A. creatopus and A. carneipes clade across data type and levels of missing data (Fig. 4 ). On the other hand, the alternative topology of a sister relationship between the two alldark species received ASTRAL-III quartet supports that were nearly as high as (PE-ddRAD data) or higher (UCE data) than the main topology, and at least 8% higher than the remaining alternative topology.
For instance, SNaQ assumes that no two edges can be part of the same cycle (Solís-Lemus and Ané 2016). This assumption precludes the possibility of inferring gene flow between two taxa in both directions in the same analysis. In addition, gene tree distributions from these alternative topologies may be indistinguishable  (Pardi and Scornavacca 2015). However, the network with P. boydi as the recipient of genetic material from P. lherminieri had much higher inheritance probabilities than those with P. lherminieri as recipient. These facts, together with increased heterozygosity in P. boydi compared to P. lherminieri and P. baroli, point to P. boydi as the recipient of genetic material. Optimised inheritance probabilities on the fixed network with gene flow from P. lherminieri to P. boydi were lower when using UCE data than when using PE-ddRAD data.
ASTRAL-III quartet supports for the main and the two alternative quartet topologies and polytomy test p-values for the challenging branches highlighted in Figure2. The null hypothesis is polytomy and p-values < 0. 05 reject a real polytomy. Polytomy test p-values are shown in parentheses when the alternative topology was recovered. Note that only the TENT dataset is able to reject real polytomies at the three challenging branches. Downloaded from https://academic. oup. com/sysbio/advance-article/doi/10. 1093/sysbio/syaa101/6123758 by guest on 09 February 2021.
ACKNOWLEDGMENTS We would like to thank the many institutions that provided tissue loans for this research:  Smithsonian National Museum of Natural History, University of Washington Burke Museum, American Museum of Natural History, Louisiana State University Museum of Natural Science, the University of Kansas Biodiversity Institute , and the  Muséum National d'Histoire Naturelle. We are grateful to  Gary Nunn ,  Jeremy J. Austin ,  Chris Gaskin ,  Kazuto Kawakami ,  Juan E. Martínez-Gómez , and  Maite Louzao  for collecting fresh or dry material from the field and/or providing tissue samples from museum skins for this study. We thank http://mc. manuscriptcentral. com/systbiol.
This paper presents and discusses the "PesticideLoadIndicator" package, a new R-package to compute potential environmental and health effects of pesticide applications using the Danish 'Pesticide Load' indicator. The implementation in the R Statistical Language makes it easy for researchers, practitioners and institutions to compare potential pesticide risks for a wide range of applications and compute risk indicators at field-, crop-, farm-, regional-or national level. The tool is publicly available. It provides a possibility for the direct integration of pesticide properties needed for indicator computation (ecotoxicity, environmental fate and human health), from the Pesticide Properties Database or other sources and allows users to change default reference values and weights.
To make the 'Pesticide Load' indicator easily accessible to a large range of stakeholders, e. g. for research, education and policy evaluation, we here present and discuss a package in the R Statistical Language (R Core Team, 2020) for indicator computation  (Möhring et al. , 2021). The package allows users to compute the Pesticide Load and its subindicators for potential Ecotoxicity, Environmental Fate and Human Health effects based on pesticide application data. It further provides a possibility for the direct integration of pesticide properties, needed for the computation, from the Pesticide Properties Database  (Lewis et al. , 2016). Finally, it is flexible to integrate other data sources and/or adjust default (Danish) reference values and weights for indicator computation. Following, we briefly describe the indicator, its background and limitations, provide a description of the package and its functions and give an example of its application using integrated test data.
Description of the pesticide load indicator The Pesticide Load and its sub indicators are computed for each pesticide product (eq. 1 & 2). They may then be aggregated at a field and crop-, farm, regional or national level. We further provide the Load Index in the package, given applied standard dosages are provided. (1) Pesticide Load [Load units per kg/l product] = Ecotoxicity Load + Environmental Fate Load + Human Health Load (2) Load Index [Load per standard dosage] = Pesticide Load * Applied standard dosages of product (standard dosage of product defined in national pesticide databases) The Ecotoxicity Load is calculated, based on 11 LC/LD/EC50 (median lethal dose/concentration) and NOEC (No Observed Effect Concentration) parameters. Eight of these parameters relate to acute ecotoxicity/ short-term effects of active substances on birds, mammals, fish, daphnia, algae, aquatic plants, earthworms, bees, while three relate to chronic ecotoxicity/long-term effects on fish, daphnia and earthworms, respectively.
Values for chronic ecotoxicity are adjusted for degradation rate of the active ingredient in case half-life values in soil and water are available. The Environmental Fate Load is based on three parameters: the halflife in soil (DT 50 ), the bioaccumulation factor (or the log P ow value if no BCF value is reported) and the SCI-GROW index that reflects the mobility and risk of leaching to the groundwater of the active ingredient and the major metabolites listed in the Pesticide Properties Database. Whereas the calculation of the Ecotoxicity and Environmental Fate Loads are based on properties of the active ingredients in the product and are then aggregated on a product level, the Human Health Load is directly calculated based on the risk phrases and symbols on the product label and the formulation of the products.
(2018)  and  Möhring et al. (2019). It should be stressed that the "PesticideLoadIndicator" tool made available here is not identical to the Danish Pesticide Load used by the Danish authorities and the calculated pesticide loads may therefore be different. This is because pesticide products and their formulation may differ and change over time and country and pesticide properties data from the Pesticide Properties Database may have been amended and completed. The here computed indicators should further not be confused with the so called "Pesticide Load Indicator (PBI)", which indicates Load on a national level divided by a reference area (conventionally grown and sprayed crops in 2010) and is sometimes used in Denmark.
Data availability, requirements, and preparation Computation of the Pesticide Load and the Environmental Fate, Ecotoxicity and Human Health Load with the package requires pesticide application data and human health properties (risk statements) on a product level, as well as information on active ingredients of the applied pesticides, their concentration and fate, as well as ecotoxicity properties, as indicated above. Data on standard dosages and human health properties (risk statements) of products is publicly accessible in national pesticide registers. Substance properties data used for the computation of the 'Pesticide Load' indicator is publicly accessible, e. g. through EU dossiers on the authorization of pesticide products. A central collection of this information is, for example, provided by the Pesticide Properties database (PPDB,  Lewis et al.
, 2016). The here presented tool allows users to integrate a database version of the PPDB (licensed) and automatically draw required information for indicator computation 1 or provide information from other information sources themselves. Due to the heterogeneity (in quality) of pesticide use data and the low accessibility of national databases on registered pesticides, the computation of the 'Pesticide Load' indicator often requires substantial work to collect and prepare data. Further tools for data scraping and manipulation, complementing the here presented tool for indicator computation should be developed to render pesticide risk indicators more accessible.
Reproducible example Following, we provide a reproducible R code on how to compute the pesticide indicators and discuss the outcome. We compute the indicators for four pesticide applications in Switzerland (example data provided in the package): # Install and load the PesticideLoadIndicator package install. packages("PesticideLoadIndicator") library(PesticideLoadIndicator) # Load the dataframe containing the pesticide use data. 1 Further, note that the package is regularly updated to account for potential changes in the PPDB so that workload for users with an updated version (i. e. licensed users) is significantly reduced. CRAN further provides the possibility to use previous version of the package, which might be more compatible if older versions of the PPDB are used. Please see the package description for a reproducible example of how to compute the indicators only based on simple pesticide use data, i.