No signal of deleterious mutation accumulation in conserved gene sequences of extant asexual hexapods

Loss of sex and recombination is generally assumed to impede the effectiveness of purifying selection and to result in the accumulation of slightly deleterious mutations. Empirical evidence for this has come from several studies investigating mutational load in a small number of individual genes. However, recent whole transcriptome based studies have yielded inconsistent results, hence questioning the validity of the assumption of mutational meltdown in asexual populations. Here, we study the effectiveness of purifying selection in eight asexual hexapod lineages and their sexual relatives, as present in the 1 K Insect Transcriptome Evolution (1KITE) project, covering eight hexapod groups. We analyse the accumulation of slightly deleterious nonsynonymous and synonymous point mutations in 99 single copy orthologue protein-coding loci shared among the investigated taxa. While accumulation rates of nonsynonymous mutations differed between genes and hexapod groups, we found no effect of reproductive mode on the effectiveness of purifying selection acting at nonsynonymous and synonymous sites. Although the setup of this study does not fully rule out nondetection of subtle effects, our data does not support the established consensus of asexual lineages undergoing ‘mutational meltdown’.

Oenothera evening primroses and Boechera rockcress, it was absent in Lineus ribbon worms as well as four aphid species and, opposite to predictions, reduced in asexual as compared to sexual oribatid mites [14][15][16][17][18][19] . Moreover, all whole transcriptome based studies found excessive variation among genes, and hence doubts have been raised about the robustness of inferences drawn from single gene analyses 20 . These conflicting results highlight the need for more studies of larger gene sets along with broader taxonomic sampling to infer whether or not accumulation of deleterious mutations is indeed a consequence of asexual reproduction.
Here, we study the effectiveness of purifying selection in obligately asexual and related sexual hexapod lineages covering eight hexapod groups using transcriptome data generated by the 1KITE project (1 K Insect Transcriptome Evolution, www.1kite.org). As parthenogenesis is a lineage-level trait, we use the term 'lineage' , whenever to distinguish between the two reproductive modes (sexual and asexual) within a hexapod group, from here on. With the term 'hexapod group' we refer to one of the following analysed taxa: Collembola, Zygentoma, Phasmatodea, Mantodea, Thysanoptera, Sternorrhyncha, Hymenoptera and Psocodea, respectively.
We compared the accumulation of slightly deleterious mutations in nuclear protein-coding genes between eight asexual lineages and their sexual relatives as available from 1KITE (see Fig. 1). To this end, we first inferred divergence at nonsynonymous sites normalised for background mutation rates (dN/dS), and the potential 'deleteriousness' of nonsynonymous substitutions. Second, we investigated selection on Codon Usage Bias (CDC), because selection also acts at synonymous sites 21 . We based all analyses on 99 single copy orthologues which we found to be under purifying selection in the analysed lineages. We found extensive variation in dN/dS and CDC among genes and between hexapod groups, but no overall difference between reproductive modes. Methods species selection. We searched for obligately asexual lineages within the 1KITE species list based on primarily van der Kooi & Schwander 2014 22 33 . The taxon sampling includes eight asexual lineages along with their closest sexual relatives as present in 1KITE, covering eight hexapod groups. Sexual and asexual lineages are depicted in red and blue, respectively. Silhouettes courtesy of Hans Pohl.
www.nature.com/scientificreports www.nature.com/scientificreports/ were used as input for identification of orthologous sequences among all studied species. For information on References, BioSample IDs and Bioproject IDs of assembled transcriptomes, see Supplementary Table S1. orthologue detection and alignment processing. To infer orthologue protein-coding genes among the 16 hexapod species, we used the Orthologous MAtrix (OMA) version 2.1.1 installed on the Vital-IT computing cluster 36 . To generate the input for OMA, we extracted long Open Reading Frames (ORFs) and predicted likely coding regions with Transdecoder version 2.0.1 using default options for each transcript and species 37 . Amino acid sequences of likely coding regions of each of the 16 species were passed to OMA together with an unrooted cladogram built by reducing a published phylogeny of hexapods to the eight hexapod groups analysed in this study 33 , which resulted in 286 orthologues shared among all analysed species. We aligned the amino acid sequences (ORFs) of each of the 286 orthologues with M-Coffee using a combination of different alignment methods creating a consensus alignment (option clustalw2_msa muscle_msa kalign_msa mafftgins_msa t_coffee_msa) 38 . Based on these, we generated corresponding codon alignments from the original nucleotide sequences with T-Coffee version 11.00.8 38 . As orthologue detection based on naturally incomplete datasets such as transcriptome data can lead to the detection of false positives, e.g. due to loci with paralogous sequences, we tested all protein sequences of the 286 orthologues of all 16 species for their presence in a precompiled set of orthologues of insects (insecta_hmmer3.1; www.deep-phylogeny.org/hamstr/download/datasets/hmmer3) using HaMStR version 13.2.6 39 . We only kept multiple sequence alignments of orthologues for which all species had an equivalent in the precompiled orthologue set leaving a set of 153 loci for further analyses. We curated these multiple sequence alignments using Gblocks version 0.91b with sequence type set to codons (t = c) and minimum block length set to 4 (b4 = 4) 40 .
Accumulation of nonsynonymous mutations. As a measure of purifying selection, we analysed the divergence at nonsynonymous sites normalised for background substitution rates (dN/dS) using CodeML as implemented in PAML version 4.9 41 . For this, we first manually constructed an unrooted cladogram of the 16 analysed species based on a published hexapod phylogeny 33 . To exclude orthologues that were under positive selection, we ran BUSTED as implemented in HYPHY version 2.3.10 42 providing as input the unrooted cladogram described above and the multiple sequence alignments of the 153 orthologues. We found 54 loci showing signatures of episodes of positive selection, which left 99 loci for subsequent analyses (for GenBank Accession numbers, see Supplementary Table S2). CodeML relies on a Maximum Likelihood framework to estimate the goodness of fit of a codon substitution model to a sequence alignment and an unrooted species tree with gene-specific branch lengths for inference of branch-specific dN/dS ratios. We calculated per-gene branch lengths for the unrooted cladogram based on each of the 99 codon alignments, accordingly, using RAxML version 8.2.8 43 with GTRGAMMAI (with four GAMMA rate categories) set as model of sequence evolution. We modified a custom script used by Brandt et al. 15 (see Supplementary information) to pass the loci-specific branch lengths and fixed species tree together with each codon alignment to CodeML for divergence rate estimations. Due to the 1KITE taxon sampling, the asexual and sexual lineages of each hexapod group used in our analysis are likely not natural sister lineages (i.e. not the closest extant relatives) and therefore the time of transition to asexuality (or the split from the closest sexual relative) was unknown. As this could potentially lead to an overestimation of dN/dS under a model restricted to only two rates (one for asexual and one for sexual branches), we chose a free model allowing for different dN/dS ratios, one for each branch in the tree. For between-species comparisons of dN/dS ratios, we excluded all dN/dS ratios of internal branches and four terminal branch dN/ dS ratios that were > 1, indicating positive selection acting at one branch, from statistical analyses. We then tested whether branch-specific dN/dS ratios differed according to (I) gene (II) reproductive mode, or (III) hexapod group using a permutation ANOVA with 5,000 bootstrap replicates (available at https://gist.github.com/ KamilSJaron/358c997698b67486be47d4e8eef2921d) 44,45 . Differences in dN/dS can be driven by differing synonymous substitution rates and, in the long-term, different levels of saturation at synonymous sites. Given the old age of the splits between sexual and asexual lineages (~40 myo for Phasmatodea -~160 myo for Zygentoma) 33,46 we tested for differences in branch-specific dS as described above. To infer whether or not sexual and asexual lineages within individual hexapod groups differed in dN/dS ratios, we compared between reproductive modes using Wilcoxon signed-rank tests. All statistical analyses were done in R version 3.4.4 47 .
'Deleteriousness' of nonsynonymous mutations. To infer the 'deleteriousness' of nonsynonymous substitutions we analysed hydrophobicity changes from ancestral to replacement amino acids along the terminal branches of the phylogenetic tree. Hydrophobic interactions are the main determinants of the 3D conformation of proteins and thus an indicator of protein stability 48 . Inference of ancestral amino acids relies on the presence of an outgroup sequence included in the input amino acid alignments and the phylogenetic tree (cladogram with loci-specific branch lengths) used in analyses with CodeML. Therefore, we first searched for orthologues shared among all 16 hexapod species plus the crustacean Xibalbanus cf. tulumensis (previously Speleonectes cf. tulumensis), as a representative of Remipedia, the sister-group of hexapods 32,33 . For this, we predicted ORFs from the assembled transcriptome of X. cf. tulumensis and checked the ORFs of X. cf. tulumensis for presence of orthologues in the precompiled orthologue set of insects as described above. We found 73 ORFs of X. cf. tulumensis, each of them orthologous to one of the previously inferred 99 clusters of orthologues of the 16 hexapod species (for GenBank Accession numbers, see Supplementary Table S2). We again aligned the amino acid sequences (ORFs, including now the sequences of X. cf. tulumensis) of the 73 clusters of orthologues and subsequently generated corresponding codon alignments from the original nucleotide sequences as described above. Further, we added X. cf. tulumensis to the unrooted manually constructed cladogram used for the calculation of branch lengths for dN/dS ratio analyses as outgroup to all hexapods. Using this tree as a fixed topology, we calculated branch lengths for each orthologue locus individually as described above for analyses of nonsynonymous www.nature.com/scientificreports www.nature.com/scientificreports/ mutation accumulation, and then translated the curated alignments into amino acids using EMBOSS version 6.6.0 49 . To predict ancestral amino acid sequences for each internal node in each inferred tree from each gene, we passed the 73 amino acid codon alignments individually with its respective species trees (and estimated loci-specific branch lengths) to CodeML using the modified custom script mentioned above (see Supplementary information). We determined the strength of hydrophobicity changes (Hydrophobicity Scores; HS) for each of the amino acid transitions along the terminal branches of the species tree using a hydrophobicity scoring (HS) matrix 44,50 . HS indicates the 'deleteriousness' of a nonsynonymous mutation by measuring the strength in hydrophobicity change from ancestral to replacement amino acid. The lower the HS, the stronger is the change in hydrophobicity and, hence, the deleteriousness of the underlying nonsynonymous mutation. We compared values of HS between the two reproductive modes using Generalised Linear Mixed Models (GLMM) implemented in the R package lme4 with gene nested in species set as random effect, correction for overdispersion and Poisson distribution fitting 51 .
Accumulation of synonymous mutations. Synonymous mutations are generally assumed to be neutral but can be subject to purifying selection because different codons can influence the speed and accuracy of translation 21 . Hence, we also analysed the effectiveness of selection acting on Codon Usage Bias (CUB). First, we inferred the existence of CUB for each of the 99 orthologues of the 16 species using the Effective Number of Codons (Nc) as measure with the software codonW version 1.4 52,53 . The Nc specifies the deviation of observed codon usage from equal usage of all codons ranging from 20 (each amino acid is encoded by one codon only; strong CUB) to 61 (equal use of all possible codons; no CUB). Next, we inferred selection on CUB using the Codon Deviation Coefficient (CDC) 54 . Measurement of CDC allows for cross-species comparisons by correcting for background nucleotide composition and is particularly robust because, unlike dN/dS ratio analyses, it does not rely on likelihood and branch length estimates. The CDC represents the deviation of expected CUB based on observed positional GC and purine contents from observed CUB, ranging from 0 (no deviation; no detectable selection on CUB) to 1 (maximum deviation; strong selection on CUB). We estimated CDC for the processed alignments using Composition Analysis Toolkit version 1.3 54 and analysed it for an effect of I) gene, II) reproductive mode, and III) hexapod group and inferred within-group differences as described above for statistical analyses of dN/dS ratios.

