{'fraction_non_alphanumeric': 0.024725943970767356, 'fraction_numerical': 0.03215590742996346, 'mean_word_length': 5.342603321745848, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Predicting genome sizes and restriction enzyme recognition-sequence probabilities across the eukaryotic tree of life

High-throughput sequencing of reduced representation libraries obtained through digestion with restriction enzymes – generically known as restriction-site associated DNA sequencing (RAD-seq) – is a common strategy to generate genome-wide genotypic and sequence data from eukaryotes. A critical design element of any RAD-seq study is a knowledge of the approximate number of genetic markers that can be obtained for a taxon using different restriction enzymes, as this number determines the scope of a project, and ultimately defines its success. This number can only be directly determined if a reference genome sequence is available, or it can be estimated if the genome size and restriction recognition sequence probabilities are known. However, both scenarios are uncommon for non-model species. Here, we performed systematic in silico surveys of recognition sequences, for diverse and commonly used type II restriction enzymes across the eukaryotic tree of life. Our observations reveal that recognition-sequence frequencies for a given restriction enzyme are strikingly variable among broad eukaryotic taxonomic groups, being largely determined by phylogenetic relatedness. We demonstrate that genome sizes can be predicted from cleavage frequency data obtained with restriction enzymes targeting ‘neutral’ elements. Models based on genomic compositions are also effective tools to accurately calculate probabilities of recognition sequences across taxa, and can be applied to species for which reduced-representation data is available (including transcriptomes and ‘neutral’ RAD-seq datasets). The analytical pipeline developed in this study, PredRAD (https://github.com/phrh/PredRAD), and the resulting databases constitute valuable resources that will help guide the design of any study using RAD-seq or related methods.


45
The use of type II restriction enzymes to obtain reduced representation libraries from nuclear 46 genomes, combined with the power of next-generation sequencing technologies, is rapidly becoming one of the most-used commonly strategies to generate genome-wide genotypic and sequence data in both 103 deviations across the eukaryotic tree of life remains unknown. Better models to calculate restriction 104 recognition-sequence probabilities across taxonomic groups are needed to improve the accuracy of 105 predictions of cleavage site frequencies in species without sequenced genomes. These models could be 106 applied using non-genomic datasets (e.g., transcriptomes) to obtain recognition-sequence probability 107 estimates, thus aiding the applicability of RAD-seq methods in non-model organisms.

109
Eukaryotic genomes have heterogeneous compositions with characteristic signatures at the level 110 of di-and trinucleotides that are largely independent of coding status or function (Karlin & Mrázek 1997; 111 Karlin et al. 1998;Gentles & Karlin 2001). Thus, it is possible that genome composition at these levels 112 has a large influence on the abundance of short sequence patterns, such as recognition sequences of 113 restriction enzymes. Models incorporating the information from these genomic compositional signatures 114 should improve the accuracy of restriction recognition-sequence probability calculations.

116
Here we performed systematic in silico genome-wide surveys of genome compositions and 117 recognition sequences, for diverse and commonly used type II restriction enzymes, in 434 eukaryotic 118 whole and draft genomes (Supplementary Table 1) to: (1) characterize restriction recognition-sequence 119 frequencies across the eukaryotic tree of life; (2) explore the potential for predicting genome sizes from 120 restriction recognition-sequence frequency data; (3) develop stochastic models based on genomic 121 compositions to calculate probabilities of recognition sequences across taxa; and (4) evaluate the 122 applicability of these models to species for which only non-genomic data is available (i.e., not whole or 123 draft genome assemblies), such as transcriptomes or restriction-site associated DNA sequence data. The

124
PredRAD analytical pipeline developed in this study (https://github.com/phrh/PredRAD), and the 125 resulting databases constitute a valuable reference resource that will help guide restriction enzyme choice 126 infuture studies using broadly applicable RAD-related methods.

Frequencies of recognition sequences are highly variable across taxa 131
To characterize cleavage site frequencies across the Eukaryotic tree of life, we surveyed 132 restriction recognition-sequences for 18 commonly used palindromic type II restriction enzymes in 434 133 whole and draft genomes. Observed relative frequencies of recognition sequences were highly variable 134 among broad taxonomic groups for the set of restriction enzymes here examined (Table 1)

148
To explore the potential for predicting genome sizes from restriction recognition-sequence 149 frequency data, we modeled their relationship using data from the 434 genomes and 18 restriction

153
Predicted genome sizes, calculated using the linear models with estimated beta parameters for these five 154 enzymes (Table 2)

158
To generate better models for cleavage site frequency calculation in species without sequenced 159 genomes, we developed stochastic models based on the GC content of each genome, as well as the 160 mononucleotide, dinucleotide and trinucleotide compositions to predict the expected frequency of 161 recognition sequences for each restriction enzyme. We evaluated the fit of each model by comparing the 162 in silico observed frequencies of cleavage sites to the expected frequencies predicted by the models using 163 a similarity index (SI), defined as the quotient of the number of observed and expected cleavage sites, 164 minus one. A positive SI indicates that the number of observed cleavage sites is greater than the expected, 165 whereas a negative SI indicates a smaller number of observed sites than expected. If SI is equal to 0, then 166 the number of observed sites is equal to the expectation. For example, a SI = 1 indicates that the number 167 of observed cleavage sites for a particular enzyme in a given genome is twice the number of expected 168 sites predicted by a particular model. Trinucleotide composition models were in general a better predictor, 169 in terms of their accuracy and precision, of the expected number of cleavage sites than any of the other 170 models (Figures 3, 4 and 5). The mononucleotide and GC content models produced relatively poorer 171 predictions that were indistinguishable from one another (Figures 3, 4 and 5). In a few cases the other 172 models outperformed the trinucleotide model, e.g., EcoRI (Figures 3, 4 and 5). The fit of the predictions 173 was highly variable among broad taxonomic groups but generally similar within, e.g., in Neopterigii

177
Recognition-sequence probability can be calculated from non-genomic datasets 178 Genomic resources (whole or draft genomes) are unavailable for most species (Dunn & Ryan

179
[https://bitbucket.org/caseywdunn/animal-genomes] estimate that ~0.015% of species have a published 180 genome to date). However, reduced representation datasets that capture a small fraction of a genome, such 181 as RNA-seq or RAD-seq datasets, can now be easily and economically developed. We investigated the 182 potential use of these datasets to estimate genome composition parameters for our predictive models and 183 calculate recognition-sequence probabilities of any given restriction enzyme. For this we selected a set of 184 27 species out of the 434 examined eukaryotic species with whole and draft genomes, which also have 185 publically available transcriptome data. The restriction-sequence probabilities calculated for the same 186 panel of 18 restriction enzymes, as above, were remarkably similar between those calculated using known 187 composition parameters from the whole and draft genomes and those calculated using estimated 188 composition parameters from transcriptome datasets ( Figure 6). Interestingly, the overall similarity 189 between the two kinds of calculated probabilities (measured as the mean squared error -MSE -190 calculated across all species) was greatest when probabilities were calculated using a mononucleotide 191 composition model (0.046; when MSE=0 the probabilities are identical; MSE value increases as 192 similarity decreases), and decreased when dinucleotide and trinucleotide models were used (0.06 and 193 0.07, respectively). As expected, the species' specific MSE values were variable, and tended to decrease 194 as the propotion of genome represented by the transcriptome increased ( Figure 6).

196
We also calculated recognition-sequence probabilities using parameters estimated from the in  Observed recognition-sequence frequencies for a given restriction enzyme are strikingly variable 224 across broad eukaryotic taxonomic groups, but are similar among closely-related species. This finding is 225 consistent with the hypothesis that the abundance of cleavage sites is largely determined by phylogenetic 226 relatedness. This pattern is most evident in groups that have a larger taxonomic representation, such as 227 mammals. As more genome assemblies become available, patterns within many other underrepresented 228 taxonomic groups will be clarified. Through the use of comparative methods in a robust phylogenetic 229 framework, it will be possible to establish taxon-specific divergence thresholds that could diagnostic of 230 significant evolutionary changes in genome architecture.

232
As expected, observed relative frequencies of cleavage sites with shorter recognition sequences 233 are on average higher than the observed frequencies with longer recognition sequences. However, this 234 pattern in not universal. There are several instances in which the relative frequency of cleavage sites for a 235 high-denomination cutter is higher than that for a low-denomination cutter. For example, in primates the 236 relative frequency of the 8-cutter SbfI (24.6 ± 1.7 RS/Mb) is significantly higher than the frequency of the

244
Predictability of genome sizes 245 For many of the examined type II restriction enzymes examined (e.g. EcoRI, FatI, NsilI, NspI, 246 PciI), the observed frequencies of recognition sequences in eukaryotic genomes are consistent with the 247 idea that they behave neutrally, evolutionary speaking, and therefore can be readily used as parameters in 248 linear models to estimate genome sizes (Figures 1 and 2). In contrast, the observed frequencies of

279
It is not surprising that the fit of the predictions made by the models is highly variable across

286
The remarkable similarity between probabilities calculated using parameters estimated from non-

306
For the design of a study using RAD-seq, or a related methodology, there are two fundamental 307 questions that researchers commonly face: (1) what is the best restriction enzyme to use to obtain a 308 desired number of RAD tags in the organism of interest? and (2) how many markers can be obtained with 309 a particular enzyme in the organism of interest? The results from this study coupled with the developed 310 software pipeline PredRAD, will allow any researcher to obtain an approximate answer to these 311 questions. The flow diagram in Figure 9 illustrates a suggested workflow.

313
In a best-case scenario for the practical design of a study using RAD-seq, or a related 314 methodology, the species of interest is already included in the database presented here. In this case, the 315 best proxy for the estimated number of RAD tags that could be obtained empirically would be twice the

331
In the scenario that the genome sequence of the species of interest is not available, the genome 332 size and restriction recognition-sequence probability for the enzyme(s) of interest can be estimated to 333 obtain an approximation of recognition-sequence frequencies (absolute numbers). Our observations 334 demonstrate that a genome-size range can be estimated by applying linear regression models to the 335 number of markers obtained in an empirical RAD-seq experiment using a restriction enzyme targeting a 336 'neutral' recognition sequence (e.g., EcoRI, NsilI, NspI, PciI; we advise caution using 4-cutter enzymes 337 as in some taxa they can have cleavage frequencies that may effectively lead to sequencing the whole

360
CONCLUSIONS 361 362 In this study, we performed systematic in silico genome-wide surveys of genome compositions 363 and recognition sequences, for diverse and commonly used type II restriction enzymes across the 364 eukaryotic tree of life. Our observations reveal that recognition-sequence frequencies for a given 365 restriction enzyme are strikingly variable among broad eukaryotic taxonomic groups, being largely 366 determined by phylogenetic relatedness. We demonstrate that genome sizes can be predicted from 367 cleavage frequency data obtained with restriction enzymes targeting 'neutral' recognition sequences.

368
Stochastic models based on genomic compositions are also effective tools to accurately calculate 369 probabilities of recognition sequences across taxa, and can be applied to species for which reduced-370 representation genomic data is available (including transcriptomes and 'neutral' RAD-seq datasets). The 371 results from this study and the software developed from it will help guide the design of any study using 372 RAD sequencing and related methods. As more genome assemblies become available in underrepresented 373 taxonomic groups, the patterns of compositional biases and restriction site frequencies across the 374 eukaryotic tree of life will become clearer and will improve our understanding of genome evolution.  (1) Where x is the number of restriction recognition-sequences in the genome, and β 0 and β 1 are the estimated 419 parameters. Table 2

432
Mononucleotide and GC content sequence models were used to estimate the probability of a 433 particular recognition sequence assuming that each nucleotide is independent of the others and of its    Table 2). We also used the data from the in silico RAD sequencing experiments (described above) as 465 reduced representation datasets for these species. We estimated genome composition parameters from 466 transcriptome and RAD-seq datasets, and calculated recognition-sequence probabilities using the models.

469
To assess the effectiveness of the predictive recognition sequence models, we compared the

478
To measure the overall similarity between the restriction recognition-sequence probabilities 479 calculated using known composition parameters from the genome and those calculated using estimated 480 composition parameters from reduced-representation transcriptome and genome datasets, we calculated 481 the mean squared error (MSE) per species.

483
(4) Where !! ! is the probability of a restriction recognition-sequence (of an enzyme i) calculated 486 using composition parameters estimated from reduced-representation datasets, !" ! is probability of a 487 restriction recognition-sequence calculated using known composition parameters from genome datasets.

488
Each enzyme was assigned an arbitrary number from 1 to 18 (n). When MSE=0 the probabilities are 489 identical. MSE value increases as similarity decreases.

492
To evaluate the possibility that recognition-sequence frequency patterns inconsistent with

509
Supplementary figures 1-9, and tables 1 and 2 are available as supplementary material.

606
Biases in most of these dinucleotides are likely linked to important biological processes. Notably the 607 underrepresented dinucleotide CG is a widely known target for methylation related to transcriptional 608 regulation (Bird 1980) and retrotransposon inactivation (Yoder et al. 1997) in vertebrates and eudicots.

609
The corresponding overrepresentation of AG/CT fits the classic model of "methylation-deamination-

632
The largest and more widespread overrepresentation biases were for the trinucleotide pair AAA/TTT, 633 being significant in most eukaryotes, except for the majority of the Dikarya fungi ( !"# * =1.18 ± 0.07). 634 The trinucleotide pairs TAA/TTA and AAT/ATT were significantly overrepresented in many metazoan