Gene transcriptional profiles in gonads of Bacillus taxa (Phasmida) with different cytological mechanisms of automictic parthenogenesis

Abstract

The evolution of automixis - i.e., meiotic parthenogenesis - requires several features, including ploidy restoration after meiosis and maintenance of fertility. Characterizing the relative contribution of novel versus pre-existing genes and the similarities in their expression and sequence evolution is fundamental to understand the evolution of reproductive novelties. Here we identify gonads-biased genes in two Bacillus automictic stick-insects and compare their expression profile and sequence evolution with a bisexual congeneric species. The two parthenogens restore ploidy through different cytological mechanisms: in Bacillus atticus, nuclei derived from the first meiotic division fuse to restore a diploid egg nucleus, while in Bacillus rossius, diploidization occurs in some cells of the haploid blastula through anaphase restitution. Parthenogens' gonads transcriptional program is found to be largely assembled from genes that were already present before the establishment of automixis. The three species transcriptional profiles largely reflect their phyletic relationships, yet we identify a shared core of genes with gonad-biased patterns of expression in parthenogens which are either male gonads-biased in the sexual species or are not differentially expressed there. At the sequence level, just a handful of gonads-biased genes were inferred to have undergone instances of positive selection exclusively in the parthenogen species. This work is the first to explore the molecular underpinnings of automixis in a comparative framework: it delineates how reproductive novelties can be sustained by genes whose origin precedes the establishment of the novelty itself and shows that different meiotic mechanisms of reproduction can be associated with a shared molecular ground plan.

Keywords: Gene expression; Parthenogenesis; Phylostratigraphy; RNA-seq; dN/dS.

Fig. 1

Overview of the oogenesis processes in a the bisexual Bacillus grandii maretimi (blue) and in the two parthenogens, b Bacillus atticus (yellow) and c Bacillus rossius rossius (red)

Fig. 2

Phylogenetic relationships of Bacillus taxa and phylostratigraphy of gonads-biased genes. a Schematic drawing of the species tree obtained from Mantovani et al. (2001), with indication of the reproductive strategy for each taxon. Underlined species are those considered in the present study. b Phylostratigraphy of gonads upregulated genes; numbers on the right represent the number of transcripts present in each subset / assembly

Fig. 3

PCA analysis on the TMM normalized gene expression values for a all genes in gonads samples, b gonads-biased genes in gonads samples, c all genes in legs samples, (b) legs-biased genes in legs samples

Fig. 4

Gene expression patterns of parthenogens gonads compared to the bisexual species ones: a comparison between the automictic Bacillus rossius rossius and the bisexual Bacillus grandii maretimi, b comparison between the automictic Bacillus atticus and the bisexual Bacillus grandii maretimi. Numbers in parentheses are the genes with gonads-biased or non-biased expression regulation, respectively

Fig. 5

Venn-diagram representation of gonads genes expression patterns in comparison to the bisexual species. Overlaps represent parthenogens’ genes with a shared pattern in the two parthenogens, for which Fisher’s exact test has been used to determine p and odds ratio in comparison to the genomic background (n = 15,972)

Fig. 6

Parallel instances of positive selection: a on the left, number of sites inferred to have undergone positive selection (BEB > 95) in the two parthenogens; on the right, TMM-normalized expression values for the same genes. Genes with < 10 TPMs in parthenogen gonads were not considered. On the y axis, the orthogroups in which the genes are found, and the putative homologs (in bold) are given. b Venn-diagrams representation of the shared changes across the two parthenogens, for which Fisher’s exact test has been used to determine p and odds ratio in comparison to the genomic background (n = 15,972)