dataesr/scientific-paragraphs-categorization huggingface dataset:
Funding: This research was funded by  DFG  through  SPP1175 (D. L. LI 621/4-1 ,  LI 621/4-2 ; D. C. L. ) and  Individual Grant  ( D. L. LI 621/3-3 ,  LI 621/10-1; D. C. L. LA 1971/6-1 ), the  German Excellence Initiative via "Nanosystems Initiative Munich (NIM) " (D. C. L. ) and the  Ludwig-Maximilians-University Munich  ( LMU Innovativ BioImaging Network , D. C. L. ). I. M. G. gratefully acknowledges the financial support of the  Studienstiftung des deutschen Volkes  for a PhD fellowship.
It is well known that geometric morphometric data cannot be directly used for testing general linear models because the covariance matrix of shape variables is singular (which, in turn, is due to the loss of dimensions when removing the effect of rotation, translation and scale in GPA). To circumvent this problem, we used both a traditional, and conservative, approach and a recently developed approach. In particular, the traditional approach consists in performing a principal component analysis (PCA) and retaining only axes with nonzero eigenvalues, using the scores on these axes as dependent variables in the linear model. The second approach  (Clavel et al. 2019) , implemented in the R package mvMORPH  (Clavel et al. 2015) , uses a penalized likelihood framework to circumvent the problem of singularity of covariance matrices.
All the analyses of disparity were performed using the R package GeometricMorphometricsMix (Fruciano 2020). In particular, we used permutation tests to compare the disparity between groups employing two test statistics: multivariate variance and mean pairwise Euclidean distances. These are statistics widely used to characterize the amount of variation in both traditional and geometric morphometric samples  (Ciampaglio et al. 2001; Kaplan 2004; Fruciano et al. 2014; Fruciano et al. 2016). In addition to these, we also computed in Geometric Morpho-metricsMix rarefied estimates of the following statistics: multivariate variance, mean pairwise Euclidean distances, proper variance  (Claramunt 2010) , and multidimensional convex hull volume. All of these, except proper variance, are widely used measures of disparity and computing race-specific rarefied estimates and confidence intervals allow further comparison between groups of the amount of dispersion around mean shape.
e. for each pair of axes, correlation between left PLS scores and right PLS scores) and singular values  (Rohlf and Corti 2000). Next, we tested whether the patterns of covariation detected by PLS (i. e. how shape of the two parts co-vary) were different between the two races. To this aim, as PLS identifies pairs of axes (one per module; called "singular axes"), we performed in GeometricMorphometricsMix an angular comparison of the first PLS axis (singular axis) obtained for the 2n ¼ 40 race with the corresponding axis obtained for the 2n ¼ 22 race. We did this analysis for both the data with a single GPA and the data with separate GPAs and used for the testing the method in Klingenberg and Marugan-Lobon (2013), which uses the formulas in  Li (2011). With this approach, a significant test implies that the angle is small enough that the two multivariate vectors being tested are "significantly similar," whereas a non-significant test suggests that the vectors have different directions (i.
As Xi is not a symmetric coefficient (Chatterjee 2020), we computed it in the R package XICOR (Chatterjee 2020) considering Procrustes distance from standard integration as dependent on diploid number (as per our biological hypothesis). To test whether the association detected by PLS analysis could be simply explained by allometry, we performed a further angular test. For this test, we considered the case with separate GPAs and for each module we tested whether the vector of coefficients of a multivariate regression of shape on size for all specimens (including also individuals with intermediate karyotype) was different from the corresponding vector of PLS coefficients for the standard race. This means asking whether the shape change due to allometry is approximately the same as the shape change due to "normal" co-variation between modules.
These E1BSs are often regarded as the 'origin of DNA replication'. The E2 viral protein is an essential transcription regulator that binds specifically to 12 bp motifs-E2binding sites (E2BSs)-located mostly within the URR (McBride 2013). In addition, E2 modulates the shift from early to late transcript production, acting independently on E2BS outside the URR  (Johansson et al. 2012). In the coding regions, the CUPrefs of a virus may shift after drastic genomic changes. Since PVs depend on the host translation machinery and on the available host tRNAs, one can expect that viruses would evolve to match their CUPrefs to those of the host. Therefore, proteins required in large amounts are usually encoded by genes optimized to the host's CUPrefs while a poor match of CUPrefs generally results in lower protein production  (Bahir et al.
PV genome sequences and their characteristics We collected the complete genomes of PVs infecting Cetartiodactyla from the Papillomavirus Episteme database (PaVE: https://pave. niaid. nih. gov/) and GenBank (https://www. ncbi. nlm. nih. gov/genbank/) between March and May 2018. The ORFs (E10, E6, E7, E1, E2, E4, E5, L2, and L1) and the URR of 58 reference PV genomes were extracted for subsequent analyses. For each PV genome, we collected information on the corresponding host species, the clinical presentation of the infection, the anatomical location, and viral taxonomy, as reported by the authors in the corresponding PaVE and GenBank entries or publications (Supplementary Table  S1 ). The PVs in this study were sampled from twenty different host species that belong to seven distinct host families (Bovidae, n ¼ 30; Camelidae, n ¼ 2; Cervidae, n ¼ 11; Delphinidae, n ¼ 8; Giraffidae, n ¼ 1; Phocoenidae, n ¼ 4; Suidae, n ¼ 2).
The detected motifs were also compared with the known regulatory motifs in PVs  (Bergvall, Melendy, and Archambault 2013; Bernard 2013; McBride 2013) , as well as with those in the online databases TOMTOM  (Gupta et al. 2007 ) and TRANSFAC  (Wingender 2000). For certain PV genomes for which important motifs were not detected in the URR, we used FIMO implemented in the MEME Suite to scan for the presence of these motifs elsewhere in the genome.
The algorithm in this program allows us to compare the CUPrefs of a query (ORFs of PV genomes) to those of a reference data set (ORFs of host genomes) and outputs a normalized value. The COUSIN score can be interpreted as follows: COUSIN ¼ 1, the CUPrefs of the PV ORFs are similar to those of the corresponding host; COUSIN ¼ 0, the CUPrefs of the PV ORFs are similar to a random usage of synonymous codons; COUSIN < 0, the CUPrefs of the PV ORFs are opposite to those of the corresponding host (i. e. the less used codons in the host reference are used more often in the query than in the null hypothesis of equal frequency), and COUSIN > 1, the CUPrefs of the PV ORFs are superior to those in the reference (i. e. the more frequent codons in the host reference are even more frequently used in the query)  (Bourret et al.
2019). To calculate the CUPrefs of the hosts, a representative genome for each host family was chosen and the respective codon usage tables were calculated. The representatives used are: Bovidae-Bos Taurus (accession: AC_000158), Camelidae-Camelus dromedarius (accession: NW_011590949), Cervidae-Odocoileus virginianus (accession: NW_018326927), and Suidae-Sus scrofa (accession: NC_010443). For PVs infecting Delphinidae and Phocoenidae, we chose a common representative, Tursiops truncatus (accession: NW_017842062), as both host families are closely related. We did not calculate the CUPrefs for the Giraffidae family as the available giraffe genomes (GenBank accessions: LVKQ00000000. 1 and SJXV00000000. 1) are not annotated, hence GcPV1 was removed from the COUSIN analysis.
Phylogenetic reconstruction of PVs infecting Cetartiodactyla We collected 58 PVs infecting Cetartiodactyla from the PaVE and GenBank databases (Supplementary Table  S1 ). ML phylogenetic trees of the concatenated early genes (E1-E2) and the concatenated late genes (L2-L1) were constructed at the nucleotide (Fig. 1 ) and amino acid (Supplementary Fig. S1 ) levels. The constructed trees are well supported with high bootstrap values, although few inner branches have low (>30 and <50) bootstrap values. In both E1-E2 and L2-L1 phylogenetic trees at nucleotide and amino acid levels, the Delta-Zeta crown group (blue in Fig. 1  and Supplementary Fig. S1 ) forms a monophyletic clade. The other crown groups, Alpha-Omikron (coloured orange/red) and Beta-Xi (coloured green), and unclassified PVs (coloured purple), form monophyletic clades in the early gene trees.
S2 ). The E2BS was detected 298 times in 56 out of 58 sequences, as we could not detect this motif in the URR of BPV19 and BPV21. The E1BS was detected sixty times in all fifty-eight sequences. Since the E1-and E2BSs are pivotal for the PV life cycle, we suspect that the E2BS is located elsewhere in the genome of PVs lacking these in the URR. Indeed, for both BPV19 and BPV21, we detected an E2BS within the L2 gene. Moreover, it is likely that additional E1-and/or E2BSs were not detected due to sequence divergence from the consensus motif sequence of PVs included in this study. Apart from the E2BS and E1BS, we detected twenty other motifs. However, we were not able to match these motifs with other known PV regulatory motifs or with those known in the online databases (TOMTOM and TRANSFAC). Certain of the URR motifs seemed to be exclusive to specific PVs; motifs M8, M9, M10, and M15 are present only in recombinant PVs, while motif M6 is solely found in a smaller group of Beta-XiPVs (Fig.
Mitochondrial diversity within Carnivora The first dataset, composed of complete carnivoran mitogenomes available in GenBank combined with the newly generated sequences of the two subspecies of P. cristata, the two subspecies of O. megalotis, Parahyaena brunnea, Speothos venaticus and Vulpes vulpes, plus the sequences extracted from UCE libraries for Bdeogale nigripes, Fossa fossana, and Viverra tangalunga, consists of 142 species or subspecies representing all families of Carnivora. Maximum likelihood (ML) analyses reconstructed a robust mitogenomic phylogeny, with 91. 4% of the nodes (128 out of 140) recovered with bootstrap support higher than 95% (Fig. 2a ). The patristic distances based on complete mitogenomes between the allopatric subspecies of aardwolf and bat-eared fox were 0. 045 and 0.
020 substitutions per site, respectively (Table  S1 ). These genetic distances are comparable to those observed between different well-defined species of Carnivora such as the red fox (Vulpes vulpes) and the fennec (Vulpes zerda) (0. 029) or the Steppe polecat (Mustela eversmannii) and the Siberian weasel (Mustela sibirica) (0. 034) (see Table  S1 ). To further assess the genetic distances between the two pairs of subspecies and compare them to both polymorphism and divergence values observed across Carnivora, two supplemental datasets including at least two individuals per species were assembled by retrieving all COX1 and CYTB sequences, which are the two widely sequenced mitochondrial markers for carnivores, available on GenBank. These datasets include 3,657 COX1 sequences for 150 species and 6,159 CYTB sequences for 203 species of Carnivora.
After adding the corresponding sequences from the newly assembled mitogenomes, ML phylogenetic inference was conducted on each dataset (Supplementary materials). The patristic distances between all tips of the resulting phylogenetic trees were measured and classified into two categories: (i) intraspecific variation (polymorphism) for distances inferred among individuals of the same species and (ii) interspecific divergence for distances inferred among individuals of different species. Despite an overlap between polymorphism and divergence in both mitochondrial genes, this analysis revealed a threshold between polymorphism and divergence of apprrowimately 0. 02 substitutions per site for Carnivora (Fig. 2b ). With a nucleotide distance of 0. 054 for both COX1 and CYTB, the genetic distance observed between the two subspecies of aardwolf (Proteles ssp.
Phylogenomics of Carnivora Phylogenetic relationships within Carnivora were inferred from a phylogenomic dataset comprising 52 carnivoran species (including the likely new Proteles septentrionalis species) representing all but two families of Carnivora (Nandiniidae and Prionodontidae). The nonannotated genome assemblies of these different species were annotated with a median of 18,131 functional protein-coding genes recovered for each species. Then, single-copy orthologous gene identification resulted in a median of 12,062 out of the 14,509 single-copy orthologues extracted from the OrthoMaM database for each species, ranging from a minimum of 6,305 genes for the California sea lion (Zalophus californianus) and a maximum of 13,808 for the dog (Canis lupus familiaris) (Table  S5 ). Our new hybrid assemblies allowed the recovery of 12,062 genes for the Southern aardwolf (P.
Indeed, estimation of the genetic differentiation either within well-defined species (polymorphism) or between two closely related species (divergence) allowed us to define a threshold between genetic polymorphism and genetic divergence across Carnivora (Fig. 5 ). Given these estimates, and in accordance with mitochondrial data, the two subspecies of P. cristata (1) present more genetic differentiation between each other than the two well-defined species of golden jackal (Canis aureus) and wolf (Canis lupus), and (2) present more genetic differentiation than the more polymorphic species of the dataset, the lion (P. leo). Despite known cases of natural hybridization reported between C. aureus and C. lupus  (Galov et al. , 2015; Gopalakrishnan et al. , 2018) , the taxonomic rank of these two species is well accepted.
In that sense, given the species used as a reference, the two subspecies of P. cristata seem to deserve to be elevated to the species level. The situation is less clear regarding the subspecies of O. megalotis. Indeed, while the genetic differentiation observed between the two subspecies is significantly higher than the polymorphic distances observed for all the well-defined species of the dataset, there is no species in our dataset that exhibits equivalent or lower genetic divergence than a closely related species. This illustrates the limits of delineating closely related species due to the continuous nature of the divergence process  (De Queiroz, 2007). The subspecies of O. megalotis fall into the "grey zone" of the speciation continuum  (De Queiroz, 2007; Roux et al. , 2016)  and are likely undergoing speciation due to their vicariant distributions.
To be congruent with the genetic divergence observed across closely related species of Carnivora (according to our dataset), we thus propose that (1) the taxonomic level of the P. cristata subspecies be reconsidered by elevating the two subspecies P. c. cristata and P. c. septentrionalis to the species level, and (  2 ) the taxonomic level for the two subspecies of O. megalotis be maintained. These new taxonomic results should prompt a deeper investigation of morphological and behavioural differences that have been reported between the two proposed subspecies of aardwolf to formally validate our newly proposed taxonomic arrangement. They also have conservation implications, as the status of the two distinct aardwolf species will have to be re-evaluated separately in the International Union for Conservation of Nature (IUCN) Red List of Threatened Species  (IUCN 2020 (IUCN , 2020)).
Mitogenomic dataset construction In order to assemble a mitogenomic data set for assessing mitochondrial diversity among P. cristata and O. megalotis subspecies, we generated seven new Carnivora mitogenomes using Illumina shotgun sequencing (Table  S6 ). Briefly, we extracted total genomic DNA total using the DNeasy Blood and Tissue Kit (Qiagen) for P. c. cristata (NMB TS307), P. c. septentrionalis  (NMS Z. 2018. 54), O. m. megalotis (NMB TS305), O. m. virgatus (FMNH 158128) , Speothos venaticus (ISEM T1624), Vulpes vulpes (ISEM T3611), and Parahyaena brunnea (ISEM FD126), prepared Illumina libraries following the protocol of Tilak et al. (Tilak et al. , 2015) , and sent libraries to the Montpellier GenomiX platform for single-end 100 bp sequencing on a Illumina HiSeq 2500 instrument to obtain about 5 to 10 million reads per sample.
We then assembled and annotated mitogenomes from these single-read shotgun sequencing data with MitoFinder v1. 0. 2  (Allio et al. , 2020a)  using default parameters. We also used MitoFinder to extract three additional mitogenomes from paired-end Illumina capture libraries of ultra-conserved elements (UCEs) and available from the Short Read Archive (SRA) of NCBI for Viverra tangalunga, Bdeogale nigripes, and Fossa fossana Additional read mappings were done with Geneious  (Kearse et al. , 2012)  to close gaps when the mitochondrial genome was fragmented. Finally, we downloaded all RefSeq carnivoran mitogenomes available in Genbank (135 species as of July 1 st , 2019) and the mitogenome of the Malayan pangolin (Manis javanica) to use as outgroup.
Mitogenomic phylogenetics and dating Mitochondrial protein coding genes were individually aligned using MACSE v2  (Ranwez et al. , 2018)  with default parameters, and ribosomal RNA genes using MAFFT  (Katoh and Standley, 2013)  algorithm FFT-NS-2 with option --adjustdirection. A nucleotide supermatrix was created by concatenating protein-coding and ribosomal RNA genes for the 142 taxa (140 species and 2 subspecies). Phylogenetic inferences were performed with Maximum likelihood (ML) as implemented in IQ-TREE 1. 6. 8 (Nguyen et al. , 2014)  with the GTR+G4+F model. Using the resulting topology, divergence time estimation was performed using Phylobayes v4. 1c  (Lartillot et al. , 2013)  with strict clock (CL), autocorrelated (LN or TK02), and uncorrelated (UGAM or UCLM) models combined with 18 fossil calibrations (Table  S7 ).
Three independent Markov chains Monte Carlo (MCMC) analyses starting from a random tree were run until 10,000 generated cycles with trees and associated model parameters sampled every cycle. A burn-in of 25% was applied before constructing the majority-rule Bayesian consensus tree with the readdiv subprogram. Finally, to determine the best-fitting clock model, cross-validation analyses were performed with Phylobayes by splitting the dataset randomly into two parts. Then, parameters of one model were estimated on the first part of the dataset (here representing 90%) and the parameter values were used to compute the likelihood of the second part of the dataset (10%). This procedure was repeated ten times for each model. Finally, the likelihood of each repeated test was computed and summed for each model with the readcv and sumcv subprograms, respectively.
Mitochondrial diversity and barcoding gap analyses To check if a threshold between intraspecific variation and interspecific divergence could be determined across Carnivora  (Meyer and Paulay, 2005) , two mitochondrial barcoding datasets were assembled from all COX1 and CYTB sequences available for Carnivora plus the corresponding sequences for the two subspecies of O. megalotis and P. cristata, respectively. After aligning each barcoding dataset with MACSE v2, ML phylogenetic inferences were performed with IQ-TREE 1. 6. 6 using the optimal substitution model as determined by ModelFinder  (Kalyaanamoorthy et al. , 2017). Then, pairwise patristic distances between all individuals were calculated from the resulting ML phylogram. Finally, based on the actual taxonomic assignment, patristic distances were considered as intraspecific variation between two individuals belonging to the same species and as interspecific divergence between individuals of different species.
Hybrid assembly of short and long reads Short reads were cleaned using Trimmomatic 0. 33  (Bolger et al. , 2014)  by removing low quality bases from their beginning (LEADING:3) and end (TRAILING:3), by removing reads shorter than 50 bp (MINLEN:50). Quality was measured for sliding windows of four base pairs and had to be greater than 15 on average (SLIDINGWINDOW:4:15). For MinION sequencing, base calling of fast5 files were performed using Guppy v3. 1. 5 (developed by ONT) with the high accuracy option, which is longer but more accurate than the standard fast model (Fig. S1 ). Long read adapters were removed using Porechop v0. 2. 3 (https://github. com/rrwick/Porechop). To take advantage of both the high accuracy of Illumina short reads sequencing and the size of MinION long reads, assemblies were performed using the MaSuRCA hybrid genome assembler  (Zimin et al.
Then a global GDI comparing the private polymorphisms of individuals from the two populations (p within AB) and the total polymorphism of the species (the two populations, p between AB) was estimated with one individual from each population (Fig. S2 ). Finally, the two GDI were compared to check if the Southern populations were more structured than the entire populations. To contextualize these results, the same GDI measures were estimated for welldefined species of Carnivora. The species pairs used to make the comparison and thus help gauging the taxonomic status of the bat-eared fox and aardwolf subspecies were selected according to the following criteria: (1) the two species had to be as close as possible, (2) they had both reference genomes and short reads available, (3) their estimated coverage for the two species had to be greater than 20x, and (4) short read sequencing data had to be available for two individuals for one species of the pair.
Demographic analyses Historical demographic variations in effective population size were estimated using the Pairwise Sequentially Markovian Coalescent (PSMC) model implemented in the software PSMC (https://github. com/lh3/psmc) (Li and Durbin, 2011). As described above, cleaned short reads were mapped against the corresponding reference genome using BWA-MEM  (Li, 2013)  and genotyping was performed using Freebayes v1. 3. 1-16 (git commit id: g85d7bfc)  (Garrison and Marth, 2012)  for the three individuals of each species. VCF files were converted to fasta format using a custom python script, excluding positions with quality below 20 and a depth of coverage below 10x or higher than 200x. Diploid sequences in fasta format were converted into PSMC fasta format using a C++ program written using the BIO++ library  (Guéguen et al.
, 2013)  with a block length of 100bp and excluding blocks containing more than 20% missing data as implemented in "fq2psmcfa" (https://github. com/lh3/psmc). PSMC analyses were run for all other populations testing several -t and -p parameters including -p "4+30*2+4+6+10"  (Nadachowska-Brzyska et al. , 2013)  and -p "4+25*2+4+6"  (Kim et al. , 2016)  but also -p "4+10*3+4", -p "4+20*2+4" and -p "4+20*3+4". Overall, the tendencies were similar but some parameters led to unrealistic differences between the two individuals from the South African population of Otocyon megalotis. We chose to present the results obtained using the parameters -t15 -r4 -p "4+10*3+4". For this parameter setting, the variance in ancestral effective population size was estimated by bootstrapping the scaffolds 100 times.
To scale PSMC results, based on several previous studies on large mammals, a mutation rate of 10 -8 mutation/site/generation  (Ekblom et al. , 2018; Gopalakrishnan et al. , 2017)  and a generation time of two years  (Clark, 2005; Koehler and Richardson, 1990; van Jaarsveld, 1993)  were selected. Results were plotted in R v3. 63 (Team, 2020) using the function "psmc. results" (https://doi. org/10. 5061/dryad. 0618v/4) (Liu and Hansen, 2017)  modifed using ggplot2 (Wickham, 2016) and cowplot  (Wilke, 2016).