Results
Accumulation of nonsynonymous mutations. We estimated nonsynonymous to synonymous divergence (dN/dS) along individual branches of a phylogenetic tree comprising eight asexual and eight sexual hexapod species for 99 single-copy orthologous protein-coding genes under purifying selection (see Fig. 1; Methods). In genes under purifying selection, nonsynonymous mutations have likely deleterious effects, hence, a higher dN/dS ratio indicates less effective purifying selection 55 . Consistent with the expectation for loci under purifying selection, dN/dS ratios at terminal branches were generally low (mean dN/dS = 0.032). Contrasting the established consensus on deleterious mutation accumulation in asexual lineages, there was no difference in dN/ dS when compared between sexual and asexual branches over all hexapod groups (gene effect P < 0.001, reproductive mode effect P = 0.488, hexapod group effect P = 0.048, interaction reproductive mode * hexapod group P = 0.145; Permutation ANOVA). The absence of a difference between reproductive modes was not driven by a difference in dS (gene effect P < 0.001, reproductive mode effect P = 0.278, hexapod group effect P < 0.001, interaction reproductive mode * hexapod group P = 0.004; Permutation ANOVA). There was significant among-gene variation in dN/dS (range 0-0.454; variance 1.04 * 10 −3 ) and a significant difference in dN/dS among hexapod groups. To detect whether or not the effectiveness of purifying selection differed between reproductive modes within each hexapod group, we compared dN/dS between reproductives modes in each hexapod group on a pergene basis. For Zygentoma, we found significantly lower per-gene dN/dS in the sexual as compared to the asexual terminal branch indicating more effective purifying selection for the sexual branch (see Table 1; Fig. 2a red box).
'Deleteriousness' of nonsynonymous mutations. In addition to nonsynonymous mutation accumulation in asexual hexapod lineages, purifying selection is expected to lead to more deleterious amino acid transitions in asexual lineages as reflected by stronger hydrophobicity changes from ancestral to replacement amino acids along asexual and sexual terminal branches 56 . Contrasting this expectation, but consistent with our dN/dS estimates, HS were similar between asexual and sexual branches (z = −0.152; P = 0.879; GLMM; see P < 0.001 P < 0.001 P = 0.025 P = 0.682 P = 0.129 P < 0.001 P = 0.99 P = 0.361 Table 1. V and P values of within-hexapod group comparisons of dN/dS and CDC between reproductive modes. Values were inferred by comparing dN/dS and CDC per gene between the sexual and asexual lineages in each of eight hexapod groups using Wilcoxon signed-rank tests. Underlined and bold measures indicate more effective purifying selection in sexual species (underlined) and asexual species (bold), respectively (see Fig. 2a,b) www.nature.com/scientificreports www.nature.com/scientificreports/ Supplementary Fig. S1). Percentages of ancestral to replacement amino acid transitions with more dissimilar hydrophobicity (HS < 90) were similar between asexual and sexual branches (44.936% and 44.893% of all asexual and sexual transitions, respectively).

