dataesr/scientific-paragraphs-categorization huggingface dataset:
The PylC, PylD, and PylS datasets gave consistent results (additional files, Figures S1 to S3 available online at http://dx. doi. org/10. 1155/2014/374146) and were therefore concatenated to provide increased phylogenetic signal (Figure  3 ). "Ca. M. alvus", "Ca. M. intestinalis, " and "M. luminyensis" form a robustly supported monophyletic cluster, indicating a common origin of the Pyl system in the 7th methanogenic order. The products of the two pyl gene clusters of M. luminyensis (including the third extra PylC) are more closely related among them than to the other sequences, indicating that they originated from a specific duplication. Methanosarcinaceae and bacterial sequences also form two distinct and robustly supported monophyletic clusters. Moreover, the sequences from the 7th methanogenic order appear to be more closely related to their bacterial counterparts than to Methanosarcinaceae (Figure  3 ).
Acknowledgments This work was funded by PhD Scholarship supports from the  French "Ministère de l'Enseignement Supérieur et de la Recherche " to  Nadia Gaci Simonetta Gribaldo  was supported by the  French Agence Nationale de la Recherche (ANR)  through Grant  ANR-10-BINF-01-01  " Ancestrome. " Paul W. O'Toole  was supported by  Science Foundation Ireland  through a  Principal Investigator award  and by an  FHRI award  to the  ELDERMET  project by the  Department of Agriculture, Fisheries and Marine of the Government of Ireland.
Analysis. Phylogenetic analysis was carried out with each of the Pyl gene products by recovering all homologues available in current sequence databases. It appears that bacterial members belong to related Firmicutes species, and a single deltaproteobacterium, Bilophila wadsworthia, a common resident of the gut. For PylC, PylD, and PylS, no evident closely related homologues Methanolobus psychrophilus 100/1 100/1 0. 2 87/0. 99 Ace toha lobi um arab atic um D Thermincola potens D es ul fo to m ac ul um id an ac et ox Th e r m o a c e t o g e n i u m p h a e u m B il o p h il a w a d s w o r t h ia e D s u lf o t o m a c u lu m g ib s o n ia e s.
frugiperda, we report here the generation of an NGS-based sequencing of RNA libraries obtained from a large variety and number of tissues and developmental time-points. These libraries have been assembled together in a reference transcriptome, dubbed Sf_TR2012b, comprising around 55,000 sequences available for search through a dedicated database. We further used the Sf_TR2012b assembly to annotate several gene families and study their expression profile. We focused especially on the genes involved in immunity and in olfaction. We validated some of our developmental genes predictions by qPCR to demonstrate that in Sf_TR2012b, some regulated transcripts without clear orthology are bona fide functional genes in S. frugiperda.
Construction of a reference transcriptome The reference transcriptome presented here is named Sf_TR2012b. We incorporated various sources of RNA sequences from 454, Illumina and Sanger sequencing and had to develop a custom pipeline. This assembly is described in details in the Methods section and illustrated in Additional file 1: Figure  S1. Its assembly has been performed by the software MIRA  [4]  and incorporated various sources of RNA sequences obtained from a laboratory strain of S. frugiperda. The main source of RNA molecules has been obtained by 454 RNAseq of a library containing cDNA extracted from 27 different samples (Table  1 ) comprising 14 developmental timepoints samples (Table  1 , column A) and 13 dissected tissues (Table  1 , column B), in order to cover the majority of the mRNA produced by the FAW.
Evaluation of the Sf_TR2012b assembly First, we wanted to evaluate whether the last step of assembly, combining the short reads with the long reads sequences, was improving the quality of our transcriptome. We evaluated this improvement by aligning independent Illumina libraries on the different reference transcriptome assemblies using the Bowtie software  [7]. We used 4. fastq files containing sequences for whole larvae RNA extracted in experimental conditions, which were unrelated to our transcriptome project (N. Volkoff, personal communication), to align the reads against Sf_TR2012b and Sf_GATC_clusters. For all 4 samples, the average of unmapped reads is around 20% when aligned on both assemblies (Sf_GATC_Clusters and Sf_TR2012b) (Figure  1A ). However we greatly improved the percentage of uniquely mapped reads from 22% to 45% while decreasing the number of multiple reads from 58% to 32% (Figure  1A ).
This indicates that our second step of assembly was particularly effective at eliminating most of the redundancy that one could expect from an organism that is far from being isogenic. Then, we wanted compare our approach with other insects datasets for which the transcriptome was obtained differently. We thus compared the number of contigs (55,000) and the total size (37 M nt) of the Sf_TR2012b assembly with 2 other insect models, Drosophila melanogaster and Bombyx mori. The Drosophila melanogaster transcriptome contains more than 28,000 sequences for a total length of 81 M nt (Table  2 ). The Drosophila transcripts set is built primarily from computer predictions of the genome sequence and also from permanent curations due to the large Drosophila research community and more than a hundred years of genetics research on this model.
0%) contigs are similar to a Heliconius protein (Hmel1-1_Release_20120601), while 9,709 Heliconius proteins (75. 7%) are similar to a TR2012b contig. Similarly, 21,439 (92. 7%) contigs are similar to a Monarch protein (Dp_geneset_OGS2), while 10,887 Danaus Plexippus (72. 0%) proteins are similar to a TR2012b contig. Both comparisons have been made using blastx with a threshold of 1e-10 (p-value) without complexity filter. We found only 200 contigs (0. 4%) matched prokaryotes (archea, bacteria and viruses), making large scale contaminations of our samples highly improbable. What was more surprising was that 57. 6% of our RNA sequences didn't match any known protein sequences. Among these 31,650 sequences without hit against the NR databank, 491 include at least one protein domain found with the following algorithms (BlastProDom, FPrintScan, Gene3D, HMMPanther, HMMPfam, HMMSmart, PatternScan, ProfileScan and superfamily).
These comparisons have been done using Interproscan (v4. 8) directly on the transcript sequences. This low number of hit does not necessarily reflect that the other sequences are spurious, but they may derived from UTR, or the predicted ORF (Interproscan uses getorf from emboss) might be uncompleted or the sequences may contain a frameshift. This is enhanced by the observation that the transcripts without hit are shorter than the transcripts with hits (mean : 526. 5, median : 446 versus global mean : 868. 7, and global median : 694). Finally, some of the transcripts without any match against NR databank might correspond to lncRNA, as well as transcribed repeat elements. For the sequences matching eukaryota, we performed blastx against eukaryotic core gene sets such as the CEGMA (v2. 5) geneset (http://korflab.
ucdavis. edu/datasets/cegma/) and the BUSCO proteins set (ftp://cegg. unige. ch/OrthoDB7/ BUSCO). We found that 452 among 457 CEGMA proteins (98. 9%) are similar to 1831 Sf_TR2012b contigs. Similarly 2961 among 3369 BUSCO drosophila proteins (87. 9%) and 3025 among the 3299 BUSCO Danaus plexippus  proteins (91. 7%) were present in the Sf_TR2012b assembly, suggesting that the core components of the FAW transcriptome were present in our assembly.
Conversely, we searched for more constrained genes belonging to 3 functional families : homeobox-domain genes (Hox), odorant and pheromone-binding proteins and immune related genes. Hox proteins can be identified with certainty thanks to the conserved signature of their homeodomains even between distantly related species  [8]. Within species however, many paralogs can be identified. Contrary to the ribosomal proteins, the expression of genes encoding Hox proteins in D. melanogaster is usually temporally and spatially restricted. Thus we expected these particular proteins to be more difficult to find in our assembled transcriptome due to their underrepresentation in the RNA samples collected, compared to other abundant transcripts such as ribosomal proteins. We used a collection of 105 D. melanogaster homeodomain protein sequences from the Homeodomain Resource Database  [9]  as a tblastn query of our Sf_TR2012b transcriptome assembly.
Access to Sf_TR2012b through Lepidodb The Sf_TR2012b assembly sequences can be downloaded through the Lepidodb database (www6. inra. fr/lepidodb/Private and http://www. inra. fr/lepidodb/downloads/ TR2012b) [login: lepiduser password: papillon]. In this database, the transcripts can be queried with their identifier or with the name of an ortholog. BLASTx against the nr database have been performed for all transcripts in the assembly and the database contain the best 10 results for each transcript. BLAST searches can also be performed from the database. In addition a GBrowse system has been added that presents a set 44 BAC sequences (25 BAC sequences have been added in LepidoDB to the 19 already published  [31] ) representing each around 100 kb of contiguous genomic sequence, on which a computed gene prediction has been performed by KAIKOGAAS (http://kaikogaas.
dna. affrc. go. jp/). We aligned the Sf_TR2012b transcripts on this genomic reference  [31]. The alignment results by blast have been converted into a. gff3 file directly viewable in a GBrowse genomic viewer instance  [32]  hosted through the Lepidodb portal (http://www6. inra. fr/lepidodb/Private) (Figure  3A ). Throughout the BACs, we have 49. 4% of KAIKOGAAS gene predictions that overlap with an Sf_TR2012b transcript, while 80. 3% of Sf_TR2012b overlap with a gene prediction (Figure  3B ). We also provide for each transcript its level of expression for 10 different samples. These levels of expression are also represented on the BACs present in LepidoDB as tracks (Figure  3A )  [33].
Measure of transcripts expression by RNAseq We were interested in using the reference transcriptome to study families of genes involved in larval development, in chemosensory reception and in immune response. These categories of genes are all necessary to the adaptation of the feeding larva to its environment, whether be its host plant or its panoply of pathogens. Specifically we were interested in investigating the repertoire of genes and their level of expression for each of these categories. Thus, we extracted total RNA corresponding to 6 developmental time-points and 4 experimental dissected tissues to generate a quantitative gene expression dataset by Illumina sequencing (see Methods and Table  1 ). Total RNA has been extracted from the selected samples and sequenced to produce 10 M single-end 50 bp reads (see Additional file 7: Figure  S2 ), representing approximately 1.
Common and specific developmental genes When designing the reference transcriptome, we emphasized the production of 6 developmental time-points specific RNA libraries to be sequenced by Illumina. We were interested in the specificity of development of S. frugiperda compared to other insects. Indeed, S. frugiperda is a pest in its larval stage and resistance to common pesticides has become a particularly prevalent issue  [37]. We reasoned that identifying developmental genes specific to S. frugiperda would provide the community with specific targets for the development of new strategies of pest control. To identify genes involved in the embryonic development of S. frugiperda, we focused on the comparison of 2 Illumina RNAseq libraries sequenced from eggs and L2 stage larvae RNA extracts. We used the R framework package DEseq to identify genes that are overexpressed in eggs and no longer in L2 and that are thus required only during embryonic development.
Conclusions In conclusion, we provide in this study a reference transcriptome for S. frugiperda, available through a dedicated database, along with measures of differential expression across 10 different samples. We found this resource invaluable to annotate and study the expression of different families of genes. In particular, we were able to identify and validate Lepidoptera specific genes involved in development. We also analyzed a set of genes involved in olfaction across two Spodoptera species. And finally, we annotated an almost complete set of immune related genes and observed in particular that some anti-microbial peptides are highly expressed in chemosensory organs, even in absence of induction, raising the possibility that antennae and palps can naturally act as primary organs of immune response, since they are in open contact with the natural environment.
Assembly and alignments The flowchart of the assembly process is presented in Additional file 1: Figure  S1. A first step of assembly of the 454 reads has been performed by GATC Company using the CD-HIT software  [45]  and resulted in 183,373 clusters. The 1,042,944 454 reads clipped from adaptator (GATC) were compared to Univec (https://www. ncbi. nlm. nih. gov/ tools/vecscreen/univec/, version of march 2011) leading to the removal of 3,297 reads similar to known vector sequences. The 79,148 ESTs previously sequenced by Sanger method and coming from 8 different libraries (Sf9 cell lines, Sf21 cell lines, hemocytes, induced hemocytes, midgut, induced midgut, fat body and a tissue mix)  [6]  were also compared to UniVec resulting in the removal of 1132 ESTs. We performed an assembly of the 78,016 ESTs and 1,039,647 454 reads using the MIRA software  [4].
This step resulted in 52,865 contigs with an N50 of 794 bp. Then we clipped the sequences at both ends by 80 bp and mapped the combined 90,454,901 reads from the 10 Illumina libraries onto this reference with Gassst  [46]. The resulting 17,724,510 short reads that were unmapped underwent a subsequent step of short reads assembly using the velvet and Oases softwares  [47]. We finally used MIRA to assemble together the 24,505 contigs from the Velvet/Oases assembly and the 52,865 contigs from the previous MIRA assembly of the 454 sequences and EST sequences. This final assembly is Sf_TR2012b. To produce the time-course expression datasets, the reads from the 10 Illumina libraries have been aligned against Sf_TR2012b and the BACs set using bowtie  [7]. Alignment files were further processed by samtools  [48].
Analysis of expression data The Illumina 50 bp reads were mapped on the transcriptome with bowtie using the options (-a -m 1 -best -strata), reporting the best alignment (i. e. having the least number of mismatches). The raw counts by library were calculated by contig using a home-made program (based on the perl Bio::DB::Sam library). The raw counts have been divided by the total number of aligned reads in order to obtain the RPM normalized values. As well, for each library we computed the normalizations by quantile normalize BetweenArrays function of the limma R library  [32]  and the by size factors (functions estimateSizeFactors and sizeFactors from the DESeq R package  [33] ). BACs 25 BACs from Spodoptera frugiperda genome have been isolated as previously described  [31]. Their sequence has been deposited at the European Bioinformatics Institute nucleotide archive : http://www.
47% to 99. 61% for each sample). Within the 2. 4 Mb critical region on 2p23, we identified 168 variants with 99. 4% coverage of coding exons at $10X coverage. We did not find any copy number variant (CNV) shared among RA cases in this region. We applied stringent filtering criteria to select for rare pathogenic variants present in this family but not in any public database with non-reference allele frequency $0. 01 (dbSNP v132, 1000 Genomes Project Phase I data  [24] , and NHLBI Grand Opportunity Exome Sequencing Project [ESP] 5400  [25] ). Of the 168 protein-coding variants within the critical region on 2p23, only one variant was identified by our filtering approach, a missense mutation in the phospholipase B1 (PLB1) gene (Figure  3B  right panel). For completeness, we also evaluated the other 13 IBD stretches shared in at least one copy among all 5 RA cases and 1 ACPA-positive unaffected subject, and found one additional missense single-nucleotide variant (SNV) and 1 insertion-deletion (Indel) included in the IBD stretches (Table  1 ).
GWAS meta-analysis of the PLB1 locus and RA risk of noncoding variants To provide additional support for the role of PLB1 in risk of RA, we evaluated common (minor allele frequency [MAF]. 5%) and low-frequency variants (MAF 1-5%) near the PLB1 gene locus in an outbred population of European ancestry (Figure  1D ). We used a GWAS of 8,875 seropositive RA cases and 29,367 controls from 11 studies (Table  S4 ). Five of the GWAS datasets were previously unpublished (3,427 cases and 6,837 controls from ReAct, Dutch, anti-TNF response to therapy collection [ACR-REF], COR-RONA, and Vanderbilt)  [29] [30] [31] [32] [33] , which increased the power of our published dataset (Figure  4  and Table  2 )  [2, 4, 7]. We observed no evidence of systematic bias genome-wide (inflation factor l GC = 1. 034; Table  S4 ).
The H3K4me3 mark is particularly informative for cell-type specific overlap with traitassociated variants, and RA risk variants showed significant enrichment in T reg primary cells  [34]. The RA risk SNP of rs116018341 (and the SNPs in absolute LD with it; r 2 = 1. 00), was within one of the H3K4me3 peaks of T reg primary cells (P = 0. 043), while the other risk SNP, rs116541814, was not (Figure  4C ). In a search of public eQTL databases (eQTL Browser and Blood eQTL Browser  [35] , see URL), we found no evidence that either SNP (or SNPs in LD with them, r 2. 0. 8) influenced PLB1 gene expression.
Samples We report a 4-generation, consanguineous pedigree in which 8 of 49 individuals are affected with RA (Figure  2  and Table  S1 ). The pedigree members were recruited by a board-certified rheumatologist from Tishreen Hospital, Damascus, Syria. All   2 ). For the case-control association analysis of RA in the European populations, we studied 8,875 RA cases and 29,367 matched controls for GWAS meta-analysis (Table  S4 ). Six GWAS data has been previously published  [2] ; additional five GWAS datasets are previously unpublished  [29, 32, 33] , as described in more detail below. The 1,088 RA cases and 1,088 matched controls for exon sequencing of PLB1 represent a subset of patients with GWAS data (Table  S4 ). All the subjects were confirmed to be of European origins using both self-reported ethnicities and the results of principal component analysis (PCA).
We assessed significance of the probability to observe the IBD stretch shared among 5 ACPA-positive RA cases and 1 ACPA-positive unaffected subject while only 6 of 16 ACPAnegative unaffected subjects by a permutation procedure (610,000 iterations). For each of the iteration steps, we randomly selected 6 members from the pedigree as ''affected'' subjects, and assessed whether at least one of the IBD stretches observed among these 6 ''affected'' subjects were observed in #6/18 of the other ''unaffected'' subjects. Java TM software for this novel IBD mapping method is available at http://plaza. umin. ac. jp/ ,yokada/datasource/software. htm. We did not apply parametric linkage analysis methods for SNP genotype data such as Merlin  [49] , since the software did not work properly due to the complex pedigree structure including multiple loops.
Calling of CNV was conducted by using the ExomeDepth software version 0. 9. 7  [54]. Whole-exome sequencing data of the pedigree is available to other researchers upon request. Filtering of the identified variants was conducted according to the following processes: (i) variants likely to be pathogenic (missense, nonsense, frameshift Indels, or splice-site acceptor/ donor); (ii) variants not registered in the databases (dbSNP v132, 1000 Genomes Project Phase I data  [24] , and ESP5400  [25] ) with non-reference allele frequency $0. 01; (iii) for Indels, ones not located 65 bp of known variants; and (iv) variants of which $1 non-reference alleles were observed in all the exome sequenced 4 RA cases.
For the subset of reads which mapped to the pathogen genome, further PCR duplicate detection was performed. PCR-duplicated reads are known to arise during the final amplification steps of library preparation, and are related to the reduced complexity of captured libraries. In our case, 148,695,448 duplicated reads were found in the dataset mapping to A. phagocytophilum, which corresponded to a percentage of 66. 6%. This value is relatively high as compared to those obtained by other authors in similar contexts. For instance,  Bright et al. (2012)  reported to have obtained percentages of duplicated reads ranging from 12. 1% to 46. 2% in Plasmodium vivax field samples, with an average of 26. 1%  [26]. The properly paired A. phagocytophilum reads retained after the removal of duplicates  (74, 464, 184)  were combined with the properly paired unmapped reads (24,532,400, as previously stated) to obtain a final dataset of 98,996,584 reads, roughly corresponding to 9.
Sequencing, genome assembly and annotation The sequencing of the captured library was performed on an entire lane of a flow cell on a HiSeq2000 (Illumina) sequencer as paired-end 108-bp reads. Base calls were performed using RTA software. After removing low quality bases (PHRED < 10) from 3' read extremities, the reads were mapped to the genomes of B. taurus (UMD3. 1) and A. phagocytophilum HZ (NC_007797. 1) using bwa (v 0. 6. 1) [73]  with default parameters. After this step, the reads mapping to the B. taurus genome where filtered out. Therefore, only the properly paired reads mapping to the A. phagocytophilum genome and the properly paired reads which where unmapped on both genomes, were retained. PCR duplicates were detected on the subset of reads which were mapped to the pathogen genome and then subtracted using samtools rmdup  [74].
On the other hand, the filtering of duplicated reads was not carried out on unmapped reads. The A. phagocytophilum properly paired reads retained after the removal of duplicates were combined to the properly paired unmapped reads to obtain a final dataset. Before the assembly step, the reads were digitally normalized using Khmer  [27, 28]  with a kmer size of 20, and a cutoff of 20. 17% of paired-reads were conserved. Remaining reads were then assembled using Velvet 1. 2. 07  [75] , with a kmer value between 51 and 91, and default parameters. Genome annotation was performed using Agmial  [76]. Functional annotation was performed by preferentially transferring annotations from the A. phagocytophilum HZ strain where appropriate. The 169 scaffolds have been deposited in the European Nucleotide Archive (CCXQ01000001-CCXQ01000169).
Genome alignment The 169 assembled scaffolds were aligned against A. phagocytophilum strain HZ (NC_007797. 1) using MUMmer 3. 0 software  [77]. All perfectly matching regions were determined using the mummer module with default parameters, and were used to determine coverage statistics. Then the Nucmer module was used with default parameters to cluster and extend exact matches in order to create the largest similarity regions. These large regions were then used to generate dot-plots.
Core genome determination The BOV 10_179 genome sequenced in this study (CCXQ01000001-CCXQ01000169), and the genomes of the previously sequenced Human strains HZ (NC_007797. 1), HZ2 (NC_021879. 1), HGE1 (APHH0 1000001. 1 and APHH01000002. 1), the dog strain Dog 2 (NC_021881. 1), the horse strain MRK (PRJNA216999), the rodent (Zapus hudsonius) strain JM (NC_021880. 1), the sheep Norway Variant 2 (PRJNA217033) and the Ixodes scapularis, collected from white-tailed deer, strains CRT35 (PRJNA217037) and CRT38 (APHI01000001. 2 and APHI01000002. 1) were compared for core genome determination. Complete protein sequences encoded by non-pseudogenes in the ten different genomes were compared using protein BLAST to define groups of orthologous proteins by single linkage clustering (e-value <10-3; >70 % identity over >98 % of the longest sequence length).