dataesr/scientific-paragraphs-categorization huggingface dataset:
Transcriptome Taxonomic Assigment Taxonomic assigment of the whole collection of the de novoassembled transcripts was performed with two different tools: (i) with MEGAN version 6. 20. 18 (Huson et al. , 2016)  after a DIAMOND-blastx (version 2. 0. 6) interrogation of the mRNA transcripts with parameters: -sensitive -max-target-seqs 100 -e value 1e-3 -F 15 -range-culling against the RefSeq protein database, accessed in February 2021  (Buchfink et al. , 2021) ; and (ii) with Blobtool2 implemented in BlobToolKit version 1. 3. 6  (Challis et al. , 2020)  after a BLASTn (version 2. 9. 0) interrogation of the mRNA transcripts against the nt database accessed in February 2021 with the following parameters: -max_hsps 1 -max_target_seqs 10 -evalue 1e-25 and a DIAMOND-blastx interrogation of the mRNA transcripts against the RefSeq protein database with parameters: -sensitive -max-target-seqs 1 -evalue 1e-25 (version 2.
ONT Bioinformatics Analyses Raw fast5 Nanopore reads were base-called with Guppy (version v4. 5. 3), then base-called fastq reads were searched for residual Nanopore adapter and PCR primer with PoreChop (version 0. 2. 4) (https://github. com/rrwick/Porechop) with the following parameters: -adapter_threshold 80 -middle_threshold 80 -extra_middle _trim_bad_side 10. Clean fastq reads were then aligned with minimap2 (version: 2. 18) to the full-length B. ostreae rRNA operon (accession: MZ305451) assembled and identified in the present study. Finally, the alignment file was finally processed with SAMtools to produce a sorted and indexed bam file, which was subjected to variant calling with the tools PEPPER (https://github. com/ kishwarshafin/pepper) (version 0. 4) with parameters -ont -gvcf.
Assembly of rRNA Reads We decided to run the de novo assembly of the rRNA reads with the Transabyss assembler because it has been shown to be more accurate at assembling non-coding RNA reads  (Lau et al. , 2018). We eventually assembled 3189 rRNA contigs that we have clustered and manually curated in Geneious to finally obtain 468 nonredundant rRNA transcripts (Table  1 ). We considered those 468 rRNA transcripts as highly representative of the rRNA population of our transcriptome since 96. 72% of our identified rRNA reads were aligned back to our de novo-assembled rRNA transcripts (Table  1 ). To accurately determine the taxonomic classification of these rRNA transcripts, we then predicted and extracted both SSU and LSU from all rRNA transcripts and submitted each collection of SSU or LSU to SILVA ATC for taxonomic classification.
Since the SILVA SSU database is more complete and accurate than the LSU database, we have focused our SILVA taxonomic classification results on the SSU transcripts. This SILVA classification, along with BlastN top hit identification, showed that more than 90% of the SSU rRNA reads aligned to SSU rRNA transcripts belonging to the order Haplosporida and more precisely to B. ostreae (Figure  3A ). We then wanted to obtain a complete overview of the B. ostreae diversity in our parasite-enriched sample. To this end, we first identified the putative full-length B. ostreae rRNA operon (see Materials and Methods) and used it as a reference for read mapping and variant calling (see Materials and Methods). Our B. ostreae rRNA operon diversity analysis shows that the majority of the rRNA operon sequences are either identical or show minor diversity in our samples.
A total of 10,100 bases called nanopore reads, corresponding to 15. 4 million bases, were sequenced. Among them, 1,573 reads with an average size of 4,048 bp, a GC% of 52. 4%, and composed of 6 million sequenced bases were aligned to the B. ostreae operon, leading to a mean coverage of 880 reads. In contrast, the 8,067 unaligned reads showed an average size of 1,089 bp, a GC% of 35%, and a total of 9. 4 million sequenced bases. Those unaligned reads were mostly assigned to the non-rRNA region of O. edulis (using the draft genome of O. edulis data not shown). This result can be explained by the fact that a total of 200 ng of old O. edulis DNA (12 years stored at -20°C) was used for the four LR-PCR reactions used to amplify the B. ostreae rRNA operon. As nanopore sequencing is known to promote the sequencing of small DNA molecules, this could have resulted in a large number of small O.
Assembly of mRNA Reads Assembly of the mRNA reads was done with the assembly pipeline DRAP-Trinity (see Materials and Methods). This pipeline allows us to increase the compactness of the assembly, reduce the chimerism of the de novon-assembled contigs, and improve the error rates of the whole transcriptome  (Cabau et al. , 2017). Finally, we were able to de novo assemble 15,618 highquality transcripts with an average size of 1,425. 5 bp and a GC% of 37. 1% (Table  1 ). As observed for the rRNA assembly, we are confident that our mRNA assembly is highly representative of our transcriptome since 96. 97% of our identified mRNA reads were aligned back to our de novo-assembled mRNA transcripts (Table  1 ). Furthermore, BUSCO analysis of our de novoassembled transcriptome identified the group Alveolata as the most complete BUSCO.
Since the group Rhizaria to which B. ostreae belongs does not have its own BUSCO database and that Alveolata is the closest group of the Rhizaria inside the SAR clade, we believe that our transcriptome is mainly represented by transcripts closely related to protozoa belonging to SAR clade and very likely to B. ostreae (Figure  4 ). To obtain a global overview of the taxonomic range of our transcriptomic dataset, we then performed a taxonomic classification of the 15,618 identified transcripts with two different tools: MEGAN and BlobToolsKit (see Materials and Methods). The MEGAN taxonomic classification identified that more than 90% of the whole cumulative expression of our transcriptome measured as transcripts per million (TPM) is represented by 11,682 transcripts distributed as follows: 4,338 transcripts with no hits in the RefSeq protein database, 5,572 transcripts belonging to Eukaryota, 1,376 transcripts belonging to a cellular organism, 198 transcripts that could not be assigned to any taxonomic group by MEGAN (usually when a transcript have several Blast hits with the same evalue) and 198 transcripts assigned to Bilateria (Figures  5A, B ).
The BlobToolKit taxonomic classification showed a similar trend with most of the cumulative expression (>60%) represented by transcripts with no hit in the interrogated databases (Figure  5C ). Interestingly, this blobplot showed that the transcripts with no hit tend to have a lower GC% than the transcripts harboring a hit with a taxonomy assignment (Figure  5C ). These results support the hypothesis that our B. ostreae transcriptome is mainly represented by transcripts without homology to other sequences present in the interrogated database, and this emphasizes the need to obtain more reference sequences for B. ostreae and other protozoa from the Rhizarian clade. polycistronic transcripts. Additionally, 2,018 transcripts do not contain any ORFs (Figure  6A ), these orphan transcripts likely represent other populations of noncoding RNA different from the rRNAs (i.
Annotation of mRNA Transcripts The 24,995 identified proteins have then been extensively annotated through several database interrogations (see Materials and Methods). Among them, 6,673 protein-coding transcripts were mapped to GO terms from InterProScan similarity hits, with the most represented terms belonging to "catalytic activity," "binding," "metabolic process," and "cellular process" (Figure  6B ). Furthermore, 9,950 proteins were assigned to putative functional categories using the Clusters of Orthologous Groups of proteins (COG) database from the EggNOG mapper (Figure  6C ). The most representative COG function includes "A-RNA processing," "J-Translation," and "O-Post-translational modification," while the less representative COG function includes "W-Extracellular structure," "Y-Nuclear structure," and "N-Cell motility" (Figure  6C ).
This functional annotation gives some insight into the metabolic activity of B. ostreae and confirms the active function related to its intracellular lifestyle. The observed functional annotation was further validated by a fine analysis of the top 30 expressed transcripts, which represent more than 20% of the cumulative expression measured as TPM (Table  2 ). Indeed, among those 30 transcripts, 6 transcripts do not harbor any ORF and do not have similarities with the Rfam database, 6 transcripts possess a predicted ORF but without identified functional annotation, and 12 transcripts code for a protein with functions involved in transcription and translation, which reflects an important biological activity (Table  2 ).
ostreae. However, the presence of transcripts from other organisms reflects the difficulty of obtaining a pure B. ostreae sample. Thus, we decided to process this transcriptome as a meta-transcriptome and divide the analysis into two main parts: (i) the evaluation of the rRNA population and (ii) the evaluation of mRNA. Indeed, even if the molecular biology procedure used here was not designed to sequence rRNA molecules since the cDNAs were produced by polyT-based library preparation, we were able to identify 2. 3% of rRNA reads in our dataset. The evaluation of this rRNA population allowed us to confirm that the great majority of rRNAs present in our dataset corresponds to B. ostreae [>80% of the total rRNA expression (TPM)]. Furthermore, the extensive functional analysis performed on the putative protein-coding-containing transcripts assembled from mRNA reads showed that the majority of this metatranscriptome codes for proteins with no homology with proteins belonging to the current major protein database.
Moreover, most of the cumulative expression observed in this meta-transcriptome corresponds to transcripts with no hits in those databases or to transcripts without precise taxonomic assignment. This observation emphasizes the lack of -omic data regarding the Rhizarian clade and, more specifically, regarding B. ostreae in a public biological database. This transcriptome is a highly valuable data source to start characterizing the complex biology of this intracellular parasite. The transcriptome has some limitations associated with difficulties in obtaining pure parasite, but the most complete dataset available is likely to prove extremely useful. Additionally, those data will help to develop a new diagnostic technique to improve the identification of B. ostreae in complex mollusc samples. Thanks to this dataset, we were able to assemble the putative full-length B.
ACKNOWLEDGMENTS The authors would like to thank  Jean-Yves Daniel  and  Arnaud Huvet  ( Ifremer, PFOM-LPI, Brest, France ) for their valuable advice on the development of RNA extraction and control. The authors are also very grateful to  Christophe Klopp  ( INRAE, Plate-forme Genotoul Bioinfo, Auzeville, France ) for the preliminary read analyses and his advice for this RNASeq approach. Illumina sequencing was performed at the  Platforme GeT-PlaGe, Genotoul Centre INRA, Toulouse, France.
DATA AVAILABILITY STATEMENT The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www. ncbi. nlm. nih. gov/, PRJNA731671 https://www. ncbi. nlm. nih. gov/ genbank/, MZ305451 https://figshare. com/, https://doi. org/10. 6084/m9. figshare. 16592480. v1 https://figshare. com/, https://doi. org/10. 6084/m9. figshare. 16598633. v1.
CFIGURE 5 | 5 FIGURE 4 | BUSCO transcriptome evaluation. (A) Simplified Stramenopila, Alveolata, and Rhizaria (SAR) phylogeny adapted from Sierra et al. (2016). The Rhizaria clade does not have a proper BUSCO database due to the lack of genomic data in this clade, whereas the closely related clades Alveolata and Stramenopila are represented by proper BUSCO databases. (B) BUSCO evaluation of the B. ostreae transcriptome. This evaluation has been performed by interrogating a large panel of different database clades related or not to B. ostreae and its host, O. edulis.
We hypothesize that the strong structure of O. neglectus populations across a variable environment 5. Our study supports the expectation that genome scans may be prone to errors in highly structured populations. Nonetheless, it also empirically demonstrates that careful statistical controls enable the identification of candidate loci that invite more detailed investigations. Ultimately, genome scans are valuable tools to help guide further research aiming to determine the potential of non-model species to adapt to future environments. + + + + (c) (d) (e) whether population connectivity driven by kelp rafting counteracts local adaptation. To address these aims, we first undertook demographic analyses to explore how variation in the different datasets is structured. Second, we compared geographical patterns of variation from 'outlier' sets of markers (derived from significant results of two methods of genome scans for outliers [GSO]: PCAdapt and Bayescan) with those from the remaining loci, to identify potentially different evolutionary histories for these sets of loci.
| Sampling, DNA extraction and genotypingby-sequencing A total of 188 samples of O. neglectus from populations distributed across New Zealand and its subantarctic islands were available from previous work  (Nikula et al. , 2012; Salloum et al. , 2020). In all, 16 populations were included in the NZ-wide dataset (spanning 17 latitudinal degrees, and 167 individuals after filtering, Figure  1a , Table  1 , Table  S1 ). A subset of 10 of these populations was included in the southern dataset, representing the previously identified South clade  (Salloum et al. , 2020). One population (Christchurch) was excluded from the southern dataset due to excess missing genotype data (see below; final southern dataset spanning six latitudinal degrees, and 88 individuals, Table  1 , Table  S1 ). The upper North Island (NI) dataset has 45 individuals (two latitudinal degrees), and the South Island (SI) dataset has 43 individuals (one latitudinal degree).
| Demographic patterns Population differentiation and genomic diversity (observed and expected heterozygosity, overall F ST , and Nei's F IS ) were estimated for all datasets, with the r package hierfstat v. 0. 5-7  (Goudet & Jombart, 2020). The same package was used to estimate overall observed heterozygosity within each population in each dataset, and to check heterozygosity distribution across loci within populations. Principal component analyses were done with the r package lea v. 2. 6. 0  (Frichot & François, 2015)  for all datasets (Table  S2 ). Allele frequencies were compared for all datasets with a custom R script. In addition, Weir and Cockerham's F ST was estimated for all loci with VCFtools  (Danecek et al. , 2011)  for pairwise comparisons between clades in the NZ-wide dataset (North-South, North-Central and Central-South).
The allele frequency spectrum was derived for each population NZ-wide, using vcf2sfs  (Liu, 2020). In addition, δaδi  (Gutenkunst et al. , 2009)  was used to plot 2D allele frequency spectra between North-Central, North-South and Central-South clades. To further analyse population structure in the NZ-wide dataset, a population assignment plot was done with the r package Adegenet v. 2. 1. 1  (Jombart, 2008; Jombart & Ahmed, 2011) , as this enables the optimal number of ancestral populations and the admixture proportions of these different ancestries within each individual to be inferred  (Jombart & Ahmed, 2011). Finally, directional relative migration was estimated for the NZ-wide dataset with the DivMigrate function in the DiveRsity r package  (Keenan et al. , 2013) , using the G ST method  (Sundqvist et al.
| Genome scans for outliers To identify 'outlier' markers, two GSO were run on each dataset separately, the PCA-based method PCAdapt  (Luu et al. , 2017)  and the Bayesian method Bayescan  (Foll & Gaggiotti, 2008). For PCAdapt, principal component analyses were initially undertaken for 50 principal components (20 for the NI and SI datasets), with the default 5% threshold for minor allele frequencies and no LD thinning (LD was controlled by the random choice of one SNP per rad-tag when calling SNPs). The default Mahalanobis distance was used as a test statistic, and p-values were corrected using the genomic inflation factor. The proportion of variance explained by the principal components was checked (Figure  S4 ) and the analysis was re-run with three principal components for the NZ-wide dataset, two principal components for the NI dataset, one for the SI dataset, and six for the southern dataset, using the same settings as above.
After checking their distribution (Figures  S5-S11 ), the corrected p-values were transformed into q-values with the qvalue package v. 2. 16 in r  (Storey et al. , 2019). For Bayescan, runs were performed with the default configurations (5,000 iterations, thinning interval = 10, 20 pilot runs of 5,000 length, burn-in = 50,000, prior odds for the neutral model = 10). The results of both methods were combined by calculating geometric means of the q-values of PCAdapt and Bayescan for each SNP, that is, qval combined = 2 √ qval A × qval B , where qval A and qval B correspond to the locus q-value resulting from the two different methods. The geometric mean can be seen as a form of arithmetic mean on the logarithmic scale and, compared to the absolute scale, is therefore more closely linked to the natural interpretation of the q-value.
| Genotype-environment association The PCA-based method RDA and the Bayesian method BayeScEnv were used to search for loci in association with environmental variables (GEA-loci). Environmental variables (described below) were mean-centred and scaled to a variance of one. RDA implemented in the r package vegan v. 2. 5. 6  (Oksanen et al. , 2019)  does not allow missing data; therefore, missing loci were imputed (imputing with the most common genotype within the clade). In addition, the environmental variables used as predictors in the RDA should not be correlated  (Dormann et al. , 2013) ; thus, only variables with a squared correlation coefficient smaller than or equal to 0. 7 with each other were included, leaving six semi-independent variables in the NZwide and southern datasets, and three semi-independent variables in the NI and SI datasets (Tables  S4  and S5 ).
The RDA performed was set to scale SNP loadings to unit variance. Analyses of variance (ANOVA) with 999 permutations were undertaken to test the RDA full model for significance. To obtain q-values, RDA loadings treated as z-scores were mean-centred and scaled to unit variance, which were transformed into two-tailed p-values that were converted to q-values using the qvalue r package  (Storey et al. , 2019) , following checks of their distributions (Figures  S13-S16 ). For BayeScEnv  (de Villemereuil & Gaggiotti, 2015) , correlation between environmental variables is not of concern because they are tested independently in multiple runs of a univariate method; thus, all environmental variables were included. In addition, no imputation is required for BayeScEnv analyses. Runs were performed with the default configurations (5,000 iterations, thinning interval = 10, 20 pilot runs of 5,000 length, burn-in = 50,000, upper bound for the Uniform prior of parameter g = 10, mean alpha prior = -1, prior probability for non-neutral models = 0.
1, prior preference for the locus-specific model = 0. 5) in the NZ-wide, NI, SI and southern datasets. The environmental variables available to test in the GEA analyses were latitude, longitude and long-term means of monthly sea-level pressure  (Kanamitsu et al. , 2002) , net shortwave, precipitation rate, surface pressure, total cloud cover, zonal wind velocity, meridional wind velocity, air temperature  (Kalnay et al. , 1996) , ocean temperature, ocean salinity  (Monterey & Levitus, 1997) , salinity, sea surface height relative to geoid  (Behringer & Xue, 2004 ) and sea surface temperature  (Ishii et al. , 2005). With the exception of latitude and longitude, environmental variables were downloaded from NOAA (data provided by the NOAA ESRL Physical Sciences Division, Boulder, Colorado, USA, from their website at http://www.
esrl. noaa. gov/psd/) as netcdf files, and converted to a data frame with the ncdf4 v. 1. 16. 1 r package  (Pierce, 2019). The geographical coordinates of the populations sampled and of the points where environmental data are available are not exactly matched in NOAA's satellite data; thus, the value of the closest point to the actual sample coordinates was used. In addition, RDA and BayeScEnv were also run using a set of 100 randomly generated variables, normally distributed with a mean of zero and standard deviation of one. Our rationale behind these runs was to identify loci that generated spurious genotype-environment association, suggesting that significant results with the true environmental data for these loci are more likely to be spurious association. Both methods were run using the same parameters as for the real environmental variables.
The results of the two GEA methods were combined using the same method as for GSO: calculating geometric means of the qvalues of RDA and BayeScEnv for each locus and environmental variable pair, in the real and random datasets. Loci identified as significantly associated with real variables (combined q-value <0. 05) were then examined for association with the randomly generated variables. Those loci that also significantly associated with fewer than five random variables (i. e. 5% of the random runs) were treated as high-confidence associations. Since BayeScEnv runs included the environmental variables correlated with latitude (r 2 > 0. 7), q-values of all latitude-correlated variables were combined using the geometric mean, and significant associations to latitude in the 'broad sense' are reported.
High-confidence GEA loci were checked for observed heterozygosity and missing data across populations. Sequences of high-confidence GEA loci were submitted in a BLAST search against a custom database containing the Acanthopleura granulata genome and transcriptome  (Varney et al. , 2021) , using the blastn algorithm with default configurations  (Altschul et al. , 1990). To identify genes and functions, the longer A. granulata transcriptome sequences matching our short O. neglectus queries (92-base fragments) were extracted and run in another BLAST search against other molluscs (taxid:6447) using the nucleotide (nr/nt) database, and the program selection blastn, optimizing for somewhat dissimilar sequences.
| SNP datasets The SNP calling pipeline resulted in a NZ-wide dataset consisting of all 16 populations and 10,987 loci, with population-average 0. 2% missing data (Figure  S3A , Table  1 , Table  S2 ). SNPs were called with a population map grouping the three clades because it resulted in a more complete dataset with less missing data than grouping by population only (average 11% missing data; Figure  S3B ). The NI dataset consists of 12,012 loci with 0. 0% missing data (Figure  S3C , Table  1 , Table  S2 ). The SI dataset consists of 7,476 loci with population-average 0. 1% missing data (Figure  S3D , Table  1 , Table  S2 ). The southern dataset consists of 9 populations and 13,004 loci, with population-average 8. 7% missing data (Figure  S3E , Table  1 , Table  S2 ). One population (Christchurch; CR) was removed from the SI and southern datasets due to excess missing data; note that it had less missing data in the NZ-wide dataset and hence is included there.