dataesr/scientific-paragraphs-categorization huggingface dataset:
2012; Smith et al. 2014). These two methods have recovered concordant phylogenetic relationships when using large enough datasets  (Leaché et al. 2015; Harvey et al. 2016; Manthey et al. 2016; Collins and Hrbek 2018). RAD-Seq and UCEs have also been compared in terms of divergence time estimation http://mc. manuscriptcentral. com/systbiol using fossil-calibrated molecular-clock models and have been shown to accurately estimate divergence times at recent timescales  (Collins and Hrbek 2018). Integrating these two approaches may provide a powerful tool for disentangling the roles of ILS and introgression in rapid diversification events; however, empirical studies to test this assumption are generally lacking. Modern birds provide many case studies, as their diversification has been characterised by a succession of rapid radiations, posing a significant challenge in resolving the phylogenetic relationships of several avian clades  (Jarvis et al.
Hybridisation has been documented in several species of Procellariiformes, including shearwaters  (Genovart et al. 2012; Booth Jones et al. 2017; Masello et al. 2019 ) despite strong philopatry to breeding colonies and, in most cases, a lack of overlap between breeding areas of closely related species. Thus, shearwaters demonstrate the challenges typical of many other taxonomic groups, where rapid diversification and introgression hinder attempts to confidently resolve their evolutionary history, and provide a good case study for how combining genomic datasets can aid in clarifying relationships at recalcitrant nodes. Here, we generate the first phylogenomic datasets for shearwaters. We use paired-end ddRADSeq (PE-ddRAD) and UCE datasets to explore the role of ILS and introgression as causes of phylogenetic discordance during rapid diversification.
We adopt a thorough and integrative approach, comparing and combining the two datasets, and applying state-of-theart concatenated and coalescent-based approaches using shearwaters as a case study. We expand previous comparisons of UCE and RAD-Seq datasets for phylogenetics, divergence time estimation and inference of introgression, and we propose a strategy to optimise RAD-Seq data for phylogenetic analyses. Our approach allows us to completely resolve the phylogenetic relationships of shearwaters. We detect high levels of gene tree discordance associated with short internodes, mainly caused by ILS, and report a potential case of historical introgression. We show that our integrative approach provides a good framework for exploring the processes underlying phylogenetic incongruence during rapid diversification events.
Sampling and Sequence Data Generation We obtained blood (n = 45), high-quality tissue (n = 19) or dry tissue (n = 3) samples for 30 taxa representing 26 of the 30 recognised species of shearwaters  (Carboneras and  Bonan 2019) (Supplementary Table  S1 ). We also sampled three species of Procellariiformes as outgroups: Fulmarus glacialis, a species from the same family as the shearwaters (Procellariidae); and two more distantly related outgroups, Thalassarche chlororhynchos (Diomedeidae) and Oceanites oceanicus (Oceanitidae). Species that could not be included (Puffinus heinrothi, P. bannermani, P. persicus and P. subalaris) are mostly very localised or critically endangered. We used the Qiagen DNeasy Blood and Tissue Kit to extract genomic DNA according to the manufacturer's instructions  (Qiagen GmbH, Hilden, Germany).
For dry tissue samples, we extracted DNA using the  Dabney et al. (2013)  ancient DNA extraction protocol. We used a Qubit Fluorometer (Life Technologies) to quantify and standardise DNA concentrations of all samples. PE-ddRAD data for the outgroups were retrieved from whole genome assemblies available from the Bird 10K project  (Feng et al. 2020 ; see description of the in silico digestion protocol in the Data Assembly -PE-ddRAD-Seq dataset section below). Library preparation, capture enrichment and sequencing of the UCEs was conducted by RAPiD Genomics, LLC (Gainesville, FL, USA). Briefly, genomic DNA was fragmented, and sequencing libraries were prepared. Indexed samples were subjected to PCR for 7 cycles http://mc. manuscriptcentral. com/systbiol prior to pooling. UCE loci were captured using the Tetrapods UCE-5Kv1 probe set (available at ultraconserved.
py script. The remaining UCE sequences were aligned using MAFFT v. 7. 130b  (Katoh and Standley 2013)  and internally trimmed using Gblocks v. 0. 91b  (Castresana 2000). To retrieve polymorphism information lost when collapsing multiple reads into a single contig sequence, we used the phasing protocol described by  Andermann et al. (2018)  to obtain International Union of Pure and Applied Chemistry (IUPAC) consensus sequence alignments comparable to the PE-ddRAD alignments. Finally, we assembled datasets containing UCE loci that were present in at least 75% and 95% of the taxa using contig alignments (UCE 75 contig and UCE 95 contig) and using IUPAC consensus sequence alignments (UCE 75 IUPAC and UCE 95 IUPAC) (Fig. 1 ). PE-ddRAD-Seq dataset. -We quality-filtered and demultiplexed reads using process_radtags in Stacks v2.
41  (Rochette et al. 2019). To obtain PE-ddRAD datasets optimised for phylogenomic analyses, we performed a two-step approach to assembling RAD loci (Fig. 1 ). Here we provide an overview of this approach, which is described in detail in the Supplementary Information. First, we used a method for optimising de novo assembly of loci using Stacks  (Paris et al. 2017; Rochette and Catchen 2017)  which consists of varying each of the two key parameters (M: within-individual distance parameter, and n: betweenindividual distance parameter;  Catchen et al. 2011 ) separately using denovo_map. pl, and selecting values of M and n under which the tendency of new polymorphic loci when increasing M or n becomes linear. We also assessed the frequency of putative sequencing errors in the data by selecting the value of m (stack-depth parameter) where the proportion of singletons in the dataset stabilised  (Harvey et al.
2016). Secondly, we used two metrics to assess the effect of different pipelines and different parameterisations when assembling PE-ddRAD datasets: the number of parsimony informative sites (PIS) per locus relative to total locus length (pPIS); and the sum of bootstrap branch support (BS) obtained from maximum likelihood analyses (ML). For the second step, we used seven parameterisations from two assembly pipelines, Stacks and pyRAD 3. 0. 66 (Eaton 2014 ) (Fig. 1  and Supplementary Table  S2 ). For all analyses, we extracted PE-ddRAD loci in silico from the outgroups using the python script Digital_RADs. py (DaCosta and Sorenson 2014). We used PE-ddRAD IUPAC consensus sequence alignments from the selected parameterisation (Stacks higher m) to assemble two datasets containing loci present in at least 75% and 95% of the samples (PE-ddRAD 75 and PE-ddRAD 95) for downstream analyses.
Total evidence dataset. -To reduce the effect of data type on phylogenetic inference, we combined UCE and PE-ddRAD IUPAC consensus sequence alignments. We obtained a total evidence dataset with those UCE and PE-ddRAD loci present in at least 75% of the taxa (TENT 75). PE-ddRAD and UCE markers may overlap. To avoid including a sequence twice in our analyses, we used Blastn  (Altschul et al. 1997)  to map representative PE-ddRAD loci sequences to the UCE loci. We removed any PE-ddRAD loci that mapped to a UCE locus prior to concatenation.
Marker Distribution and Genomic Context Because the genomic distribution of phylogenomic markers can assist with understanding their informativeness across phylogenetic scales, we compared the distributions and genomic context of PE-ddRAD and UCE loci. We used Blastn  (Altschul et al. 1997)  to map both sets of loci to the Balearic Shearwater (BaSh; Puffinus mauretanicus) draft genome assembly  (Cuevas et al. 2020 ) and to the most closely related chromosomelevel genome assembly, the Anna's Hummingbird (AnHu; Calypte anna;  (Korlach et al. 2017) , which diverged between 62. 7 to 71. 1 Ma  (Jarvis et al. 2015). Due to the large divergence time between shearwaters and AnHu, we also mapped the PE-ddRAD markers to the AnHu genome using a liftover approach (see Supplementary Information). To determine the level of clustering of UCE and PE-ddRAD loci, and to determine their degree of association with protein-coding genes, we applied a permutation procedure (Bioconductor package regioneR;  Gel et al.
Phylogenetic Analyses We used unpartitioned UCE and PE-ddRAD concatenated datasets (UCE IUPAC 75, UCE IUPAC 95, PE-ddRAD 75, PE-ddRAD 95) and a total evidence concatenated dataset (TENT 75) partitioned by data type (UCE and PE-ddRAD) to estimate Bayesian and maximum-likelihood phylogenies using the MPI version of ExaBayes v. 1. 5 (Aberer et al. 2014 ) and raxml-ng v. 0. 6. 0  (Kozlov et al. 2019) , respectively. Additionally, we estimated PE-ddRAD ML trees including the ingroup taxon Puffinus assimilis haurakiensis, for which no UCE data were available due to sampling constraints. For each dataset, we ran two independent ExaBayes runs with four coupled chains for 1,000,000 generations. We assessed runs for stationarity in Tracer v. 1. 7  (Rambaut et al. 2018 ) by checking for effective sample sizes > 300 for all model parameters.
To verify that our datasets were not affected by these issues, we also performed analyses of UCE partitioned IUPAC consensus alignments (75% complete). We used the Sliding-Window Site Characteristics (SWSC-EN) method described in  Tagliacollo and Lanfear (2018)  and PartitionFinder 2  (Lanfear et al. 2017 ) to partition the data, which yielded 131 partitions. These analyses are explained in detail in the Supplementary Information. To account for coalescent stochasticity among individual loci, we inferred species trees with IUPAC UCE, PE-ddRAD and TENT datasets using two multispecies coalescent methods: the quartet-based method SVDQuartets  (Chifman and Kubatko 2014) , and the summary method ASTRAL-III  (Zhang et al. 2018). SVDQuartets can handle both unlinked single nucleotide polymorphisms (SNPs) and multi-locus data.
To allow a better comparison with analyses based on concatenation we ran SVDQuartets on the multi-locus sequence alignments. Analyses were run in PAUP* v. 4 (Swofford 2002) evaluating all possible quartets. For each matrix, we conducted 100 bootstrap replicates, and results were summarized in a 50% majority-rule consensus tree. For ASTRAL-III, we used RAxML v. 8  (Stamatakis 2014)  to estimate gene trees for each PE-ddRAD and UCE locus in the IUPAC 75% complete datasets. We ran 500 rapid bootstrap replicates for each individual gene followed by a thorough ML search and we estimated species trees from the best-scoring ML gene trees and bootstrap replicates using ASTRAL-III. Two analyses were run for each dataset: one using the original gene trees and one using gene trees with very low support branches (BS < 10) contracted as this procedure can improve tree accuracy  (Zhang et al.
2018). Branch support values were inferred using local posterior probabilities (PP; Sayyari and Mirarab 2016). To avoid the negative impacts of fragmentary gene sequences on gene tree and species tree reconstruction  (Sayyari et al. 2017) , we removed three samples with mean missing data values per locus higher than 5% (A. carneipes 1, A. grisea 1 and C. diomedea 2) for ASTRAL-III analyses of UCE datasets. We annotated ASTRAL-III trees with local quartet supports for the main topology, and the first and second alternatives (ASTRAL-III option -t 2) to further investigate regions of the species tree that are potentially in the anomaly zone  (Degnan and Rosenberg 2006). Finally, we performed a polytomy test  (Sayyari and Mirarab 2018)  to evaluate whether hard polytomies could be rejected at short internodes.
Divergence Time Estimation In divergence date estimation analyses, it has been common practice to reduce dataset size by selecting clock-like genes  (Smith et al. 2018). Nonetheless, we decided to perform the analyses using the three 95% complete datasets (UCE, PE-ddRAD and TENT), because recent research has shown that divergence time analyses using complete phylogenomic datasets consistently show less variance in divergence times than estimates using subsets of clock-like genes  (McGowen et al. 2019; Oliveros et al. 2019 ). We acknowledge that divergence time estimation methods based on concatenation may lead to branch-length bias and potentially misleading age estimates, particularly for younger divergence times  (McCormack et al. 2011; Angelis and Dos Reis 2015). However, as our objective was not necessarily to calculate accurate estimates but rather to compare estimates based on different phylogenomic markers, we decided to use a concatenation-based method that allowed us to use complete phylogenomic datasets.
Divergence time analyses were performed using MCMCTree v. 4. 9 from the PAML package  (Yang 2007 ). MCMCTree allows Bayesian divergence time inference of phylogenomic datasets (dos Reis and Yang 2011) by implementing approximate likelihood calculation. We pruned the topology of our ExaBayes TENT 75% tree and the TENT 75 dataset so that they contained one individual per taxon to be used as input for divergence dating analyses, retaining the most complete individual. To decide the best-fitting clock model, we used the stepping-stones method  (Xie et al. 2011)  as implemented in the mcmc3r R package  (dos Reis et al. 2018)  to calculate marginal likelihoods for relaxed-clock models using the computationally expensive exact likelihood method because the approximate likelihood method cannot be used for marginal likelihood calculation  (dos Reis et al.
We followed the two-step procedure outlined in dos Reis and  Yang (2011)  to infer divergence times using approximate likelihood calculation. For each analysis, we ran two independent Markov chain Monte Carlo (MCMC) chains, collecting 10,000 samples after a burn-in of 5,000 and a sample frequency of 500. We assessed likelihood convergence and parameters by examining trace plots in Tracer and checking that the estimated sample size (ESS) for each parameter was not smaller than 300, and by comparing results between independent runs. We also ran MCMCs with no data to generate joint prior distributions. Finally, we generated infinite-sites plots to assess how uncertainty in time estimates differed between analysis of the three datasets.
GC-biased Gene Conversion GC-biased gene conversion (gBGC) is known to strongly affect several features of avian genomes  (Nabholz et al. 2011; Weber et al. 2014). To investigate potential signatures of GC-biased gene conversion (gBGC) in base composition and substitution rates in shearwaters, SNP data from the PE-ddRAD 75 dataset were output in VCF format using the populations program in Stacks. For biallelic variant sites, we computed reference and minor allele frequencies using VCFtools v0. 1. 15 (Danecek et al. 2011). We assigned variant sites into one of the following mutation categories: strong-to-strong (S-to-S), strong-to-weak (Sto-W), weak-to-strong (W-to-S), and weak-to-weak (W-to-W), where C and G are strong bases (3 hydrogen bonds) and A and T are weak bases (2 hydrogen bonds).
Patterson's D-statistic (ABBA-BABA test). -To further explore whether tree discordances are due to past introgression or other forms of model misspecification, we quantified the Patterson's D-statistic  (Green et al. 2010; Patterson et al. 2012 ) for all species quartets compatible with the time-calibrated topology. Calculations were performed using Dsuite Dtrios  (Malinsky et al. 2020 ) with a PE-ddRAD-derived SNP dataset that included loci with data in at least 5 taxa (ddRAD_min5, 295,779 SNPs). We performed analyses for the 3 genera separately. In each analysis, the sister species to the rest of the genus was used as the outgroup (i. e. P. nativitatis, C. leucomelas, and A. bulleri and A. pacifica). Blockjackknife resampling was used to evaluate significant deviations from zero in Patterson's Dstatistic (P < 0.
We ran phylogenetic networks for each of the datasets (PE-ddRAD, UCE and TENT) and independently for Puffinus and Ardenna, to reduce computation time and to improve the accuracy of the inferred networks (Solís-Lemus and Ané 2016; Solís-Lemus et al. 2017). In both cases, we used the time-calibrated topology, pruned to only contain taxa in the genus under study, as a starting topology. For each dataset, we conducted 10 independent runs with random seeds of SNaQ to infer the optimal coalescent tree with no hybridisation edges (h0). We then performed network searches from 1 to 4 (1 to 2 for Ardenna datasets) hybridisation edges providing, in each case, the optimal network with hmax-1 hybridisation edges. The preferred number of hybridisations was selected based on the analysis of the slope of a plot of log-pseudolikelihood against the number of hybridisations (Solís-Lemus and Ané 2016; Solís-Lemus et al.
2017). We expected a sharp improvement until the number of hybridisation edges reached the best value. We did not evaluate more than four (or two in the case of Ardenna) hybridisation edges because of the lack of change in the slope heuristic. To assess whether the MSC adequately explained gene-tree discordance to our coalescent trees with no hybridisation edges (h0), we used the TICR test (Tree Incongruence Checking in R;  Stenz et al. 2015) , using the phylolm R package (Tung Ho and Ané 2014). A chi-squared test was used to compare observed concordance factors (CF) with expected CF calculated from the h0 species trees under the MSC. We further looked for taxa that did not fit the tree model with ILS retrieving the outlier 4-taxon sets. Because searches for an optimal network produced inconsistent results in independent runs and in different datasets, we focused on evaluating candidate reticulation events based on the results of the Patterson's D-statistics.
Briefly, the tendency of new polymorphic loci when increasing M or n parameters in Stacks, became linear at M=5 and n=8 (Supplementary Fig. S2 ) and the proportion of singletons in the dataset stabilised at m=7 (Supplementary Fig. S3 ). For the second optimisation step, different parameterisations had a minor effect on the pPIS with the exception of the default parameters in Stacks that yielded loci with a much lower number of PIS. With similar parameterisations, Stacks yielded approximately twice the number of loci than pyRAD for each level of missing data, but pyRAD loci had a slightly higher pPIS (Supplementary Fig. S4a  and Table  S2 ). Phylogenetic analyses using the different datasets and levels of missing data resulted in overall highly resolved and congruent phylogenies. However, the general trend was a slight increase in resolution when reducing missing data from a maximum of 35% to a maximum of 25% and thereafter, a slight decrease when reducing it to a maximum of 5% (Supplementary Fig.
S4b ). Phylogenetic analyses using datasets with a higher amount of missing data (35% and 25%) tended to yield higher bootstrap supports on recent splits, whereas analyses using datasets with a low level of missing data (5%) yielded higher bootstrap supports on more ancient splits (Supplementary Fig. S5 ). pyRAD datasets with a maximum of 25% missing data yielded the highest overall resolution but Stacks higher m datasets were the most consistent across different levels of missing data. Thus, we selected the latter parameterisation for downstream analyses. For the same taxon coverage, PE-ddRAD concatenated alignments were both longer (i. e. 2,156,937 bp for the PE-ddRAD 75 matrix vs. 1,732,076 bp for the UCE 75 matrix) and contained more than double the number of loci than UCE alignments (Table  1 ).
gavia and P. huttoni); a Subantarctic and New Zealand clade (P. elegans and P. assimilis haurakiensis); a North Pacific clade (P. newelli and P. opisthomelas); a Tropical Indian and South Pacific clade (P. bailloni); and a Caribbean, North Atlantic and Mediterranean clade (P. puffinus, P. mauretanicus, P. yelkouan, P. lherminieri, P. boydi and P. baroli) (Fig. 3a  and Supplementary Fig. S12 ). Only three phylogenetic relationships, all localised to short internodes, were challenging to resolve due to discordances between analyses using different genomic markers, levels of missing data, or phylogenetic methods (Fig. 3b-d ). The first discordance was the relationship among the three shearwater genera (triangle). All analyses using the TENT dataset and most analyses using the UCE and the PE-ddRAD datasets recovered Calonectris as sister to Ardenna (Fig.
4 ). However, raxml-ng and ASTRAL-III trees based on the PE-ddRAD 75 dataset and the SVDQuartets tree based on the UCE 95 dataset recovered the alternative topology (Ardenna as sister to Puffinus). A polytomy test based on local quartet supports (Fig. 3b ) using the TENT dataset marginally ruled out that this branch should be replaced by a polytomy (P = 0. 0324). Under a true polytomy we would expect local quartet supports for the three alternative topologies to be equal to 1/3. This test was also significant when using only UCE data but failed to discard the likelihood of a polytomy when using PE-ddRAD data (Table  2 ). The second discordance was the relationship between the three North Atlantic lineages of Puffinus (Fig. 3c ). In this case, all phylogenetic analyses recovered P. puffinus as the sister species to the Mediterranean P.
mauretanicus and P. yelkouan, with the exception of the PE-ddRAD 95 ExaBayes tree, which recovered P. puffinus as the sister to P. lherminieri, P. baroli and P. boydi. Despite consistency in the recovered relationships for this case, only the polytomy tests based on the PE-ddRAD 75 and TENT 75 datasets were able to reject the null hypothesis that the branch should be replaced by a polytomy. The last case was the relationship between the two all dark species of Ardenna (A. tenuirostris and A. grisea) (Fig. 3d ). The only analyses that recovered A. grisea and A. tenuirostris as sister species used coalescent-based methods and UCE datasets (Fig. 4 ), although it should be noted that ASTRAL-III analyses using UCEs were performed with only one A. grisea individual due to missing data filtering (See Phylogenomic Analyses in Materials and Methods).
Interestingly, all tests rejected a polytomy although UCE datasets supported one topology and PE-ddRAD datasets another (Table  2 ), showing different phylogenetic signals between both types of markers. It is also noteworthy that the monophyly of the most recently diverged taxa P. mauretanicus -P. yelkouan and A. creatopus -A. carneipes was not supported in several phylogenetic analyses (Supplementary Table  S3 ).
Divergence Dating Analysis For all sampled alignments except for one, the IR model had the highest posterior probability (Supplementary Table  S4 ). This model was interpreted as the best-fitting model and used in downstream MCMCTree analyses. Mean posterior time estimates using PE-ddRAD, UCE or TENT datasets differed only marginally (Fig. 5  and Supplementary Fig. S13 ). Analyses using the TENT dataset with two partitions produced posterior estimates with the narrowest credibility intervals while analyses using the UCE dataset produced the largest credibility intervals. The slope of the regression line in the infinite-sites plot was highest for UCE data (0. 51), was slightly lower for PE-ddRAD data (0. 49), and dropped to 0. 43 when using the TENT dataset, meaning that 0. 43 Ma of uncertainty was added to the 95% CI for every 1 Ma of divergence (Supplementary Fig.