dataesr/scientific-paragraphs-categorization huggingface dataset:
We assigned a function to each gene following a BLAST (Basic Local Alignment Search Tool) analysis using known databases and comparison with the literature. We focused on the genes that followed the main expression pattern in response to poly(I•C), i. e. , a sustained expression profile. We restricted this analysis to genes specific to poly(I•C) response in excluding genes whose expression was also found sustained under the control condition [Table  S5D ]). We found that 43% of those genes could be assigned to an antiviral response ("V"), 5% to an antibacterial response ("B"), and 20% potentially to both ("V/B"), while 32% could not be directly associated with known immune pathways ("other") (Fig. 6C ). In the sustained expression profile, we then assigned the genes from the "V" and "V/B" categories to known antiviral pathways (Fig.
5 months; shell length, 20. 6 Ϯ 0. 1 mm; total weight, 1. 62 Ϯ 0. 1 g; means Ϯ standard errors [SE]) of the Pacific oyster Crassostrea gigas (or Magallana gigas  [65] , Thunberg 1793). Oysters were offspring of the F 11 family of the DECIPHER program (ANR-14-CE19-0023) produced in July 2016. These oysters were chosen for their high susceptibility to juvenile oyster syndrome  (14). As priming solution, poly(I•C) high molecular weight (HMW) (InVivogen; catalog code tlrl-pic) was prepared following the manufacturer's instructions and diluted in filtered seawater (FSW) to a concentration of 475 ng l Ϫ1. As a control, a solution of autoclaved and filtered seawater at 0. 02 m was used. For the priming experiment, oysters were anesthetized in hexahydrate MgCl 2 (ACROS; catalog number 197530250, 50 g liter Ϫ1 , 100 oysters liter Ϫ1 ) according to the method of Suquet et al.
Virus-specific primers targeting a region of the OsHV-1 genome predicted to encode a DNA polymerase catalytic subunit (DP) were chosen as previously described (ORF100; nucleotides 147655 to 153291 of the OsHV-1 genome [GenBank accession number AY509253]: OsHVDPFor, ATTGATGATGTGG ATAATCTGTG; OsHVDPRev, GGTAAATACCATTGGTCTTGTTCC)  (67). For absolute quantification, DP amplification products were cloned into the pCR4-Topo vector and replicated in Escherichia coli DH5␣ (Invitrogen). Plasmids were extracted using the Wizard Plus SV miniprep DNA purification system (Promega) and standard curves of known concentration of plasmid generated according to the Applied Biosystems manual of absolute real-time RT-PCR quantification  (77). Absolute quantification of OsHV-1 genome copies in oyster samples was then estimated by comparing the observed Cp values to known plasmid standards from 10 3 (limit of detection) to 10 9 copies of DP.
Viral load data, expressed in copies of DP per nanogram of total gDNA, were transformed with the neperian logarithm function y ϭ ln ͑x ϩ 1͒ and normality (D'Agostino and Pearson omnibus normality test, P value Ͻ 0. 05) and variance homoscedasticity (Brown-Forsythe test, P value Ͻ 0. 05) were tested. The effect of priming on OsHV-1 load was determined using one-way analysis of variance (ANOVA) using the computer software package GraphPad Prism version 6. 0. Bonferroni's correction for multiple-comparison tests was used to compare the viral load from each postchallenge sample to 10 DPP. RNA sequencing. Quality check of total RNA was performed by capillary electrophoresis (Agilent BioAnalyzer 2100) and showed no degradation. mRNA purification and library construction were performed by MGX (Montpellier Genomix, Institut de Génomique fonctionnelle, Montpellier, France).
Equimolar pools including 5 or 6 libraries were generated, and sequencing was performed on a HiSeq2500 using the single-read protocol (50 nucleotides [nt]) on a total of 9 lanes using v4 reagents. RNA-seq analysis. All data treatments were carried out on a local Galaxy server (https://bioinfo. univ -perp. fr) (68). Read quality was checked by FastQC software (version 0. 11. 5), and the high quality observed (Phred score Ͼ32) allowed us to keep all sequenced reads for further analyses. The sequencing depth was between 19 and 40. 3 million reads per sample (Table  S1 ). For OsHV-1 transcriptome analysis, reads were aligned to the genomes of OsHV-1 variant A (GenBank accession number KY242785) and variant B (GenBank accession number KY271630) (69) using Bowtie2 (Galaxy version 2. 3. 0. 1) (70) (Table  S1 ).
For oyster transcriptome analysis, mapping on the C. gigas reference genome (version V9  [20] ; (GenBank accession number AFTI00000000) was performed using RNA STAR (Galaxy version 2. 4. 0d-2) under default parameters  (72). The read count was performed using HTSeq-count (Galaxy version 0. 6. 1), in union mode, with a minimum alignment quality of 10 and feature type of mRNA  (71). The gff file used as input was provided by Zhang et al. (20). DEseq2 (Galaxy version 2. 11. 39) was used to analyze the differential gene expression level between treatments  (73). Fold changes between treatments were considered significant when the adjusted P value with the false-discovery rate (FDR) (Benjamini-Hochberg procedure  [74] ) was Ͻ0. 05. Gene ontology analysis. In order to identify biological functions affected by poly(I•C) priming, enrichment analyses were performed by a gene ontology (GO) analysis using a rank-based statistical test (Mann-Whitney U test) combined to an adaptive clustering  (75).
Analyses were done with an R and Perl script available at https://github. com/z0on/GO_MWU (76). The following parameters were used: largest ϭ 0. 1, smallest ϭ 10, and clusterCutHeight ϭ 0. 25. To generate the table of continous measures, from the 28. 027 genes of the oyster genome, the log 2 fold change was assigned to genes significantly differently expressed (DEseq2; FDR Ͻ 0. 05) and "0" to genes that were not significantly differentially expressed. A GO category was considered enriched for an FDR of Ͻ0. 01. To synthetically represent the GO results, a ratio was calculated for upregulated and downregulated enriched categories with the following calculations: (i) for upregulated enriched categories, (number of genes significantly upregulated in the category/total number of genes in this category) or (ii) for downregulated enriched categories, Ϫ1 ϫ (number of genes significantly downregulated in this category/total number of genes in this category).
Results were represented on a heat map using a hierarchical clustering with uncentered Pearson correlation (MeV 4_9_0 software). Analyses of transcript differential expressions along the priming and challenge kinetics. In order to characterize the effect of priming and challenge on the expression pattern of candidate genes through time, DEseq2 was used to calculate log 2 fold change expression on the set of genes present in the enriched GO categories previously obtained. Regarding postpriming [poly(I•C) or FSW], each time point (i. e. , 0. 5, 1, and 10 days) was compared to T0. Regarding postchallenge with OsHV-1, each time point (i. e. , 0. 5 and 1 day) was compared to 10 DPP. To precisely determine the effect of OsHV-1, a third DESeq2 analysis was performed postchallenge by comparison between OsHV-1 and control inoculum treatments under poly(I•C)-and FSW-primed conditions.
Results were represented on a heat map using Microsoft Excel 2010. To investigate the putative functions of the genes falling into the different expression patterns, we used the functional annotations of those genes. To work with current functional annotations of the C. gigas gene set, we performed a de novo functional annotation. BLASTX comparison against the NR database was performed for the 28,027 genes annotated in the genome, with a maximum number of target hits of 20 and a minimum E value of 0. 001. The putative role of immune genes was assessed by investigating related literature. RT-qPCR analyses. Total RNA (1 g) extracted from ground oyster powder was reverse transcribed in 20 l using the Moloney murine leukemia virus reverse transcriptase (MMLV-RT; Invitrogen) according to the manufacturer's instructions using random primers (Invitrogen).
Pipetting and amplification into 384-well plates were performed with a Labcyte acoustic automated liquid handling platform (ECHO) and a Roche LightCycler 480, respectively. The total RT-qPCR volume was 1. 5 l and consisted of 0. 5 l of cDNA (dilution, 1 ⁄4) and 1 l of LightCycler 480 SYBR green Master I mix (Roche) containing 0. 5 M PCR primers (Eurogentec SA; primer list in Table  S2 ). The amplification efficiency of each pair of primers was validated by serial dilution of a pool of all the cDNAs. Relative expression quantification was calculated according to the quantification cycle (2 ϪΔΔCq ) method with the mean of two reference genes (Cg-RPL40, GenBank accession number FP004478; Cg-RPS6, GenBank accession number HS119070). Statistics. Statistical data analysis was carried out in GraphPad Prism (version 6.
ACKNOWLEDGMENTS We thank  M. Leroy  and the  IHPE MIMM team for laboratory assistance  and help during the multiple injection trials. We acknowledge the use of the facilities of the experimental platform of the  Ifremer Palavas Station , the  Ifremer nursery team in Bouin (France) , the  FINA Ifremer program , and the  DECIPHER  project ( ANR-14-CE19-0023 ) for SPF oyster production. C. M. thanks  C. Cosseau  for helpful advice on epigenetic genes. We are grateful to  Philippe Clair  from the qPHD-Montpellier GenomiX platform for useful advice. M. L. was supported by grant from the  Perpignan University (UPVD ) and Ifremer. This work was funded by the  PROVIGAS  project  [Protection antivirale de l'huître Crassostrea  gigas (antiviral protection of the Pacific oyster), la région  Occitanie/ Pyrénées-Méditerranée , project call Chercheuse(r) d'avenir 2013],  Ifremer ,  ANR  project  DECIPHER  ( ANR-14-CE19-0023 ) and  JCJC INVIMORY  ( ANR-13-JSV7-0009 ), and  EU  project  VIVALDI  ( H2020 program , no.
, 2018; Downing et al. , 2011). Moreover, advances in DNA sequencing technology together with the availability of reference genome assemblies for most of the clinically important species  (Downing et al. , 2011; Gonzalez-de la Fuente et al. , 2019; Peacock et al. , 2007; Real et al. , 2013; Rogers et al. , 2011)  in public databases  (Aslett et al. , 2010)  now make it feasible to sequence collections of isolates and determine genetic variants genome-wide. Several studies on the L. donovani complex have applied such an approach including foci in Nepal (16 isolates,  Downing et al. , 2011) , Turkey (12 isolates,  Rogers et al. , 2014) , the Indian subcontinent (204 isolates,  Imamura et al. , 2016) , Ethiopia (41 isolates from 16 patients,  Zackay et al. , 2018)  and  Brazil (20 and 26 isolates, respectively, Teixeira et al.
, 2018 ) and a subset of the extensive dataset available from the Indian subcontinent (N = 33,  Imamura et al. , 2016)  resulting in a total of 151 isolates (Supplementary file 1, visualised at https://microreact. org/project/_FWlYSTGf; Argimon et al. , 2016). Accurate SNP variants were identified across 87. 8% of the reference genome with a genotype quality of at least 10 (median = 99), indicating a < 0. 1 (median = ~10 À10 ) probability of an incorrect genotype call. The remaining 12. 2% could not be assayed as short reads could not be uniquely mapped to repetitive parts of the genome. This identified a total of 395,624 SNP sites out of the 32 Mb reference assembly. We also used these sequence data to infer extensive gene copy-number variation (91. 5% of genes varied in dosage; 7,625/8,330 genes) and larger genome structure variation, including copy numbers of individual chromosomes (aneuploidy) that is common in Leishmania.
Together, these data represent the most comprehensive, global database of genetic variation available for any Leishmania species. Evolution of the L. donovani complex Phylogenetic reconstruction based on whole-genome SNP variation clearly separated L. infantum from L. donovani strains. L. donovani separated into five major groups that coincide with geographic origin (Figure  1A-B , Figure  1 -figure supplement 1) and show a strong signal of isolation-by-distance (IBD) between countries (0. 76, p-value<=0. 0001, Mantel test, Supplementary file 2). While the inferred root of the phylogeny is between L. infantum and L. donovani, groups within L. donovani showed similar levels of divergence as between the two species, with the deepest branches within L. donovani in East Africa. The largest L. donovani group in our collection, Ldon1, included samples from the Indian subcontinent, and could be further divided into two subgroups that separate samples from India, Nepal and Bangladesh from three samples of Sri Lankan origin; both subgroups displayed strikingly little diversity.
All 5 'MON-1' American samples formed a monophyletic sub-clade that was most closely related to parasite strains from Portugal, Spain, Italy and a single isolate from Israel suggesting a South-Western European origin of Central and South American L. infantum (Figure  1 -figure supplement 2A). This result was still valid when including another 26 L. infantum isolates sampled from three states in Brazil  (Carnielli et al. , 2018). They all clustered in a single 'American' clade with little genetic diversity (Figure  1 -figure supplement 2B). L. infantum in Central and Southern America, however, is not generally monophyletic as also one non-MON-1 L. infantum isolated in Panama was present in our dataset (Figure  1 -figure supplement 2). For the zoonotic parasite L. infantum, 12 of our 30 MON-1 strains were isolated from dogs -previously also described as the prevalent zymodeme in dogs  (Pratlong et al.
All the isolates from the Ldon3 group, were also highly heterozygous and so potentially hybrids, but we cannot exclude other possible origins for this heterozygosity (Figures  3A  and 4  , Figure  4 -figure supplement 1, Table  1 ). While the CUK samples are known to be of hybrid origin between a JPCM5-like L. infantum isolate and an unidentified parasite from the L. donovani complex  (Rogers et al. , 2014) , our admixture results did not suggest hybridisation between genetic groups present in our dataset. This still held when varying K (the specified number of subpopulations) from 2 to 25 (Figure  4 -figure supplement 6). We therefore took a haplotype-based approach to increase the power to identify putative parents of these hybrids similar to that in  Rogers et al. (2014) , but now compared them to our larger set of isolates.
Genetic variation for known drug resistance loci We investigated how genetic variation previously associated with drug resistance is distributed across our global collection of isolates, including loci involved in resistance to or treatment failure of antimonial drugs and Miltefosine (Table  2 ). The best-known genetic variant associated with drug resistance in Leishmania is the so-called H-locus: amplification of this locus is involved in resistance to several unrelated drugs including antimonials  (Callahan and Beverley, 1991; Dias et al. , 2007; Grondin et al. , 1993; Leprohon et al. , 2009; Marchini et al. , 2003). In our dataset, the four genes at this locus had an increased gene copy number in 30% of the samples (CN +1 to +44) and a reduced copy number in 9% (CN À1; Table  2 ). 36% of all isolates had a copy number increase of varying degree with identical insertion boundaries that included the genes YIP1, MRPA and argininosuccinate synthase (Figure  9A   2 ).
Four additional isolates showed a copy number increase for only two of the genes at the locus, with different boundaries but always including the MRPA gene (Fig- ure 9-figure supplement 1B). Differential expression of the Mitogen-activated protein kinase 1 (MAPK1) has previously been associated with antimony resistance. However, while  (Singh et al. , 2010)  suggested that overexpression is associated with resistance,  (Ashutosh et al. , 2012)  suggest the opposite effect potentially implicating an impact of the genetic background. As expression in Leishmania is typically tightly linked with gene copy number (Prieto  Barja et al. , 2017; Iantorno et al. , 2017) , we summarised MAPK1 CNVs in our dataset (Table  2 ). 45% of all isolates had an amplified copy number at this locus, including all isolates of Ldon1 and Ldon3 with the highest copy number increase in Ldon1 isolates of between 12 and 41 copies (Figure  9 -figure supplement 2A, Table  2 , Supplementary file 6).
Only a single L. infantum isolate had a reduced copy number of one. Increased copy number of MAPK1 is thus associated with isolates from geographical locations with high historical antimonial drug pressures such as the Indian subcontinent and to a lesser extend Africa. Another protein, the membrane channel protein aquaglyceroporin (AQP1), is known to be involved in the uptake of pentavalent antimonials: reduced copy number and expression have been associated with drug resistance  (Andrade et al. , 2016; Gourbal et al. , 2004; Monte-Neto et al. , 2015; Mukherjee et al. , 2013) , as has other genetic variation at this locus  (Imamura et al. , 2016; Monte-Neto et al. , 2015; Uzcategui et al. , 2008). In our dataset, copy number at this locus was reduced in 6% and increased in 35% of all isolates with small effect sizes (CN À2 to À1 and +1 to +3) but at least one copy of the locus was always present (Figure  9 -figure supplement 2B, Table  2 ).
05, FDR = 1) and 11 and 12 an excess of non-synonymous differences (p-value<0. 05, FDR = 1) for L. donovani and L. infantum, respectively (Figure  9 -figure supplement 3, Supplementary file 11). Interestingly, one of the genes with putative signs of adaptive evolution in L. donovani (LINF_330040400 v41, LinJ. 33. 3220 v38) was previously associated with in vivo enhanced virulence and increased parasite burden in vitro for L. major when overexpressed  (Reiling et al. , 2010). In our dataset, this gene contained nine missense, 3 synonymous and 19 upstream/intergenic SNP-variants differentially fixed between L. donovani and L. infantum (Supplementary file 4), which might provide further candidates for differences in virulence between both species. While genetic variants can become fixed in different populations by either neutral forces (genetic drift) or positive selection, we took advantage of the genetic differentiation between groups to search for group-specific SNPs that might be of biological relevance.
Of regions where the genetic diversity had previously already been described for many samples, we chose subsets representing the known genetic diversity (i. e. Imamura et al. , 2016; Zackay et al. , 2018). In an additional analysis (Figure  1 -figure supplement 2B), we included 26 isolates from three different states in Brazil  (Carnielli et al. , 2018)  to confirm reduced genetic diversity in South America. The 97 samples sequenced for this study are deposited in ENA under the study accession numbers: PRJEB2600 (ERP000767), PRJEB2724 (ERP000966), PRJEB8947 (ERP009989) and PRJEB2115 (ERP000169) (https://www. ebi. ac. uk/ena/data/view). All metadata on the 151 isolates including ENA accession numbers of individual samples are summarized in Supplementary file 1 (see also https://microreact. org/project/_FWlYSTGf; Argimon et al.
Read mapping pipeline Reads were mapped with  SMALT (RRID:SCR_005498, v0. 7. 4, Ponstingl, 2010)  using the parameters: '-x -y 0. 9 -r 1 -i 1500' specifying independence of paired-end reads, a minimum fraction of 0. 9 of matching bases, reporting of a random best alignment if multiple are present and a maximum insert size of 1500 bp against the reference genome JPCM5 of L. infantum (MCAN/ES/98/LLM-877, Tri-TrypDB v38, RRID:SCR_007043;  Aslett et al. , 2010). Mapped reads were sorted and duplicate reads were marked with picard 'MarkDuplicates' (RRID:SCR_006525, v1. 92, https://broadinstitute. github. io/picard/). For resulting individual bam files per isolate, indels were called and local realignment was performed with GATK using the 'RealignerTargetCreator' and 'IndelRealigner' with default settings (RRID:  SCR_001876, v2.
Variant calling Variant calling was done following the Genome Analysis ToolKit (GATK, RRID:SCR_001876) bestpractice guidelines  (Van der Auwera et al. , 2013)  with modifications detailed below. Given the aneuploidy of Leishmania, we considered individual somies per chromosome and isolate: the GATK 'HaplotypeCaller'  (RRID:SCR_001876, v3. 4-0, DePristo et al. , 2011)  was used with the parameters '-sample_ploidy SOMY -dt NONE -annotateNDA' and additionally all-sites files were generated by adding the additional flag '-ERC BP_RESOLUTION' to the above HaplotypeCaller command. Individual vcf files (by chromosome and isolate) were processed, filtered and combined with custom made scripts implementing the following steps: only SNPs outside masked regions (see 'Reference Genome Masking') were extracted; SNPs were hard filtered excluding genotypes failing to pass at least one of the following criteria: DP >= 5*SOMY, DP <= 1.
Phylogenetic reconstruction For phylogenetic reconstruction from whole-genome polymorphism data, all 395,602 SNPs that are polymorphic within the species complex and have a maximum fraction of 0. 2 non-called sites across all 151 samples were considered. Nei's distances were calculated for bi-allelic sites per chromosome with the R package  StAMPP (v1. 5. 1, Pembleton et al. , 2013) , which takes into account aneuploidy across samples. Resulting distances matrices of Nei's distances per chromosome were weighted by chromosomal SNP count forming a consensus distance matrix, that was used for phylogenetic reconstruction with the Neighbor-Joining algorithm implemented in the R package APE (RRID:SCR_ 017343, v5. 2,  Saitou and Nei, 1987). For rooting of the tree, the phylogenetic reconstruction was repeated using three additional outgroup samples, of L.
Then individual g. vcf files were combined using 'GenomicsDBImport', SNPs across all samples were called using 'GenotypeGVCFs' and hard filtered using parameters "QD <2. 0, MQ <50. 0, FS >20. 0, SOR > 2. 5, BaseQRankSum < À3. 1, ClippingRankSum < À3. 1, MQRankSum <À3. 1, ReadPosRankSum <À3. 1 and DP <6'. The resulting vcf file were analysed in R: only SNPs with a missing fraction across samples < 0. 2 were retained; Nei's distances between samples were called using the R package  StAMPP (v1. 5. 1, Pembleton et al. , 2013)  and phylogenetic trees calculated with neighbour joining with the r package ape (RRID:SCR_017343, v5. 3,  Paradis et al. , 2004).
Phylogenetic reconstruction of maxicircles Sequence reads were mapped against the maxicircle DNA of the reference strain, LV9 (MHOM/ET/ 1967/HU3), of L. donovani (TriTrypDB v46, RRID:SCR_007043) with SMALT (RRID:SCR_005498, v0. 7. 4, Ponstingl, 2010) using parameters: '-x -y 0. 8 r À1 -i 1500' and duplicates were marked with picard, 'MarkDuplicates' (RRID:SCR_006525, v1. 92, https://broadinstitute. github. io/picard/). Local indel realignments were performed on the resulting alignments with GATK using the 'RealignerTargetCreator' and 'IndelRealigner' with default settings  (RRID:SCR_001876, v3. 4-0, DePristo et al. , 2011)  and subsequently filtered for a mapping quality of 20 and proper pairs using samtools, parameters '-q 20 f 0 Â 0002 F 0 Â 0004 F 0 Â 0008' (RRID:  SCR_002105, v1. 3, Li et al.
Gene-feature annotation and GO enrichment analysis All SNPs were annotated with gene features using the software SNPeff (RRID:SCR_005191, v4. 2,  Cingolani et al. , 2012). Annotations for the reference genome L. infantum, JPCM5, were downloaded from TriTrypDB (v38, RRID:SCR_007043;  Aslett et al. , 2010). Several gene sets of interest were subsequently tested for Gene ontology (GO, RRID:SCR_002811) term enrichments for the ontology 'biological process'. GO mappings for L. infantum genes were downloaded from TriTrypDB (v38, RRID:SCR_007043), where 4704 of the 8299 annotated coding genes were also associated with a GO term. Enrichment of functional categories was tested using the weightFisher algorithm in topGO (RRID:  SCR_014798, v2. 34. 0, Alexa et al. , 2006)  sing all genes annotated in the 'gene to GO' mapping file (v38).