dataesr/scientific-paragraphs-categorization huggingface dataset:
Restriction-Site-Associated DNA Sequencing Analysis 2. 5. 1. Population Differentiation Estimates among Sampling Locations Pairwise population differentiation was estimated using five different metrics: G ST , G" ST  [30] , D ST  [29] , Φ ST and F ST -like statistics based on genotype likelihood  [39]. G ST , G" ST and D ST were estimated for SNP and haplotype datasets, and Φ ST was estimated on the sequence datasets. We use G ST , G" ST and D ST due to their different sensitivity to the mutation rate, level of heterozygosity and population size  [30, 34]. G" ST , the standardized measure of G ST taking into account bias of sample size, and D ST are less sensitive to the level of heterozygosity (genetic diversity). D ST  [29]  is more sensitive to the variation of mutation rate between loci but it is considered as a better measure of differentiation in a wide range of circumstances because it is less sensitive to the population size  [29] [30] [31] 34].
To test the statistical significance of the differentiation among populations, we used Fisher's exact test, implemented in the mmod R package  [67]  using 1000 permutations for each locus. Then, the p-values of each locus were combined following Fisher's method [68] and corrected using Holm's formula for multiple testing  [69]. We also computed the pairwise genetic distances between populations by taking differences between sequences into account (Φ ST ), using K80 (Kimura 2 parameters) and TN93  [70]  substitution models. As both models provided similar results (R 2 = 0. 99), only the K80 model Φ ST s were retained. Statistical significance of the pairwise Φ ST s was evaluated using 1000 permutations for each locus and the p-values were combined using Fisher's method [68] and corrected for multiple testing according to Holm's formula  [69].
Pairwise G ST , G" ST and D ST s were estimated using the mmod R package, while Φ ST was computed with the pegas R package  [71]. The base call uncertainty caused by the low-coverage dataset was taken into account by computing a Fst-like statistic proposed by Fumagalli et al. [39] , based on genotype likelihood implemented in ANGSD using realSFS tools. The differences between various pairwise differentiation indices and different datasets in Table  1  were summarized by calculating distances and visualized with a metric multidimensional scaling approach in R. The distances were calculated as 1-r, where r is the Pearson correlation coefficient based on the pairwise distances among populations computed among indices/datasets.
Population Structure and Relationship among Groups Population structure was further inferred with both SNPs and haplotype datasets using the discriminant analysis of principal components (DAPC  [72] ) implemented in the adegenet v. 2. 0. 1 R package  [73]. The DAPCs were performed both considering the best clustering scheme inferred by the sequential K-means approach and using the sampling sites as a prior. The choice of the optimal number of clusters (K), within a range of 2 to 10, was performed according to the lowest score of the BIC criteria. The optimal number of retained principal components was performed using a cross validation procedure  [72]. The statistical support of the number of independent lineages (or populations) and their phylogenetic relationships was estimated by performing a Bayesian multi-coalescence analysis using unlinked bi-allelic SNP with SNAPP  [74]  with a Bayes Factors Delimitation (BDF*) procedure  [75, 76] , implemented in BEAST2.
SNAPP estimates the species tree using unlinked bi-allellic SNPs while by-passing the estimate of multiple gene trees. It allows estimating the statistical support of the tree topology and posterior probability of the present and ancestral effective population size assuming constant population size for each branch of the tree and sudden changes  [74]. For computational tractability, five individuals were selected per location with the least number of missing bi-allelic SNPs for the datasets with a minimum coverage of 2 and 3, respectively, and after removing loci in low frequency and potential paralogs. Three models were compared, M1: the four sampling locations as distinct entities [(Sv,Sa,Th,Kl)], M2: the mitochondrial scenario with three distinct entities [(((Sv,Sa),Th), Kl)] and M3: the ITS scenario with two distinct entities (the populations in the north vs.
We used the SNP frequency of heterozygous site to compute the D statistic considering two alternative scenarios: the mitochondrial scenario tests the introgression between Thingvallavatn (P3) and Svartárvatn (P1) or Sandur (P2) using Klapparós as an outgroup, and the ITS scenario tests the introgression between Klapparós (P3) and Svartárvatn (P1) or Sandur (P2) using Thingvallavatn as an outgroup. Switching the position of Svartárvatn and Sandur did not qualitatively affect the results. To obtain the most complete data matrix, we selected the five most informative individuals per population (maximizing the locus coverage among populations) for two datasets (with a minimal coverage of 2 and 3, respectively) and after removing low-frequency data (Cov2Hmiss_cl and, Cov3Hmiss_cl datasets). For the ingroup taxa, we selected randomly one individual for each of P1, P2 and P3 populations among the five individuals retained per population, while pooling the five individuals of the outgroup to estimate the SNP frequency  [26].
Description of the Datasets The difference in the minimal coverage used, from two to three, had a major effect on the total number of de novo loci and SNPs retained in the analysis by pyRAD shifting from 335 to 123 loci, respectively, of which 323 and 108 were polymorphic, respectively. The coverage cut-off also reduced the dataset from 3539 SNPs with two to 835 SNPs with three. The haplotype datasets contained more alleles when indels were considered as 5th state compared to when they were considered as missing data (C2, missing = 4. 5 ± 3. 3; 5th = 6. 5 ± 4. 9; paired Wilcoxon test p < 0. 001: C3, missing = 3. 8 ± 3. 2; 5th = 5. 5 ± 4. 7; paired Wilcoxon test p < 0. 001). No unlinked SNPs were identified as putatively under selection by Bayescan, while haplotype datasets revealed seven loci under putative selection using a minimum coverage of two, regardless of whether the indels are considered as missing data or not.
For a minimum coverage of three, we detected five and six loci under putative selection when indels were considered as missing data or as 5th state, respectively. Among the loci detected as putatively under selection, the Blastn approach was able to retrieve three annotated DNA fragments with an E-value statistic below 0. 001: a histone H3 nuclear gene, a mitochondrial NADH dehydrogenase subunit 5 and another mitochondrial fragment overlapping the genes ND2 and tRNA-Trp (Supplementary Materials Table  S5. 1). The NADH dehydrogenase subunit 5 was found only with the dataset with a minimal coverage of two. The proportion of loci in low frequency (with MAF < 0. 05), under putative selection and with observed heterozygosity above 0. 5 was 29. 6 ± 9. 7% for the SNPs and 28. 4 ± 3. 6% for the haplotypes.
Differentiation among Populations A large majority (i. e. , 92%) of the datasets and the differentiation metrics supported the ITS scenario (Figure  3 ). The Φ ST , in particular, was insensitive to the datasets used and supported systematically the ITS scenario (Figure  3 ). For all datasets, the significance of the pairwise G ST , G" ST and D ST between populations was always supported (p < 0. 001), even after p-value correction for multiple testing among populations. However, all the pairwise Φ ST between the southern (Thingvallavatn) and the northern populations (Sandur, Svartárvatn and Klapparós) were significant at α = 0. 05 but not among the northern populations (except between Sandur and Svartárvatn, with a minimum coverage of three, Table  S5. 2). All unlinked SNP datasets supported the ITS scenario regardless of the differentiation metrics, t coverage or the presence/absence of the loci in low-frequency or under selection (Figure  3 , Table  S5.
2). Similarly, regardless of the dataset used, the F ST -like estimates based on genotype likelihood taking into account the base call uncertainty always supported the ITS scenario. However, minimal coverage had a stronger effect on the differentiation among populations than the effect of removing the loci under low frequency at a given coverage (Figure  3 ). Furthermore, all haplotype datasets with a minimal coverage of three supported the ITS scenario regardless of the differentiation metrics used. However, with a coverage of two, considering the presence of indels as informative or not, removing loci in low frequency and under selection affected the differentiation metrics variably, as well as the inferred scenario. Two approaches showed opposite results. First, G ST , G" ST and D ST supported the ITS scenario when indels were considered as missing data and all the loci were retained, while in the second approach, considering indels as a 5th state and removing the loci in low frequency and under selection, these three metrics supported the mitochondrial scenario.
For the latter approach, 68% of the loci removed supported the ITS scenario. Finally, two haplotype datasets (C2Hmiss_cl and C2H5th_br) showed discordant results when using different differentiation metrics, with the G ST supporting the mitochondrial scenario, while the G" ST and D ST supported the ITS scenario. Such discrepancies between the G ST and the two other differentiation metrics might be explained by the greater sensitivity of G ST to the significant increase in the number of alleles and the heterozygosity within subpopulation (Hs) for all the loci supporting the discordant scenario compared to all the loci supporting the concordant scenario (Cov2H5th_br: nb loci discordant = 22, Wilcoxon test P Hs(global) = 0. 024, P nbAlleles < 0. 001; Cov2Hmiss_cl: nb loci discordant = 7, Wilcoxon test P Hs(global) = 0.
Congruence Tests among Loci Considering all datasets, 57. 4 ± 5. 0% of the loci on average supported the ITS scenario rather than the mitochondrial scenario for the G ST , G" ST and D ST. This proportion increased to 74. 9 ± 4. 0% on average when considering the Φ ST. However, no significant difference was found for the extent of differentiation between the two scenarios for any of the datasets and differentiation metrics (Wilcoxon test p > 0. 1).
Population Structure and Relationship among Groups Using SNP and Haplotype Frequencies The K-means clustering approach, based on the RADseq data, retained two clusters systematically, except for the haplotype datasets using a minimal coverage of two and after removing the loci in low frequency and under selection. For these two exceptions, the BIC inferred by sequential K-means never reached clear minima and supported the maximum number of clusters allowed (i. e. , 10). On average, the first two axes of the DAPC retained 18. 3 ± 7. 2% of the total variance. The clustering segregated the individuals into two groups, Thingvallavatn (southwestern Iceland) vs. the three populations from northern Iceland. The assignment plot shows that the level of admixture is rather weak between these two groups, and three to five individuals present in Thingvallavatn can be considered admixed depending on the dataset (Figure  4a-c ).
Congruence Tests among Loci Considering all datasets, 57. 4 ± 5. 0% of the loci on average supported the ITS scenario rather than the mitochondrial scenario for the GST, G''ST and DST. This proportion increased to 74. 9 ± 4. 0% on average when considering the ΦST. However, no significant difference was found for the extent of differentiation between the two scenarios for any of the datasets and differentiation metrics (Wilcoxon test p > 0. 1).
Regardless of the dataset, considering sampling location as the a priori number of clusters showed higher levels of admixture, especially among the three northern samples (Figure  4b-d ). (c, d ). The plots in (a,c) display the results of the optimal number of clusters inferred by a sequential K-means approach according to BIC criteria. The plots in (b,d) display the results of four sampling locations used as a priori number of clusters. The scatter plot of the DAPC is displayed as small plots and presents the main discriminants axis (1 axis for (a,c), and 2 axes for (b,d)), the ellipses represent the inertia of the clusters and the dots represent individuals. The barplots show the assignment probability of each individual to the different clusters. Th: Lake Thingvallavatn, Sv: Lake Svartárvatn, Kl: Klapparós, Sa: Sandur.
The Bayes Factor Delineation (BFD*) approach based on the bi-allelic datasets with the five most informative individuals per location favored the model of four distinct populations (M1), showing a strong support against M2 (three populations with the mitochondrial scenario) and a very strong support against M3 (two populations with the ITS scenario) for both a coverage of two and three (Table  2 ). However, the topology of the most supported four-population species tree supported clearly the scenario, with the northern populations more closely related together than with the Thingvallavatn population (Figure  5  and Figure  S5. 5). The posterior probability of the node supporting the northern clade is high (C3SNPr dataset p = 0. 97, n = 88; C2SNPr dataset p = 1, n = 223), but the node support of the clade E (from Sandur vs.
Svartárvatn) is very weak (C3SNPr dataset p = 0. 77, n = 88; C2SNPr dataset p = 0. 59, n = 223). However, the confidence intervals of the date divergence between Thingvallavatn and the three other northern populations overlap, showing a large variance in the coalescence process (Figure  5a ). The posterior distributions of θ of the populations overlapped as well and showed similar θ median between 0. 10 and 0. 13 (Thingvallavatn θ median = 0. 13, 95% HPD = 0. 08-0. 19; Klapparós θ median = 0. 10, 95% HPD = 0. 05-0. 16; Sandur θ median = 0. 12, 95% HPD = 0. 07-0. 19; Svartárvatn θ median = 0. 11, 95% HPD = 0. 06-0. 17, see Figure  5b ). [75, 76]  testing the support of three models of population structure of Crangonyx islandicus using RADseq data on unlinked bi-allelic SNPs after removing low-frequency SNPs considered as potential paralogs or putatively under selection.
Introgression Tests The introgression tests based on the D statistic showed a comparable amount of introgression among the northern populations of C. islandicus (from Sandur and Svartárvatn) with either the southern (from Thingvallavatn) or the northeastern population (from Klapparós). A test based on the conservative dataset, with a minimal coverage of three (including 21 to 45 loci) and the mitochondrial scenario with Klapparós population as an outgroup, supported an introgression between Thingvallavatn and Sandur populations, in the genetic combinations of 34 individuals for [Sv,Sa,Th][Kl] but 28 individuals when the structure of the ingroup in the test was changed [Sa,Sv,Th][Kl] (Table  3 ). Only one combination supported an introgression between Thingvallavatn and Svartárvatn populations.
The tests based on the ITS scenario (i. e. , with the Thingvallavatn population as an outgroup) showed a similar support of introgression among the northern populations with 28 or 25 combinations of individuals, indicating an introgression between Svartárvatn and Klapparós populations, but only one combination between Sandur and Klapparós populations (Table  3 ). However, when considering additional loci with a minimal coverage of two (between 62 and 128 loci were used for comparisons), a small number of combinations of individuals supported an introgression between Thingvallavatn and Sandur (five), Thing-vallavatn and Svartárvatn (four or five), Klapparós and Svartárvatn (five or seven) and Klapparós and Sandur populations (three or four, Table  3 ). population (Figures 5 and S5. 5). The posterior probability of the node supporting the northern clade is high (C3SNPr dataset p = 0.
97, n = 88; C2SNPr dataset p = 1, n = 223), but the node support of the clade E (from Sandur vs. Svartárvatn) is very weak (C3SNPr dataset p = 0. 77, n = 88; C2SNPr dataset p = 0. 59, n = 223). However, the confidence intervals of the date divergence between Thingvallavatn and the three other northern populations overlap, showing a large variance in the coalescence process (Figure  5a ). The posterior distributions of θ of the populations overlapped as well and showed similar θmedian between 0. 10 and 0. 13 (Thingvallavatn θmedian = 0. 13, 95% HPD = 0. 08-0. 19; Klapparós θmedian = 0. 10, 95% HPD = 0. 05-0. 16; Sandur θmedian = 0. 12, 95% HPD = 0. 07-0. 19; Svartárvatn θmedian = 0. 11, 95% HPD = 0. 06-0. 17, see Figure  5b ). parameters inferred by SNAPP  [73]  according to the Bayes Factor Delineation's (BFD*  [75, 76] ) best scenario for the four populations using 88 unlinked bi-allelic SNPs dataset (C3SNPr).
Node number on the MCCT displays the posterior probability of the node, and the dark blue node bar indicates the 95% HPD interval of the node age. Th: Lake Thingvallavatn, Sv: Lake Svartárvatn, Kl: Klapparós, Sa: Sandur. [73]  according to the Bayes Factor Delineation's (BFD*  [75, 76] ) best scenario for the four populations using 88 unlinked bi-allelic SNPs dataset (C3SNPr). Node number on the MCCT displays the posterior probability of the node, and the dark blue node bar indicates the 95% HPD interval of the node age. Th: Lake Thingvallavatn, Sv: Lake Svartárvatn, Kl: Klapparós, Sa: Sandur. Table  3. Results of the four-taxon D statistic testing for introgression among sampling locations of Crangonyx islandicus using RADseq data. P1, P2 and P3 are the populations under investigation for introgression and O is the outgroup.
Range Z is the range of the Z value, upper values deviated significantly from 0 (p < 0. 01). NSign/Ntot for positive D reports the number of significant replicates (combination of individuals) over the total number of tests for positive D statistic indicating an introgression between P2 and P3. NSign/Ntot for negative D indicates an introgression between P1 and P3. Range NbLoci is the range of loci used among replicates to test the D statistic. Th: Lake Thingvallavatn, Sv: Lake Svartárvatn, Kl: Klapparós, Sa: Sandur. Min. cov. presents the minimum coverage of datasets 3 (Cov3Hmiss_cl) and 2 (Cov2Hmiss_cl).
Discussion The current study showed the robustness of the mito-nuclear discordance of the phylogeographic history of C. islandicus, as newly sequenced individuals from previously described populations for mitochondrial DNA confirmed the previous mitogenomic findings by Kornobis et al. [8] , while a multitude of nuclear genomic markers sequenced through ddRADseq provided further support to the ITS scenario  [9]. A structured coalescent analysis using the mitochondrial 16S and CO1 sequences from the same four populations assessed with ddRADseq data confirmed the previous mitochondrial pattern and supports a clear pre-Ice Age divergence between the northeastern Icelandic (from Klapparós) population and the northern and southwestern populations. Despite the use of multiple ddRADseq analytical approaches designed to overcome potential misinterpretations, including different data types (SNPs, haplotypes, DNA sequences), different dataset configurations (brute vs.
[30] ). Only five out of sixty statistics calculated for the subdivision of the populations supported the mitochondrial scenario; the other supported the ITS scenario. Of the five statistics supporting the mtDNA, three were based on GST, which does not take into account the dependency of the statistic on the overall variation as G" ST and D ST do, which may thus give a biased overall estimate. Furthermore, these five methods are based on the lower coverage datasets (minimal coverage two), which included a higher number of loci but at the price of higher base call uncertainty, and four of them included indels as the 5th state. Higher coverage and omitting indels is a more conservative approach, as increasing coverage offers a more robust estimate of the DNA sequence polymorphisms and indels can lead to problems in alignment and may reflect errors in sequencing  [80].
Thus, those five approaches (dataset and metrics) supporting the mtDNA-scenario may be clearly suboptimal to assessing the population structure. Specific methods accounting for uncertainty related to low-coverage data (regardless of the minimal coverage) based on a maximum likelihood base call estimate (ANGSD) confirmed the ITS population structure scenario, also inferred by most classic approaches not accounting for such a bias. Methods relying on genotype likelihood instead of SNP calling to infer population statistics offer great promise to considerably decrease the bias caused by base call uncertainty, even for very low-coverage datasets (i. e. , x2) but remain limited to estimating basic summary statistics  [27, 39, 41, 43]. However, our results suggested that the potential noise brought by low-coverage data may not be strong enough to blur the main phylogeographic signal.
Although the genome size of C. islandicus is currently unknown, we suspect it could be large, which could partly explain the low coverage of our dataset. Indeed, crustacean and amphipod species in particular are prone to gigantism in their genome size, especially species supposedly with low metabolic rates  [82]  present in cold water environments such as the Arctic  [83, 84]  and deep waters from Lake Baikal  [85] , which share many similarities with the groundwater habitats of C. islandicus  [8, 45, 46]. Moreover, groundwater species tend to have larger genome sizes than their surface sister species due to the reduction in effective population size  [86]. When the information in a DNA sequence was reduced to a haplotype or a single SNP per locus, the results did not show any systematic bias toward a particular scenario.
The Phylogeographic Pattern of Crangonyx islandicus We have shown here that different analyses based on the extended mtDNA dataset and a larger nuclear dataset (ddRADseq) than in  [8, 9]  confirmed the different mitochondrial and nuclear phylogeographic patterns in C. islandicus. Mito-nuclear discordance among structured populations or closely related taxa has been detected in several taxonomic groups  [1, 6, 7] , including a groundwater species complex belonging to genus Crangonyx Bate, 1859  [89]. Disentangling the phylogeographic scenario and its causes remains extremely difficult, as many mechanisms can lead to the same patterns. However, from the previous and present results it is clear that this groundwater species diverged within Iceland through the Pleistocene, when Iceland was repeatedly covered by glaciers  [46, 90].
Supplementary Materials: The following supporting information can be downloaded at: https: //www. mdpi. com/article/10. 3390/d15010088/s1, Supplementary Materials S1-Sampling information; Table  S1. 1: Sampling information on 59 individuals of Crangonyx islandicus obtained for RADseq analyses and to complement mitochondrial analyses; Table  S1. 2: Information on 37 individuals of Crangonyx islandicus extracted from Kornobis et al. [1]  used to complete the MultiTypeTree structured coalescent analysis. Supplementary Materials S2-Adapter and barcode sequences; Supplementary Materials S3-Protocol of the SYBR gold fluorometric assay for dsDNA quantitation; Supplementary Materials S4-Additional analysis of the NADH and H3 loci showed indication of selection in RADseq analysis  [97] [98] [99]. Table  S4.