Accumulation of synonymous mutations.
We assessed whether or not purifying selection on synonymous sites was less effective in asexual as compared to sexual hexapod lineages by inferring selection on CUB. To infer whether or not the species investigated use some codons preferentially over others, we first analysed the Effective Number of Codons (Nc). Nc ranges from 20 (each amino acid is encoded by one codon only; strong CUB) to 61 (equal use of all possible codons; no CUB). Overall, all species showed CUB, with Nc means of species ranging from 40.6 in D. vitifoliae (Sternorrhyncha) to 54.1 in F. candida (Collembola; see Supplementary Fig. S2) which is comparable to values of Nc found in other invertebrate groups, e.g. bivalves and nematodes 57,58 . Afterwards, we directly inferred selection on CUB using CDC as a measure. CDC calculates the deviation of observed from The boxplots show the distribution of per-gene differences in dN/dS between the sexual and asexual terminal branch (Δ sex-asex dN/dS) and of per-gene differences in CDC between the sexual and asexual species (Δ sex-asex CDC) of each of eight hexapod groups for 99 orthologues. For better representation, the ordinate is restricted to between −0.06 and 0.06 and −0.1 and 0.1, respectively, and outliers outside of 1.5 times the interquartile range (whiskers) are excluded. For dN/dS seven of eight within-hexapod-group comparisons between sexual and asexual branches yielded non-significant results (see Table 1). The red plot denotes a significantly lower per-gene dN/dS in sexual as compared to asexual terminal branches of Zygentoma (V = 3786, P < 0.001***; Wilcoxon signed-rank test). For CDC four of eight within-hexapod-group comparisons between sexual and asexual branches yielded significant results (see Table 1). There was significantly higher per-gene CDC in sexual as compared to asexual Collembola and Phasmatodea (red plots; V = 1191, P < 0.001*** and V = 1832, P = 0.025*, respectively; Wilcoxon signed-rank test) contrasting a significantly lower per-gene CDC in sexual as compared to asexual Zygentoma and Sternorrhyncha species (blue plots; V = 4198, P < 0.001 and V = 3562, P < 0.001, respectively; Wilcoxon signed-rank test). (2019) 9:5338 | https://doi.org/10.1038/s41598-019-41821-x www.nature.com/scientificreports www.nature.com/scientificreports/ predicted CUB by correcting for background nucleotide composition ranging from 0 (no selection on CUB) to 1 (effective selection on CUB). Consistent with the results from analyses of nonsynonymous mutation accumulation and 'deleteriousness' of nonsynonymous mutations, there was no difference in per-gene CDC between reproductive modes (gene effect P < 0.001, reproductive mode effect P = 0.283, hexapod group effect P < 0.001, interaction reproductive mode * hexapod group P < 0.001; Permutation ANOVA). Further, there was significant among gene variation for CDC (range 0.064-0.329; variance 1.653 * 10 −3 ) and a significant difference between hexapod groups but, unlike for dN/dS estimates, there was also a significant interaction between hexapod group and reproductive mode. As for dN/dS, we compared per-gene CDC between reproductive modes in each hexapod group, individually. Four of eight within-hexapod-group comparisons between sexual and asexual species yielded significant results (see Table 1; Fig. 2b red and blue boxes, respectively). For Collembola and Phasmatodea, there was a significantly higher per-gene CDC in sexual as compared to asexual species indicating more effective selection on CUB in sexual species, whereas for Zygentoma and Sternorrhyncha there was a significantly lower per-gene CDC in sexual as compared to asexual species indicating more effective selection on CUB in asexual species.

