dataesr/scientific-paragraphs-categorization huggingface dataset:
Read Processing and Genotyping We trimmed adapters, low-quality bases, merged overlapping reads, and removed PCR duplicates from raw reads using the expHTS pipeline (v. 0. Mar112016; https://github. com/msettles/expHTS). We then used pseudo-it (v1;  Sarver et al. 2017)  to generate pseudo-reference exomes for each species by iteratively mapping (four iterations and allowing ambiguity codes) cleaned single and paired-end reads from one individual per species (Supplementary Table  S1 ) to the European rabbit reference genome (OryCun2. 0) (Carneiro et al. 2014). For snowshoe hares and black-tailed jackrabbits (L. californicus), we used pseudo-references generated by  Jones et al. (2018). The genomes of hares and European rabbits are broadly syntenic despite two known large scale fusions/fissions  (Robinson et al.
2002)  and sufficiently closely related that iterative mapping to the rabbit reference is not expected to bias genotyping  (Sarver et al. 2017; Marques et al. 2020). We then mapped data from all individuals to each species-specific pseudo-reference using bwa-mem (v. 0. 7. 12-r1039; Li 2013) with default options. Mapped reads were sorted with samtools (v1. 4; Li et al. 2009) , assigned to read groups, filtered for duplicates (Picard v1. 140; http://broadinstitute. github. io/picard/), and realigned for insertion-deletion length variation using  GATK (v3. 4. 46; Van der Auwera et al. 2013). We calculated coverage statistics and capture efficiency using CalculateHSMetrics from Picard. For each individual, we called and filtered genotypes using the bcftools (v1. 4; Li 2011) mpileup, call and filter pipeline.
2018 ) and the retention of information for phylogenetic analysis (alignment length > 1 kb). For each window, we filtered positions with missing data for > 30% of individuals using TriSeq (TriFusion 1. 0. 0; http://odiogosilva. github. io/TriFusion/), excluded windows smaller than 1 kb, and used RAxML to estimate local maximum likelihood gene trees (GTR+, 100 bootstraps). For each gene tree, we used the corresponding bootstraps to calculate a tree certainty score based on the sum of certainty scores for all internodes of a tree (-L MRE option in RAxML;  Salichos and Rokas 2013; Salichos et al. 2014; Kobert et al. 2016). The internode certainty score weighs the support of the bipartition represented by a given internode in the gene tree against the support of the second most prevalent conflicting bipartition  (Salichos et al.
2014)  present in the bootstrap trees. In our case, the maximum theoretical value of the tree certainty score is 31, or k -3 with k equal to the number of taxa  (Salichos et al. 2014 ). Only trees with certainty score above 5 were used in the species tree inference. We unrooted the gene trees using R package ape  (Paradis et al. 2004 ) and estimated a consensus species tree using ASTRAL-  III (5. 6. 3; Zhang et al. 2018). We also estimated coalescent species trees using only variable sites with SVDquartets (Chifman and Kubatko 2014) implemented in PAUP* (4a163;  Swofford 2003). For the analyses based on variable sites, we recovered single nucleotide variants (SNVs) distanced 10 kb along the genome (within targeted regions and 200 bp flanking regions) using  snp-sites (v2. 3. 3; Page et al.
2016) , a custom script and bedtools intersect, and excluded sites with missing information for >30% of the individuals using a custom script (available at https://github. com/evochange). For both SVDquartets and ASTRAL analyses, species trees were estimated with and without assigning species identities and using sites/intervals genome-wide or only from the X chromosome. The European and/or pygmy rabbits were included for all analyses requiring outgroups. Additional details on the phylogenetic analyses are provided in the Supplementary Materials and Methods.
Bayesian Divergence Time Inference We performed Bayesian inference of divergence times in the inferred species tree using an approximate maximum likelihood method and assuming an autocorrelated relaxed molecular clock, implemented in MCMCtree (PAML v. 4. 9; Yang 2007) and described in dos Reis and  Yang (2011 Yang ( , 2019)). For one individual per species (Supplementary Table  S1 ), we extracted the coding sequence for all genes included in our capture design (18,798 genes in the OryCun2. 0 ENSEMBLE 94 database) with bedtools getfasta, after selecting the longest transcript per gene using R package biomaRt (v2. 34. 2; Durinck et al. 2005 Durinck et al. , 2009)). We excluded alignments with more than 20% missing data using AMAS  (Borowiec 2016 ) (see Supplementary Materials and Methods available on Dryad for details).
Gene Tree Discordance and Phylogenetic Networks We performed a series of analyses to explore the amount and effect of gene tree discordance in our dataset. To minimize misinterpreting poor phylogenetic resolution as discordance, we retained only the filtered gene trees used as input for ASTRAL. These trees were used to estimate a SplitsTree4 network (v4. 14. 6; Huson and Bryant 2006)  with the option "Consensus Network with distances as means" and a 5% weight threshold to draw a split, and a majority rule consensus tree with RAxML (-L MRE option), including internode certainty scores  (Salichos et al. 2014). We then used DiscoVista  (Sayyari et al. 2018 ) to plot ASTRAL quartet frequencies around nodes of interest. Finally, we calculated Robinson-Foulds normalized distances between gene trees and the ASTRAL species tree where individuals are not assigned to species, and among gene trees using the function RF.
dist() from the R package phangorn (2. 4. 0; Schliep 2011). This metric varies between 0 (no discordance between trees) and 1 (complete discordance). We used PhyloNet (v3. 6. 6; Yu and Nakhleh 2015) to model species relationships under the network multispecies coalescent model, using all local genealogies (tree certainty scores > 5). Given our extensive dataset, we applied the pseudo-maximum likelihood inference of species networks  (Yu and Nakhleh 2015). We ran InferNetwork_MPL  (Yu and Nakhleh 2015)  with 0 up to 4 migration events (due to computational limitations), associating individuals to species (option -a), and optimizing branch lengths and inheritance probabilities to compute likelihoods for each proposed network (option -o). We used the best likelihoods per run to calculate BIC and AICc following  Yu et al.
(2012 Yu et al. ( , 2014) )   Given that the pseudo-maximum likelihood network inference is computationally intensive  (Yu and Nakhleh 2015) , network support could not be inferred with bootstraps. Therefore, we additionally used Treemix and 'f-branch' statistics (see below) to determine if ancestral reticulations were consistently inferred across methods. We reconstructed ancestral population graphs with TreeMix (v1. 13, options -global, -noss and -se; Pickrell and Pritchard 2012) following the steps used to generate an SNV dataset for SVDQuartets but only including Lepus species. We randomly phased the SNV alignments with a custom script (available at https://github. com/evochange) and created the TreeMix input using the script TreeMix_from_nex. py (https://github. com/mgharvey). We allowed 0 to 9 migration events and used the white-sided jackrabbit (L.
Genetic Diversity, Divergence, and Admixture We used the genomics general toolkit (https://github. com/simonhmartin/; last accessed January 14, 2019) to estimate pairwise genetic distances (d xy ) between species and nucleotide diversity () within species, and a custom script (available at https://github. com/evochange) to calculate the number of heterozygous sites per individual and the subset shared between at least one individual of each pair of species (used as proxy for shared variation). All diversity estimates were based on a genome-wide concatenated alignment, where we excluded sites with missing information for > 30% of the individuals. We then used genomics general and custom scripts (available at https://github. com/evochange) to calculate several variants of the D-statistics  (Green et al.
Whole Exome Sequencing Data Custom DNA captures showed high efficiency (32. 3 average fold-enrichment) and specificity (average 10% of sequenced bases off-target; Supplementary Table  S2 ). Mapping cleaned reads onto species-specific pseudo-references resulted in an average target sequencing coverage of 16× (5-35× on average per sample; Supplementary Table  S2 ) with 57. 8 million genotyped sites per individual (Supplementary Table  S2 ). Two lower coverage individuals (one hare individual and one pygmy rabbit individual) and one locality duplicate (one hare) were removed from the final dataset to maximize data quality and avoid geographic redundancy (Supplementary Table  S1 ). Combining new and published data, all analyses were performed on a dataset of 15 hare species (32 individuals), one pygmy rabbit, and the European rabbit reference genome, unless otherwise noted (see Supplementary Table  S1 ).
85 and 1. 99 Ma for deeper and concern species with debated species-level status (e. g. , the Corsican-broom hare and mountain-Alaskan hare complexes;  Alves et al. 2008; Melo-Ferreira et al. 2012) , while high diversity may reflect intraspecific cryptic divergence (e. g. , snowshoe hares,  Melo-Ferreira et al. 2014; Cape hares, Lado et al. 2019). On average, 49% of heterozygous sites where shared between species (Supplementary Table  S11 ). Considerable phylogenetic discordance among gene trees combined with such high levels of shared polymorphism between species could be explained by incomplete lineage sorting, secondary introgression, or a combination of both processes. We found that 88% of the minimum absolute D-statistics (D min )  (Malinsky et al. 2018 ) for all trios of species in our dataset were significantly different from zero (Bonferroni-corrected P < 0.
Genes affected by Ancestral Introgression We consistently detected ancestral introgression between the ancestors of coat color changing species in our dataset: white-tailed jackrabbits, mountain hares, Alaskan hares, and snowshoe hares (Fig. 3  and Supplementary Fig. species  (Jones et al. 2018 (Jones et al. , 2020a;; Giska et al. 2019) , we examined the contribution of gene flow to standing variation and local adaptation in these species. We detected 119 putative windows of ancient introgression across all major chromosomes, highlighting the genomewide contribution of ancestral introgression (Fig. 5 ; Supplementary Table  S13 ). The f d outlier windows of ancestral introgression contained 54 annotated genes (Supplementary Table  S14 ). This set of genes was enriched for the gene ontology term "E-Box binding" (3 of the 54 genes; Supplementary Table  S15 ), a DNA motif found in the promoters of many genes, suggesting that genomic regions affected by ancient introgression may be enriched for transcription factors involved in trans-regulation of gene expression.
2018; Giska et al. 2019 ). Here we estimated that a pulse of ancient introgression occurred at least 0. 71 Ma (Lepus fossil calibration; Supplementary Fig. S5  available on Dryad) and affected genomic regions containing genes associated with circadian rhythm regulation (ARNTL2;  Sasaki et al. 2009), pigmentation (TCF4; Furumura et al. 2001; Le Pape et al. 2009), thermoregulation (EBF2; Rajakumari et al. 2013 ) and visual perception (PDE6H;  Kohl et al. 2012). While ASIP would also be a likely candidate for adaptive introgression between these lineages  (Jones et al. 2018) , our exome sequencing coverage of this region was too sparse for detailed windowbased analysis. The functions of these introgressed genes overlap with common physiological adaptations of northern latitude animals to seasonal conditions, such as higher metabolic rates, regulation of body temperature and non-shivering thermogenesis  (Hart et al.
SUPPLEMENTARY MATERIAL Raw sequence reads used in this article are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under accession numbers SRR12020579 to SRR12020510 and BioProject number PRJNA639005. Supplementary materials, including data files and/or online-only appendices, can be found in the Dryad Digital Repository (doi:10. 5061/dryad. bzkh18967). The pipeline and custom scripts used for this study are available at https://github. com/evochange. 82x76mm (300 x 300 DPI).
A significant f b (C) value means that all descendants B of branch b share alleles with C, which is more parsimoniously explained by an event of ancestral introgression from C to b (Malinsky et al. 2018). Using the inferred species tree, we determined all conformations (A,B;C,O) needed to calculate f b (C) for all pairs of C species and b branches following Malinsky et al. (2018) with custom scripts (available at https://github. com/evochange) and R package treeman (1. 1. 3; Bennett et al. 2017). For each conformation, we calculated 'admixture proportion' (f G ) as defined in Martin et al. (2015) and Malinsky et al. (2018) and z-scores with 1 Mb block jackknife approach following Malinsky et al. (.
ACKNOWLEDGEMENTS We thank  Janet Rachlow ,  Gloria Portales ,  Armando Geraldes ,  José Carlos Brito ,  Jérôme Letty ,  Conrad Matthee ,  Alexei Kryukov ,  Ettore Randi ,  Christian Gortazar ,  Rafael Villafuerte ,  Fernando Ballesteros ,  Klaus Hackländer  and  Zbyszek Boratynski  for assistance with sample collection. The Alaskan hare samples used in this work were generously provided by the  University of Alaska Museum  (voucher codes  UAM42143  and  UAM45545 ). We thank  Fernando Seixas ,  Tiago Antão ,  Hannes Svardal ,  Nathaniel Edelman , the  EVOCHANGE lab , the  Good lab , and the  UNVEIL network  for helpful discussions on the generation and analysis of data.
Funding and support for this research was provided by  Fundação para a Ciência e a Tecnologia (FCT)  (project grant " CHANGE" -PTDC/BIA-EVF/1624/2014 ,  Portuguese National Funds)  and the  National Science Foundation (NSF)  ( EPSCoR OIA-1736249  and  DEB-1907022  grant). MSF was supported by  POPH-QREN  funds from  ESF  and  Portuguese MCTES/FCT  ( PD/BD/108131/2015  PhD grant in the scope of BIODIV PhD programme at  Faculty of Sciences, University of Porto ),  FCT  ( PTDC/BIA-EVF/1624/2014 ), and  NSF  ( OIA-1736249 ). MRJ was supported by an  NSF   Graduate Research Fellowship  ( DGE-1313190 ). JM-F was supported by an  FCT   CEEC contract  ( CEECIND/00372/2018 ). Instrumentation, laboratory, and computational support was provided by  CIBIO  NEWGEN sequencing platform, supported by  European Union 's  Seventh Framework Program  for research, technological development and demonstration under grant agreement no.
286431 , by the  University of Montana Genomics Core , supported by a grant from the  M. J. Murdock Charitable Trust. A grant from the  Eunice Kennedy Shriver National Institute of Child Health and Human Development  ( R01HD073439 ) supported the development of exome capture protocols that were utilized in the current study. Additional support was obtained from  European Union 's  Horizon 2020 Research and Innovation Programme  under grant Agreement Number  857251 , the  Laboratoire International Associé (LIA) "Biodiversity and Evolution " funded by  InEE (CNRS, France)  and  FCT (Portugal) ,  COMPETE2020 ,  PORTUGAL2020 , and  ERDF  ( POCI-01-0145-FEDER-022184 ), and from  Portugal-United States of America Downloaded  from https://academic. oup. com/sysbio/advance-article/doi/10. 1093/sysbio/syaa088/6017180 by INIST-CNRS user on 02 December 2020.
E3S was applied to MODIS/TERRA and AQUA datasets acquired during the 2005-2008 period over the mesoscale AMMA-CATCH (Analyse Multidisciplinaire de la Mousson Africaine-Couplage de l'Atmosphère Tropicale et du Cycle Hydrologique) observatory in South-West Niger. E3S estimates of instantaneous and daily available energy, evaporative fraction, and evapotranspiration were evaluated at a local scale based on two field-monitored plots representing the two main ecosystem types in the area-a millet crop and a fallow savannah bush. In addition to these ground-based observations, the local scale evaluation was performed against continuous simulations by a locally-calibrated soil-vegetation-atmosphere transfer model for the two plots. The RMSE (root mean square error) from this comparison for E3S's ET d estimates from combined AQUA/TERRA sources was 0.
Satellite and Other Mesoscale Data Table  1  describes the 1-km resolution MODIS TERRA and AQUA products used as E3S inputs, for remotely-sensed leaf area index (LAI), normalized difference vegetation index (NDVI), albedo (α), surface temperature (Ts), and emissivity (ε) of land surface (Version 5. 5, from the Land Processes Distributed Active Archive Center, https://lpdaac. usgs. gov). Temperature and emissivity were acquired every day at satellite overpass time, but they were subject to many temporal gaps due to cloud cover. To match that resolution, 4-day LAI was daily interpolated by cubic interpolation. The 8-day albedo and 16-day NDVI composites were taken as daytime values over each day within the 8-day or 16-day period. Additional information was required for the E3S method. Rainfall occurrence was obtained from the daily rainfields produced by the AMMA-CATCH observatory from its mesoscale raingauge network  [50].
Short and longwave incoming radiations were taken from the ground data at the two Wankama EC stations (next section) and applied uniformly over the mesosite. This upscaling was made possible by the limited spatial variability of incoming radiation under the predominantly clear sky conditions required for image usability. Note that the ground-based information used here could be replaced with alternative sources, such as RS or global or regional modeling. The terrain was described using the DEM GTOPO30 (http://earthexplorer. usgs. gov/). It was used to filter out MODIS pixels based on maximum slope and elevation criteria, albeit with no effect in this study area. Finally, the global land evaporation Amsterdam model (GLEAM) global ET products were used for mesoscale comparison with E3S results  [47].
Local Scale Data The two reference local plots in Wankama are equipped with EC systems, measuring in particular high frequency (20Hz) 3D wind, air temperature, and water vapor, as well as four-component radiation and soil heat flux (Table  2 ). A full description of all sensors, setup, and other measured variables can be found in Ramieret al. [52]  or Velluetet al. [4]. Turbulent fluxes, including ET, were derived as described by Ramieret al. [52] , using EdiRe software (R. Clement, University of Edinburgh) and following CarboEurope recommendations  [53]. All water and heat fluxes were aggregated to 30-min timestep. Mean energy balance residuals for raw daily daytime fluxes were 3. 6 and 14. 2 W•m -2 at the fallow and millet sites, respectively  [52]. They were attributed, at least, in part to unaccounted-for positive daytime storage in the shallow soil, for lack of shallow soil temperature measurements.
EC measurements were not corrected for energy balance closure as residuals were reasonably small and could be attributed to a variety of sources. Because these field measurements inevitably included a significant fraction of gaps in the various time series, using the datasets directly for comparison with the RS-based estimates had limitations. This was especially true for daily-integrated fluxes because of the high dependency of the EC method upon aerologic conditions over the day, which seriously reduced the number of days available for direct comparison. Therefore, to allow for a more complete evaluation, the results of the E3S method were also compared with simulations by the physically-based SiSPAT land surface model over 2005-2008. These simulations might be considered as reliable proxies of the local plots fluxes, as the model was thoroughly calibrated and validated by Velluetet al.
[4]  against the multi-variable dataset from the two EC plots for the 7 years 2005-2012. The model was shown to perform very well against the observations over the wide range of conditions encountered during those 7 years, whether in the calibration or the validation sub-period. Table  3  displays the model's skill scores calculated for the present study period  (2005) (2006) (2007) (2008) , for 30-min net radiation, soil heat, and latent heat fluxes. All these values were highly similar to those obtained for the entire model calibration and validation periods (see Table  B1  in Velluetet al. [4] ). The excellent performance of the model for the very complete set of water and energy variables supported its use as a trustworthy extrapolator of field data for E3S evaluation. Further, the model was suggested as providing a reliable alternative to measurement correction techniques (such as energy balance closure enforcement), or estimation of non-directly measured variables, such as evaporative fraction and surface soil heat flux (at zero depth).
Acknowledgments: The ground data set was obtained in the framework of the  AMMA (African Monsoon Multidisciplinary Analysis) Program  and the  AMMA-CATCH  observatory (www. amma-catch. org),with funding from  IRD ,  CNRS-INSU, OREME , and HSM. Special thanks are due to the local IRD team in Niger as well as to  F. Arpin-Pont ,  H. Barral ,  N. Boulain ,  G. Charvet , and  D. Ramier  for the ground data acquisition and processing. Provision of the MODIS products by the  Land Processes Distributed Active Archive Center (LP DAAC)  and of the GLEAM evapotranspiration products by the  Global Land Evaporation Amsterdam Model team  is gratefully acknowledged. The authors also wish to thank the three anonymous reviewers for their valuable comments.
Funding: The first author's Ph. D. was financed by a student research grant from the  GAIA Doctoral School at the University of Montpellier  (https://gaia. umontpellier. fr/). This study also received support from the  PNTS program  of the  French Centre National de Recherche Scientifique  ( PRETAO  project -PRoduits d'EvapoTranspiration en Afrique de l'Ouest) , as well as from  CNES  ( French Centre National d'Études Spatiales)  through the  TOSCA  projects  EVASPA V3. 0  and  PITEAS (Production par Interpolation Temporelle de l'Evapotranspiration journalière issue de l'infrarouge thermique à l'Aide de données auxiliaireS ), and from the  CRITEX program  of the  Agence Nationale de la Recherche  (grant  ANR-11-EQPX-0011 ).
webMCP-counter Our web interface, webMCP-counter, is composed of two main steps (Figure  1 ). The first step aims at running (m)MCP-counter per se. This step requires the user to upload a transcriptomics dataset that can be presented under various formats. The user also needs to specify several details regarding the dataset (file format, gene ID format and organism). Either HUGO gene symbol or ENSEMBL ID are valid formats for genes. Depending on the organism, either MCP-counter (human data) or mMCP-counter (mouse data) is run, the scores are presented as a table that can be visualized within the tool or downloaded for other purposes. The second step consists in a typical downstream analysis. The user is prompted to choose a number of clusters and the cell populations to include in the analysis among the populations available, which depends on the organism and the data provided.