dataesr/scientific-paragraphs-categorization huggingface dataset:
g. , it was observed the T. brucei alternative oxidase that possesses an internal MTS in addition to the classical N-terminal and cleavable one (  38 )). It is also noteworthy that MTSs do not display any strict conservation of the length, sequence, or amino acid composition, rendering their identification a nontrivial task. We used the web application MITOPROT, which calculates the N-terminal protein region that can support an MTS and a cleavage site  (39) , as a predictive tool to identify putative mt isoforms (Table  1 ). The analysis of our data from nanopore RNA-Seq revealed alternative trans-spliced mRNAs for four L. tarentolae aaRSs (Fig. 6  and Table  1 ). For both the AsnRS and the GlnRS, two transcripts were detected: a long version, encoding a predicted N-terminal MTS, corresponds to the sequence proposed in the TriTryp database, and a short version lacking the MTS concluded to encode the cytosolic isoform.
In these two cases, only one read was detected for the mt isoforms, while nearly 60 reads were counted for the cytosolic ones. Two alternatively trans-spliced transcripts were also detected for the HisRS. However, in this case, the significantly larger number of reads obtained for the longer transcript (61 versus 1 for the shorter mRNA) suggests that it corresponds to the cytosolic isoform. It contains a potential internal MTS sequence, which would be concealed in the long form and exposed in the shorter transcript. Localization experiments must be carried out to confirm this hypothesis. Finally, it is the case of IleRS, for which an alternative transsplicing event has been demonstrated in T. brucei  (20). Only two reads were counted in our data for this aaRS transcript. They are both shorter than the sequence provided by the TriTryp database and thus likely encode the cytosolic isoform, yet with a low coverage.
Sequence alignments The aaRS core MSA were previously created  (25). These alignments (one per aminoacylation system) are constituted of aaRS sequences from 92 organisms, which are representative of bacterial, archaeal, and eukaryotic phylogenetic diversity. For the present study, we introduced aaRS protein sequences from Leishmania and Trypanosoma species (the list of the organisms is provided in Fig. S1 ). Newly introduced sequences are individually aligned against the core alignment using the add option of the MAFFT program  (55). The alignments are subsequently clustered according to sequence identity using the DPC (Density Point Clustering) program  (56) , and the number of clusters is determined using the Bayesian Information Criterion  (57). Unrooted phylogenetic trees based on all the aligned sequences are built using the FastME program  (58).
Retro-transcription, sample preparation, and cDNA sequencing were performed according to nanopore Direct cDNA Sequencing instructions (SQK-DCS109). The obtained cDNA library was loaded and sequenced on MinION R9. 4. 1 flow cells (Oxford Nanopore Technologies) for 48 h. Base calling was performed using Guppy basecaller ONT (Oxford Nanopore Technologies). Reads were aligned to the reference genome (TriTrypDB-9. 0_LtarentolaeParrotTarII_Genome, obtained from the Tri-Tryp database (26)) using minimap2 (60) and visualized using IGV tool  (61). In house-made python scripts were developed to select only reads aligning with aaRS transcripts and containing the spliced leader sequence. The script has been deposited in the GitHub database https://github. com/ FlorianBrnrd/aars-nanopore-pipeline.
Three-dimensional modeling The extracted L. tarentolae aaRSs sequences from the Tri-Tryp database were used to generate homology-based models using the online web service SWISS-MODEL  (62). The used templates were selected based on their homology scores to the query sequences and the resolutions of their structures. The produced homology models only serve the purpose of illustrating the positions and sizes of the additional trypanosomatid-specific domains and insertions/extensions. For group 2 aaRSs, templates 5xix (63) and 4ye8  (64)  were used to model the C-terminal domain (CTD) and NTD of AsnRS, 6swx (  51 ) and 4kf9 to model the CTD and NTD of MetRS, 4hvc (  65 ) and 1vki to model the CTD and NTD of ProRS, and the TyrRS X-ray structure (3p0j) was already solved from L. major  (33)  and used as a template for L.
Abbreviations: C, cytosolic; M mitochondrial; nd, not detected. Sequences for Leishmania tarentolae aminoacyl-tRNA synthetases (aaRSs) obtained from our nanopore RNA-Seq experiments are compared with those proposed in the TriTrypDB (https://tritrypdb. org/tritrypdb/ (26)), for which the ID is listed in the left column. TriTrypDB IDs in gray correspond to sequences that are either misannotated or incomplete in the database. Corrected sequences are provided in Supporting information. a MTS probability is predicted on the translated protein sequences by the MITOPROT software (39). b Both alternatively spliced mRNAs were detected in our RNA-seq data with the exception of the IleRS for which only the shorter cytosolic form was sequenced. c The same mechanism as described by Rettig et al. (20) or Nilsson et al.
Acknowledgments-We deeply acknowledge  Derrick Robinson  ( MFP, Bordeaux ) for meaningful inputs on the manuscript,  Delphine Dargere  ( ARNA, Bordeaux ) for help preparing the nanopore library, and  Cameron Mackereth  ( ARNA, Bordeaux ) for linguistic improvements of the manuscript. Molecular graphics and analyses performed with  UCSF ChimeraX. This work was supported by  "Institut National de la Santé et de la Recherche Médicale" (INSERM) ,  "Centre National de la Recherche Scientifique" (CNRS) ,  Université de Bordeaux , and the  French National Program "Investissement d'Avenir"  ( Labex MitoCross ), administered by the  "Agence National de la Recherche " [ ANR-11-LABX-0057_MITOCROSS ]. Author contributions-C. P. , L. M. , Y. H. , D. D. , and M. S. formal analysis; C. P. , L. M. , Y. H.
Consensus genetic map The genetic map was built with Lep-MAP3 (Rastas 2017), as described in https://doi. org/10. 15454/QEDX2V. The resulting map had 3961 fully-informative markers (abxcd segregation) without missing data (missing marker genotypes being automatically imputed in Lep-MAP3). These data were numerically recoded in biallelic doses (0,1,2) according to the initial biallelic segregation and phase (Supplementary Table  S1 ).
Ridge regression: Ridge regression (RR,  Hoerl and Kennard 1970)  adds to the OLS a penalty on the effects using the L2 norm and solves the following equation: b RR ¼ argmin k jjY À Xbjj 2 2 þ kjjbjj 2 2. As a result, all estimated allelic effects are shrunk toward zero, yet none is exactly zero. The amount of shrinkage is controlled by a regularization parameter (k). We tuned it by cross-validation using the cv. glmnet function of the R/glmnet package version 3. 0-2  (Friedman et al. 2010)  with default parameters, except family ¼"gaussian" and a ¼ 0, keeping the k value that minimizes the mean square error (MSE). Note that effects associated with correlated predictors are averaged so that they are close to identical, for a high level of regularization.
Extreme gradient boosting: We first applied LASSO for dimension reduction and then Extreme Gradient Boosting, a popular machine learning method  (Mason et al. 1999) , to estimate marker effects. Hence, we called this method LASSO. GB. As gradient boosting is a nonlinear method, it can take into account any nonlinear interaction between markers, providing better prediction. Briefly, Extreme Gradient Boosting iteratively updates the estimation of weak predictors, in order to reduce the loss. This method requires an optimization of many parameters associated with the loss function (MSE). This optimization was done with train function from R/ caret package version 6. 0-86  (Kuhn 2008)  using the "xgbTree" method. As the optimization of numerous parameters was computationally heavy, we fixed some of them (nrounds ¼ max_depth ¼ 2, colsample_bytree ¼ 0.
Marginal false discovery rate: Marginal false discovery rate (mFDR) has been defined by  Breheny (2019)  as a modified version of the FDR in which those variables correlated with the causative features are not considered as false discovery. This study provided an accurate estimation of mFDR for a given k when using EN or LASSO, thus allowing the selection of a more conservative value of k in order to remain below a given mFDR threshold. We applied mFDR with the R/nvcreg package version 3. 12. 0 (Breheny 2019). For QTL detection on experimental data, we set mFDR to 10% for LASSO. mFDR and EN. mFDR. To our knowledge, this approach had not been adapted yet to a multivariate framework.
For penalized regression methods, parameters were optimized with specific functions such as cv. glmnet and cv. spring. • As performance metrics, we used mainly Pearson's correlation (corP) but we also calculated the root mean square predicted error (RMSPE), Spearman correlation (corS), the model efficiency  (Mayer and Butler 1993) , and test statistics on bias and slope from the linear regression of observations on predictions  (Piñeiro et al. 2008). For experimental data, the whole nested cross-validation process was repeated 10 times (r ¼ 10), whereas for simulated data it was performed only once, but on 10 different simulation replicates (r ¼ 1 and t ¼ 10). The 14 traits were analyzed jointly for MTV_RR, MTV_LASSO, and MTV_EN. But for SPRING, since analyzing all traits together was computationally too heavy, we split traits into three groups by hierarchical clustering (Supplementary Figure  S6 ) performed with R/hclust applied to C.
Simulated data: The quality of a predictor selection method is usually assessed through the relationship between statistical power (i. e. , the true positive rate, TPR) and type I error rate (i. e. the false positive rate, FPR). To compare methods, we thus used the ROC (receiver operating characteristic) curve  (Swets et al. 1979) , which is the plot of TPR as a function of FPR over a range of parameters (Table  2 ), and the pAUC (partial Area Under the Curve; McClish 1989; Dodd and Pepe 2003). Any marker selected at 62 cM of a simulated QTL was counted as a True Positive. For methods with two tuning parameters, one parameter was kept constant (a at 0. 7 for EN and EN. mFDR, and k 2 at 10e-8 for SPRING). We tested several a values for EN but it did not change much the results (data not shown).
arabidopsis. org/, last accessed on 07-21-21) databases based on homologies to access a complete description of gene function, name, and corresponding locus in Arabidopsis. In addition, we used the GREAT (GRape Expression Atlas) RNA-seq data analysis workflow (https://great. colmar. inrae. fr/app/GREAT, last accessed on 07-21-21), which gathers published expression data, to assess the level of expression of our candidate genes in grapevine leaves and shoots, relevant organs for the traits under study. RNA-seq data are normalized as detailed in the 'User manual' section of the GREAT platform: "from the raw read counts, the normalized counts (library size normalization) and the RPKM (gene size normalization) are calculated for each gene in each sample. " Data were retrieved with all filters set to "Select All" except for the organ considered that was restricted to 'Leaves' and 'Shoot'.
Traits with different heritability values: To further compare prediction accuracy between univariate and multivariate methods, we simulated two correlated traits with different heritability values, 0. 1 and 0. 5. MTV_LASSO performed slightly better than univariate LASSO for the lowest heritability trait; however, differences were not significant (Supplementary Figure  S11 ). On the opposite, prediction accuracy was lower with MTV_LASSO than with univariate LASSO for the highest heritability trait, reaching quite low values with 200 simulated QTLs. The same trends were also visible for MTV_EN and EN. MTV_RR never improved prediction over RR and spring. reg never performed better than RR. Since these results were unexpected, we also compared prediction accuracy of the above methods using the simulated data published by  Jia and Jannink (2012).
Candidate genes After merging the QTLs colocalized between traits, we obtained 12 intervals, located on chromosomes 1, 4, 10, 12, 13, 17, and 18, harboring a total of 3461 genes according to the VCost. v3 annotation  (Canaguier et al. 2017 ). Among them, 2379 had an NCBI Refseq identifier and 1757 a putative function (Supplementary Table  S20 ). We then focused our analysis on the eight "new" intervals, i. e. , those that were not overlapping with QTL intervals repeated over years by  Coupel-Ledru et al. (2014 , 2016). They encompassed 1155 genes, half of which were annotated. We were able to retrieve from TAIR and/or UniProt a more precise description of the gene function for 86% of the annotated genes (Supplementary Table  S20 ). The remaining ones either did not have any homologous gene in Arabidopsis or were not described in the above-mentioned databases.
Multivariate vs univariate: When the same heritability was simulated for both trait variables, no superiority of multivariate methods was observed, even when both traits had QTLs at the same positions (Figure  1  and Supplementary Figure  S8 ). When different heritability values were simulated for the two traits, we observed a slight superiority of MTV_LASSO (resp. MTV_EN) over LASSO (resp. EN) only in the "same" and "major" configuration (with both traits sharing the same two QTLs) for the trait with small heritability (Supplementary Figure  S11 ). Other authors who tested multivariate GP on simulated data systematically applied different heritability values; they found a superiority of multivariate methods over univariate ones for the trait with the smallest heritability  (Calus and Veerkamp 2011; Guo et al.
Code availability All software we used was free and open-source and most analyses were done with R (R Core Team 2020), notably graphs that were created using the ggplot2 package  (Wickham 2016 ), version 3. 3. 2. All R scripts used for analysis, i. e. genetic mapping, simulation, phenotypic analysis, prediction, and QTL detection, are available online at https://doi. org/10. 15454/NOUQY2. Many of the custom functions we used are available in a package for reproducibility purposes, R/rutilstimflutre  (Flutre 2019).
| Vgsc non-synonymous nucleotide variation To identify variants with a potentially functional role in pyrethroid resistance, we extracted single nucleotide polymorphisms (SNPs) that alter the amino acid sequence of the VGSC protein from the Ag1000G phase 2 data resource (The Anopheles gambiae 1000 Genomes Consortium, 2020). We then computed their allele frequencies among 16 mosquito populations defined by species and country of origin. Alleles that confer resistance are expected to increase in frequency under selective pressure; therefore, we filtered the list of potentially functional variant alleles to retain only those at or above 5% frequency in one or more populations (Table  1 ). The resulting list comprises 23 variant alleles, including the known L995F, L995S and N1570Y resistance alleles, and a further 20 alleles which prior to Ag1000G had not previously been described in anopheline mosquitoes.
We reported 12 of these novel alleles in our overall analysis of the 765 samples in the Ag1000G phase 1 data resource 1 Codon numbering is given here relative to transcript AGAP004707-RD as defined in the AgamP4. 12 gene-set annotations. A mapping of codon numbers from AGAP004707-RD to Musca domestica, the system in which kdr mutations were first described  (Williamson et al. , 1996) , is given in Table  1. Note: AO, Angola; GH, Ghana; BF, Burkina Faso; CI, Côte d'Ivoire; GN, Guinea; GW, Guinea-Bissau; GM, Gambia; CM, Cameroon; GA, Gabon; UG, Uganda; GQ, Bioko; FR, Mayotte; KE, Kenya; Ac, An. coluzzii; Ag, An. gambiae. Species status of specimens from Guinea-Bissau, Gambia and Kenya is uncertain (The Anopheles gambiae 1000 Genomes Consortium, 2020). All variants are at 5% frequency or above in one or more of the 16 Ag1000G phase 2 populations, with the exception of 2,400,071 G > T which is only found in the CMAg population at 0.
| Data We used variant calls and phased haplotype data from the Ag1000G Phase 2 AR1 data release (https://www. malar iagen. net/data/ag100 0g-phase -2-ar1). Variant calls from Ag1000G Phase 2 are also available from the European Nucleotide Archive (ENA; http://www. ebi. ac. uk/ena) under study PRJEB36277. were aligned to the An. gambiae AgamP3 reference genome assembly  (Holt et al. , 2002). Aligned bam files underwent improvement, before variants were called using GATK UnifiedGenotyper. Quality control included removal of samples with mean coverage <= 14x and filtering of variants with attributes that were correlated with Mendelian error in genetic crosses.
| Haplotype networks Haplotype networks were constructed using the median-joining algorithm  (Bandelt et al. , 1999)  as implemented in a Python module available from https://github. com/malar iagen/ ag100 0g-phase 2-vgsc-report. Haplotypes carrying either L995F or L995S mutations were analysed with a maximum edge distance of two SNPs. Networks were rendered with the Graphviz library and a composite figure constructed using Inkscape. Non-synonymous edges were highlighted using the SnpEff annotations  (Cingolani et al. , 2012).
| Positive selection Core haplotypes were defined on a 6,078-bp region spanning Vgsc codon 995, from chromosome arm 2L position 2,420,443 and ending at position 2,426,521. This region was chosen as it was the smallest region sufficient to differentiate between the ten genetic backgrounds carrying either of the known resistance alleles L995F or L995S. Extended haplotype homozygosity (EHH) was computed for all core haplotypes as described in  Sabeti et al. (2002)  using scikit-allel version 1. 1. 9 (Miles & Harding, 2016) , excluding non-synonymous and singleton SNPs. Analyses of haplotype homozygosity in moving windows (Figures  S1  and S2 ) and pairwise haplotype sharing (Figure  S3 ) were performed using custom Python code available from https://github. com/malar iagen/ ag100 0g-phase 2-vgsc-report.
3 ± 9. 1 million reads per library, corresponding to 1,088,318 unique reads. On average, about 80% of sequences could be mapped unambiguously to the cattle reference genome, while 20% were outmapped (more than five genomic locations, e. g. , mapped to repetitive sequences) or unmapped (Additional file 3: Table  S1 ). As described in "Methods" section, reads were annotated using several databases and classified as micro-RNAs (miRNAs), Piwi-interacting RNAs (piRNAs or piRNA-like), tsRNAs, rsRNAs, circular RNAs (circR-NAs), mRNAs, and other RNAs. Sperm sncRNA comprises a wide variety of sncRNA classes, mainly miRNAs, piRNAs, tsRNAs, and rsRNAs (Additional file 3: Tables  S2-S5 ). Reads annotated as mRNAs accounted for only 0. 8% (ranging from 3000 to 17,000 reads, according to the sampling region), highlighting the absence of somatic cell RNA contaminations.
The primary cattle sncRNA sperm content is dramatically remodeled as sperm mature along the epididymis The sperm sncRNAome showed great plasticity and very large differences in terms of sperm sncRNA content were observed, making it possible to discriminate between sperm samples collected at various maturation stages along the male tract. One of the most striking findings of our study is the observed enrichment of piRNAs within testis parenchyma relative to the sperm fraction isolated from more distal (corpus and cauda) regions of the epididymis as well as ejaculated sperm. Conversely, an increase in miRNA and tRNA content was observed from testis parenchyma to ejaculated sperm, while rRNAs showed a particular expression profile, peaking at epididymis corpus. Overall, our data are consistent with previous studies in mouse and recapitulate prior findings, both in terms of sncRNA content and global expression trends  [12, 19, [30] [31] [32] [33].
However, data unavailability as well as species and/or methodological differences precludes a detailed comparison with these previous findings. For instance, Sharma et al. [30]  identified tRFs as the primary small RNA population in epididymal sperm but read mappings to rRNAs were discarded before normalization. Here, rRNAs were kept in the analysis and were identified as the major sncRNA population in epididymis corpus. We also report an increase in tRFs expression in epididymis. In addition, among 1. 6 million sncRNA sequences described by Chu et al. in mice, only 56,316 perfectly match our bovine sequences (71% rRFs, 15% piRNA, 9% tRFs, and 5% miRNAs). The remaining sequences may reflect both species specificity and the plethora of isoforms produced through RNA editing mechanisms. In this respect, concordance between mouse and bovine increases when considering canonical miRNAs instead of isomiRs represented by individual reads.
Likewise, four profiles were associated with increased expression at epididymis corpus, two at epididymis cauda, and one in ejaculated sperm. Altogether, our data suggest at least three major waves of drastic changes in sperm sncRNA content, first between testis and epididymis corpus, then at epididymis cauda and later in ejaculated sperm. Indeed, the transition from testis parenchyma to epididymis corpus is characterized by a shift from piRNAs to other sncRNA classes, following a drastic decrease in expression of sncRNAs belonging to clusters 1-4, together with a slight increase in miRNA and tsRNA expression (clusters 2-11) and a transient burst in rsRNA expression at COR (cluster 6-9). The second wave of change between epididymis corpus and cauda is characterized by a decrease in rsRNA expression associated with an increase in miRNA and tsRNA expression.
Such an event cannot thus be detected in our NGS study, due to the between-sample normalization procedure. Alternatively, if cytoplasmic droplets are enriched in specific sncRNAs, ejection will result in relative proportion changes that will be observed in our dataset. Since sncRNAs downregulated from PAR to SPZ were clustered according to four different expression profiles of variable intensity of expression decrease, the second option seems more likely. In addition, while expression of sncRNA belonging to clusters 2 and 4 decreases from epididymis caput to corpus, a greater drop in expression was observed between parenchyma and epididymis caput for piRNA belonging to cluster 1 and 3. This suggests that cytoplasmic droplet ejection alone cannot account for the observed changes and other mechanisms may also be at work to account for these piRNAs belonging to clusters 1 and 3.
Interestingly, 86% of small piRNA (≤ 24nt) from Cluster 5 were found to align perfectly on larger piRNA (≥ 25nt) from Clusters 1-4, which gather most piRNAs produced at testis parenchyma. This finding suggests that these small piRNAs may feature degradation products or isopiRs produced by editing during the transit from parenchyma to epididymis caput. Alternatively, they may only reflect false positives present in primary databases  [45]. It has also been proposed that piRNA gained at epididymis may be produced in situ by sperm following a particular unknown biogenesis pathway  [44]. Given that epididymal sperm nuclear gene expression is completely repressed, this latter hypothesis seems unlikely. However, piRNAs may be produced by the mitochondrial genome  [46]  and may be derived from rsRNAs or tsRNAs  [47].