Discussion
It has become established consensus among evolutionary biologists that sex and recombination increase the effectiveness of purifying selection, based on theoretical considerations and empirical evidence derived from a multiplicity of studies 12 . Our results do not match these studies: we find no evidence for accumulation of deleterious mutations in asexual hexapod species. Overall, asexual and sexual lineages neither differed in nonsynonymous site divergences corrected for background substitution rates (dN/dS), in synonymous substitution rates (dS) potentially influencing dN/dS ratios, in the 'deleteriousness' of nonsynonymous mutations nor in selection acting on CUB (CDC). Comparisons between sexual and asexual lineages within each group differed for several hexapod groups which may hint at group-specific differences in the consequences of asexuality on effectiveness of purifying selection. However, here, we refrain from drawing conclusions based on the within-group comparisons because the taxon sampling of 1KITE did not allow for analysing multiple within-group replicates but restrict the discussion to the overall result of more effective selection being absent in our data of sexual hexapods: which (non-mutually exclusive) mechanisms might be responsible for the discrepancy between this finding and the established consensus?
First, analyses of purifying selection acting at nonsynonymous sites and on their 'deleteriousness' may have been affected due to data limitations. dN/dS ratio analyses and ancestral state reconstructions rely on branch length estimates (see Methods). Due to the limited nature of the taxon sampling of the 1KITE data set for this study, the asexual and sexual species analysed are most likely not sister species, but rather more distantly related 46,59 . Therefore, the loss of sex did not occur with the split of the sex-asex species pair as present in the given phylogenetic tree here, such that evolution over some fraction of the asexual branch was likely sexual. If mutations occuring in asexual lineages did not accumulate at greatly increased rates compared with sexuals, a change in deleterious mutation accumulation along with the transition to asexuality might be masked by the rates that occured in the sexual fraction of the branch. Thus, the power to detect an effect of reproductive mode on the effectiveness of purifying selection might be low. Further, the orthologue search among phylogenetically distantly related hexapod groups and the stringent control for false positive orthologues and loci under positive selection resulted in a rather small orthologue set for analysis (99 orthologue loci). This biases the analyses towards strongly conserved loci and excludes recently evolved orthologues which might differ in accumulation of deleterious mutations between reproductive modes. Also, the within hexapod group comparisons between reproductive modes for CDC did not resemble those for dN/dS; in the case of Zygentoma they even opposed them (see Fig. 2a,b; Table 1). This is surprising, given that translational selection acting at synonymous sites is assumed to be weak and effective purging of synonymous mutations likely occurs at lower rates as compared to that of nonsynonymous mutations. Hence, within hexapod group comparisons of purifying selection acting at nonsynonymous sites should not oppose those at synonymous sites.
Second, the analysed asexual lineages may have lost sex too recently to have fixed enough deleterious mutations to be detected. This was for example also assumed to be responsible for absence of deleterious mutation accumulation observed in whole transcriptome data of fissiparous Lineus ribbon worms and four aphid species 14,18 . In fact, the occurrence of abundant males in some locally restricted populations of N. phytophila, E. californica, and L. bostrychophila analysed in this study is in line with a rather recent loss of sex 18 . Additionally, interference of rare, furtive, or cryptic sex with asexual genome evolution may explain the observed absence of less effective selection in asexual hexapods because rare events of sex are assumed to be sufficient to compensate for predicted consequences of asexuality 60 .
A third reason for the absence of deleterious mutation accumulation in the analysed data of asexual hexapod lineages may be effective homogenising mechanisms. Gene conversion and DNA repair have been shown to maintain DNA integrity, e.g. within the human Y chromosome, higher plant chloroplasts and animal mitochondria [61][62][63] . If homogenising mechanisms play a role in the analysed species remains to be investigated.
Fourth, besides reproductive mode, population size acts as a major determinant of effectiveness of selection with the speed of mutation accumulation being inversely related to population size, as shown by modelling approaches 9,64,65 . Large population sizes have been suggested to maintain effective purifying selection in asexual organisms, e.g. in oribatid mites and polyphagous pest insects, such as scale insects 15,66,67 . Further, large population sizes have been shown to increase the effectiveness of selection acting on CUB in a variety of animal species and on nonsynonymous mutations with large deleterious effects in Caenorhabditis elegans 68,69 . As information on population sizes of the analysed species is absent, any correlation between population sizes and the observed absence of more effective purifying selection acting at nonsynonymous and synonymous sites has to remain speculative.
In conclusion, our results do not support the established consensus of reduced effectiveness of purifying selection in asexual species, contrasting earlier analyses of few individual genes in a variety of animal and plant species. However, whether or not this is due to the limited nature of our data or due to non-mutually exclusive biological www.nature.com/scientificreports www.nature.com/scientificreports/ mechanisms has to remain elusive at this point. Future studies, hence, need not only to include large gene sets but also carefully selected closely related sexual and asexual focal lineages to study the effectiveness of purifying selection in asexual organisms in more detail.