dataesr/scientific-paragraphs-categorization huggingface dataset:
Simulation of Deletion Frequencies The method was adapted from  Zhang and Gerstein (2003). Briefly, we created a random nucleotide string with the same 65 length as the total length of pseudogenes in our dataset and simulated deletions events over different generations. Every m k generations, we introduced a deletion of length k (with k ranging from 1 up to 3) into the string. By fixing m 1 and varying m 2 and m 3 , we can estimate which ratio of deletion 70 frequencies best fits with our observed data. Simulations were stopped when the number of deletions of size 1 reached the observed number. From every combination of 27 m 2 33 and 27 m 3 33 we found that m 2 ¼ 29 and m 3 ¼ 29 yielded a deletion distribution closest to our observed data (results for every combination were averaged over 100 runs).
We counted 396 deletions, for a total size of 7,421 nt and 300 insertions which summed up to 1,340 nt. This results in an insertion to deletion ratio of 0. 758, corre-15 sponding to a size ratio of 0. 181. Indels are thus strongly biased towards deletions in the genomes of leaf nodule symbionts, even taking into account that our data likely underestimates the total number of deletions since the requirement for the framing of pseudogenes by functional 20 genes excluded large deletions spanning multiple genes. For comparison, we also extracted a set of 24 pseudogenes from various free-living Burkholderia species using the same methods and criteria. From these, we identified 42 deletions (1,662 bp in total) and 33 insertions (398 bp). 25 The insertion to deletion rate in free living Burkholderia is 0.
786, for a size ratio of 0. 239, comparable to what we observed in leaf nodule symbiont genomes. These values are also consistent with the experimental values 16 insertions for 17 deletions (ratio ¼ 0. 941, size ratio ¼ 0. 321) 30 reported for a mutation accumulation experiment in B. cenocepacia by Dillon et al. (Dillon et al. 2015). Size of Small Deletions, but Not Insertions, Show a Distinct Bias towards Multiples of 3 We next exploited our high quality dataset to ask whether the 35 indel mutation rate is affected by local %G þ C. We measured indel rate as the frequency of indels per kb. We estimated the ancestral length of pseudogenes by calculating the length of the alignment of the two functional orthologs. We did not observe any significant correlation between the rate of indels and the %G þ C of the pseudogene (Pearson correlation coefficient ¼ 0.
0265, P-value ¼ 0. 74). During our analyses, however, we noticed that many of the pseudogenes in our dataset seemed to accumulate indels with sizes multiple of 3 nt (heretofore referred to as 3n indels). Because pseudogenes are not assumed to be translated into functional proteins, selective constraints resulting in the preservation of the reading frame are not expected. To investigate whether 3n indels are indeed overrepresented in our dataset, we analyzed further the distribution of indel sizes in our data. The size distribution of deletions or insertions follows a power law, as established previously for prokaryotic and eukaryotic genomes  (Cartwright 2009; Sung et al. 2016). The tail of the distribution is very long, with few insertions or deletions exceeding a size of 15 nt. We therefore focused on indels of sizes 15 nt in order to simplify the statistical analyses.
The results of 10 our simulations show that in order to reach the observed deletion distribution, single deletions of size 2 and 3 should each occur 0. 34 (10/29) times as often as deletions of size 1. These ratios are very similar to the observed ratios of deletions in our dataset (44/131 or 0. 336 and 45/131 or 15 0. 344 for sizes 1/2 and 1/3, respectively), indicating that almost all of the observed deletions of sizes 2 and 3 are derived from single deletion events, rather than combinations of smaller deletions. Size Bias of Small Deletions Is Not Dependent on the Age 20 of the Pseudogene Because insertion sizes do not show a systematic bias towards sizes 3n, and deleterious mutations are presumed to accumulate neutrally in pseudogenes, we hypothesized that the observed bias towards 3n deletions is not due to selective 25 constraints to preserve the reading frame.
We did not observe any correlation between the calculated age of pseudogenes and the ratio of 3n deletions (fig. 3 ). Similarly, no significant correlation between pseudogene age and ratio of 3n deletions was observed for the datasets of Ca. B. umbellata and Ca. B. brachyanthoides (supplementary fig. S3 , Supplementary Material online). As further confirmation, the total number of indels, which can be seen as a proxy for the age of a pseudogene, did not correlate with the proportion of 3n deletions on our entire dataset (supplementary fig. S4 , Supplementary Material online). Together with the fact that insertion sizes do not show a significant bias towards multiples of 3 nt, our data show that selective constraints for the preservation of the reading frame.
GBE early in the evolutionary history of pseudogenes are not responsible for the observed bias towards 3n deletions. Sequence Motifs Associated with Small Deletions Hint at Molecular Mechanism 5 The sequence context surrounding indels may provide clues about the mechanisms causing small indels, in particular deletions. Short deletions generally arise from slipped strand mispairing (SSM) or error-prone repair of double strand breaks (DSB)  (Garcia-Diaz and Kunkel 2006). SSM usually occurs in 10 homopolymer runs, or in short direct repeats and results in short deletions or insertions  (Levinson and Gutman 1987). To investigate whether SSM is a significant contributor to the formation of small deletions in our pseudogene dataset, we analyzed indels for the presence of homopolymeric runs and 15 tandem repeats (TR).
GBE (Proteobacteria, Firmicutes, and Actinobacteria). The identified pseudogenes contained 104 deletions (including 81 with size 15 nt) and 57 insertions (including 52 15 nt). Linear regression analysis of log-transformed count and size data did 5 not reveal a significant contribution of size 3n over the whole range of size values (P-value ¼ 0. 119), possibly because of a lack of statistical power and because deletions larger than 8 nt are sparse in this dataset. Bootstrap analysis however showed a significant overrepresentation of deletions of 3 and 6 nt 10 (fig. 5A ). Insertions of size 3 nt are no more abundant than expected, however insertions of 6 nt seem to be overrepresented (fig. 5B ). Data for insertions >8 nt are too sparse to infer more general trends. Insertions of 6 nt seem to be evenly distributed across taxa and we could not observe any partic-15 ular motif explaining their overrepresentation.
However, and because other insertions of size 3 nt are not overrepresented, this is unlikely to be due to selective constraints for preservation of gene reading frame. Together, these data suggest that background mutations are biased towards 3n deletions in a 20 wide array of bacteria. We could not rule out confidently that the ratio of 3n deletions does not correlate with pseudogene age because individual data sets were too small for robust statistical analysis. In order to establish what mechanisms caused the short deletions in this dataset, we searched the 25 deletion sites for motifs. Similar to what we observed for Ca. Burkholderia pseudogenes, slightly over 39% of the deletion sites (41 out of 104) contained microhomologies (1-11 nt). Genome-Wide Microhomology Interval Lengths Also Show a Bias towards Multiples of Three 30 Microhomologies flanking deletion sites resemble the signature of alternative end joining (A-EJ), a form of NHEJ recently characterized in E.
Discussion Although the neutrally evolving properties of pseudogenes make them valuable to the investigation of intrinsic mutation profiles and biases, their short evolutionary retention time usually limits their usefulness. Here, we identified a large set of pseudogenes having accumulated multiple indel mutations which have unusually retained closely related functional homologs in free-living bacteria. Two factors probably contribute to the higher retention rate of pseudogenes and the large size of this dataset: 1) leaf nodule symbiosis is polyphyletic and has evolved several times in Burkholderia (now Caballeronia) over the past 10-15 Ma  (Lemaire et al. 2011; Pinto-Carb o et al. 2016) , resulting in multiple free-living and host captive lineages evolving in parallel; 2) vertical transmission and repeated transmission bottlenecks result in relaxed genome-wide purifying selection on the genomes of the Ca.
This beetle was first identified as a pest in 1859 in the Midwestern United States, after it expanded from its native host plant, Solanum rostratum (Solanaceae), onto potato (S. tuberosum)  3. As testimony to the difficulty in controlling L. decemlineata, the species has the dubious honor of starting the pesticide industry, when Paris Green (copper (II) acetoarsenite) was first applied to control it in the United States in 1864  4. Leptinotarsa decemlineata is now widely recognized for its ability to rapidly evolve resistance to insecticides, as well as a wide range of abiotic and biotic stresses  5  , and for its global expansion across 16 million km  2  to cover the entire Northern Hemisphere within the 20th century  6. Over the course of 150 years of research, L. decemlineata has been the subject in more than 9,700 publications (according to the Web of Science TM Core Collection of databases) ranging from molecular to organismal biology from the fields of agriculture, entomology, molecular biology, ecology, and evolution.
We compared the size of gene families associated with particular traits against existing available genomes from related species, particularly those sequenced by the i5k project (http://i5k. github. io), an initiative to sequence 5,000 species of Arthropods. While efforts have been made to understand the genetic basis of phenotypes in L. decemlineata (for example, pesticide resistance)  32, 43, 44  , previous work has been limited to candidate gene approaches rather than comparative genomics. Genomic data can not only illuminate the genetic architecture of a number of phenotypic traits that enable L. decemlineata to continue to be an agricultural pest, but can also be used to identify new gene targets for control measures. For example, recent efforts have been made to develop RNAi-based pesticides targeting critical metabolic pathways in L.
Results and Discussion Genome Assembly, Annotation and Assessment. A single female L. decemlineata from Long Island, NY, USA, a population known to be resistant to a wide range of insecticides  47, 48  , was sequenced to a depth of ~140x coverage and assembled with ALLPATHS  49  followed by assembly improvement with ATLAS (https://www. hgsc. bcm. edu/software/). The average coleopteran genome size is 760 Mb (ranging from 160-5,020 Mb 50 ), while most of the beetle genome assemblies have been smaller (mean assembly size 286 Mb, range 160-710 Mb)  [51] [52] [53] [54] [55]. The draft genome assembly of L. decemlineata is 1. 17 Gb and consists of 24,393 scaffolds, with a N50 of 414 kb and a contig N50 of 4. 9 kb. This assembly is more than twice the estimated genome size of 460 Mb  56  , with the presence of gaps comprising 492 Mb, or 42%, of the assembly.
We assessed the completeness of both the ALLPATHS and REDUNDANS assemblies, and the OGS separately, using benchmarking sets of universal single-copy orthologs (BUSCOs) based on 35 holometabolous insect genomes  59  , as well as manually assessing the completeness and co-localization of the homeodomain transcription factor gene clusters in the ALLPATHS assembly. Using the reference set of 2,442 BUSCOs, the ALLPATHS genome assembly, REDUNDANS genome assembly, and OGS were 93. 0%, 91. 9%, and 71. 8% complete, respectively. We found an additional 4. 1%, 5. 4%, and 17. 9% of the BUSCOs present but fragmented in each dataset, respectively. For the highly conserved Hox and Iroquois Complex (Iro-C) clusters, we located and annotated complete gene models in the ALLPATHS genome assembly for all 12 expected orthologs, but these were split across six different scaffolds (Supplementary Figure  1S  and Table  1S ).
0%, was assigned to Platyhelminthes). Gene Annotation, Gene Family Evolution and Differential Expression. We estimated a phylogeny among six coleopteran genomes (A. glabripennis, Agrilus planipennis, D. ponderosae, L. decemlineata, Onthophagus taurus, and T. castaneum; unpublished genomes available at http://i5k. github. io/)  51, 52, 54 using a conserved set of single copy orthologs and compared the official gene set of each species to understand how gene families evolved the branch representing Chysomelidae. Leptinotarsa decemlineata and A. glabripennis (Cerambycidae) are sister taxa (Fig. 1 ), as expected for members of the same superfamily Chrysomeloidea. We found 166 rapidly evolving gene families along the L. decemlineata lineage (1. 4% of 11,598), 142 of which are rapid expansions and the remaining 24 rapid contractions (Table  1 ).
Population Genetic Variation and Invasion History. To understand the propensity for L. decemlineata pest populations to rapidly evolve across a range of environmental conditions, we examined geographical patterns of genomic variability and the evolutionary history of L. decemlineata (to see the genomic analysis of environmental adaptation, refer to the Diapause and Environmental Stress section of Supplementary File). High levels of standing variation are one mechanism for rapid evolutionary change  77. We identified 1. 34 million biallelic single nucleotide polymorphisms (SNPs) from pooled RNAseq datasets, or roughly 1 variable site for every 22 base pairs of coding DNA. This rate of polymorphism is exceptionally high when compared to vertebrates (e. g. ~1 per kb in humans, or ~1 per 500 bp in chickens)  78, 79  , and is 8-fold higher than other beetles (1 in 168 for D.
Raw RNAseq reads for 8 different populations were obtained from previous experiments: two Wisconsin populations (PRJNA297027)  130  , a Michigan population (PRJNA400685), a lab strain originating from a New Jersey field population (PRJNA275431), and three samples from European populations (PRJNA79581 and PRJNA236637)  43, 121. All RNAseq data came from pooled populations or were combined into a population sample from individual reads. In addition, RNA samples of an adult male and female from the same New Jersey population were sequenced separately using Illumina HiSeq 2000 as 100 bp paired end reads (deposited in the GenBank/EMBL/DDBJ database, PRJNA275662), and three samples from the mid-gut of 4 th -instar larvae were sequenced using SOLiD 5500 Genetic Analyzer as 50 bp single end reads (PRJNA400633).
Genome Sequencing, Assembly, Annotation and Assessment. Four Illumina sequencing libraries were prepared, with insert sizes of 180 bp, 500 bp, 3 kb, and 8 kb, and sequenced with 100 bp paired-end reads on the Illumina HiSeq 2000 platform at estimated 40x coverage, except for the 8kb library, which was sequenced at estimated 20x coverage. ALLPATHS-LG v35218  49  was used to assemble reads. Two approaches were used to scaffold contigs and close gaps in the genome assembly. The reference genome used in downstream analyses was generated with ATLAS-LINK v1. 0 and ATLAS GAP-FILL v2. 2 (https://www. hgsc. bcm. edu/software/). In the second approach, REDUNDANS was used  57  , as it is optimized to deal with heterozygous samples. The raw sequence data and L. decemlineata genome have been deposited in the GenBank/EMBL/DDBJ database (Bioproject accession PRJNA171749, Genome assembly GCA_000500325.
1, Sequence Read Archive accessions: SRX396450-SRX396453). This data can be visualized, along with gene models and supporting data, at the i5k Workspace@ NAL: https://i5k. nal. usda. gov/Leptinotarsa_decemlineata and https://apollo. nal. usda. gov/lepdec/jbrowse/  147. We estimated the size of our genome using a kmer distribution plot in JELLYFISH  148  , where we mapped the 100 bp paired-end reads from the 180 bp insert library, used a 19 bp kmer distribution plot, and corrected for ploidy. Automated gene prediction and annotation were performed using MAKER v2. 0 149 , using RNAseq evidence and arthropod protein databases. The MAKER predictions formed the basis of the first official gene set (OGS v0. 5. 3). To improve the structural and functional annotation of genes, these gene predictions were manually and collaboratively edited using the interactive curation software Apollo  150.
For a given gene family, known insect genes were obtained from model species, especially T. castaneum  51  and D. melanogaster  151  , and the nucleotide or amino acid sequences were used in BLASTx or tBLASTn  152, 153  to search the L. decemlineata OGS v0. 5. 3 or genome assembly, respectively, on the i5k Workspace@NAL. All available evidence (AUGUSTUS, SNAP, RNA data, etc. ), including additional RNAseq data not used in the MAKER predictions, were used to inspect and modify gene predictions. Changes were tracked to ensure quality control. Gene models were inspected for quality, incorrect splits and merges, internal stop codons, and gff3 formatting errors, and finally merged with the MAKER-predicted gene set to produce the official gene set (OGS v1. 0; merging scripts are available upon request).
For focal gene families (e. g. peptidase genes, odorant and gustatory receptors, RNAi genes, etc. ), details on how genes were identified and assigned names based on functional predictions or evolutionary relationship to known reference genes are provided in the Supplementary Material (Supplementary File). To assess the quality of our genome assemblies, we used BUSCO v2. 0 59 to determine the completeness of each genome assembly and the official gene set (OGS v1. 0), separately. We benchmarked our data against 35 insect species in the Endopterygota obd9 database, which consists of 2,442 single-copy orthologs (BUSCOs). Secondly, we annotated and examined the genomic architecture of the Hox and Iroquois Complex gene clusters. For this, tBLASTn searches were performed against the genome using orthologous Hox gene protein sequences from T.
castaneum (Tcas3. 0) and A. glabripennis. Provisional L. decemlineata models were refined, and potential gene duplications were identified, via iterative and reciprocal BLAST and by manual inspection and correction of protein alignments generated with ClustalW2 154 , using RNAseq expression evidence when available. Finally, we used BlobTools v1. 0  155  to assess the assembled genome for possible contamination by generating a Blobplot. This plot integrates guanine and cytosine (GC) content of sequences, read coverage and sequence similarity via blast searches to assess genome contamination. Putative contaminants were identified using the combined results ('bestsumorder' command) of BLASTn against NCBI's nt database  156  and DIAMOND BLAST  157  against the UniProt reference database  158  , and read coverage was assessed by mapping the 100 bp paired-end reads from the 180 bp insert library to the ALLPATHS genome with the Burrows-Wheeler Aligner (bwa-mem) version 0.
Gene Family Evolution. In order to identify rapidly evolving gene families along the L. decemlineata lineage, we obtained ~38,000 ortho-groups from 72 Arthropod species as part of the i5k pilot project (Thomas et al. in prep) from OrthoDB version 8  160. For each ortho-group, we took only those genes present in the order Coleoptera, which was represented by the following six species: A. glabripennis, A. planipennis, D. ponderosae, L. decemlineata, O. taurus, and T. castaneum (http://i5k. github. io/)  51, 52, 54. Finally, in order to make accurate inferences of ancestral states, families that were present in only one of the six species were removed. This resulted in a final count of 11,598 gene families that, among these six species, form the comparative framework that allowed us to examine rapidly evolving gene families in the L.
We used RAxML  162  with the PROTGAMMAJTTF model to make gene trees from the alignments and ASTRAL  163  to make the species tree. ASTRAL does not give branch lengths on its trees, a necessity for gene family analysis, so the species tree was again given to RAxML along with a concatenated alignment of all one-to-one orthologs for branch length estimation. Finally, to generate an ultrametric species tree with branch lengths in millions of years (my) we used the software r8s  164  , with a calibration range based on age estimates of a crown Coleopteran fossil at 208. 5-411 my  165. This calibration point itself was estimated in a similar fashion in a larger phylogenetic analysis of all 72 Arthropod species (Thomas et al. in prep). With the gene family data and ultrametric phylogeny as the input data, gene gain and loss rates (λ) were estimated with CAFE v3.
Gene Expression Analysis. RNAseq analyses were conducted to (1) establish male, female, and larva-enriched gene sets and (2) identify specific genes that are enriched within the digestive tract (mid-gut) compared to entire larva. RNAseq datasets were trimmed with CLC Genomics v. 9 (Qiagen) and quality was assessed with FastQC (http://www. bioinformatics. babraham. ac. uk/projects/fastqc). Each dataset was mapped to the predicted gene set (OGS v1. 0) using CLC Genomics. Reads were mapped with >90% similarity over 60% of length, with two mismatches allowed. The number of reads were corrected to reads per million mapped to allow comparison between RNAseq datasets that having varying coverage. Transcripts per million (TPM) was used as a proxy for gene expression and fold changes were determined as the TPM in one sample relative to the TPM of another dataset  167.
The Baggerly's test (t-type test statistic) followed by Bonferroni correction was used to identify genes with significant enrichment in a specific sample  168. Statistical values for Bonferroni correction were reported as the number of genes x α value. This stringent statistical analysis was used as only a single replicate was available for each treatment. Enriched genes were removed, and mapping and expression analyses were repeated to ensure low expressed genes were not missed. Genes were identified by BLASTx searching against the NCBI non-redundant protein databases for arthropods with an expectation value (E-value) <0. 001. Transposable Elements. We investigated the identities (family membership), diversity, and genomic distribution of active transposable elements within L. decemlineata in order to understand their contribution to genome structure and to determine their potential positional effect on genes of interest (focusing on a TE neighborhood size of 1 kb).
To identify TEs and analyze their distribution within the genome, we developed three repeat databases using: (1) RepeatMasker  169  , which uses the library repeats within Repbase (http://www. girinst. org/repbase/), (2) the program RepeatModeler  170  , which identifies de-novo repeat elements, and (3) literature searches to identify beetle transposons that were not found within Repbase. The three databases were used within RepeatMasker to determine the overall TE content in the genome. To eliminate false positives and examine the genome neighborhood surrounding active TEs, all TE candidate models were translated in 6 frames and scanned for protein domains from the Pfam and CDD database (using the software transeq from Emboss, hmmer3 and rps-blast). The protein domain annotations were manually curated in order to remove: (a) clear false positives, (b) old highly degraded copies of TEs without identifiable coding potential, and (c) the correct annotation when improper labels were given.