Multiple hybridization events and repeated evolution of homoeologue expression bias in parthenogenetic, polyploid New Zealand stick insects

Abstract

During hybrid speciation, homoeologues combine in a single genome. Homoeologue expression bias (HEB) occurs when one homoeologue has higher gene expression than another. HEB has been well characterized in plants but rarely investigated in animals, especially invertebrates. Consequently, we have little idea as to the role that HEB plays in allopolyploid invertebrate genomes. If HEB is constrained by features of the parental genomes, then we predict repeated evolution of similar HEB patterns among hybrid genomes formed from the same parental lineages. To address this, we reconstructed the history of hybridization between the New Zealand stick insect genera Acanthoxyla and Clitarchus using a high-quality genome assembly from Clitarchus hookeri to call variants and phase alleles. These analyses revealed the formation of three independent diploid and triploid hybrid lineages between these genera. RNA sequencing revealed a similar magnitude and direction of HEB among these hybrid lineages, and we observed that many enriched functions and pathways were also shared among lineages, consistent with repeated evolution due to parental genome constraints. In most hybrid lineages, a slight majority of the genes involved in mitochondrial function showed HEB towards the maternal homoeologues, consistent with only weak effects of mitonuclear incompatibility. We also observed a proteasome functional enrichment in most lineages and hypothesize this may result from the need to maintain proteostasis in hybrid genomes. Reference bias was a pervasive problem, and we caution against relying on HEB estimates from a single parental reference genome.

Keywords: Phasmatodea; asexual; gene expression; hybrid; transcriptomics.

FIGURE 1

Hi‐C contact map for all scaffolds. The red squares represent the 18 super‐scaffolds.

FIGURE 2

Gene trees reconstructed from: (a) phase block 13.0 from the serine/threonine‐protein kinase SMG1 gene; (b) phase block 352.0 from the serine proteinase stubble gene. (c) Gene tree estimated from whole mitochondrial genomes. The tips in the mitochondrial gene tree show reproductive mode, inferred ploidy and hybrid origins. Branch lengths for gene trees are drawn proportional to the expected number of substitutions per site following the scale bars. Numbers above branches are maximum likelihood bootstrap percentages where values <50% were omitted. Individual alleles sampled from Acanthoxyla individuals are coloured orange and those sampled from Clitarchus are coloured blue. Codes in parentheses at tip labels refer to individual phased alleles.

FIGURE 3

Dsuite plots showing the f _b‐ratios calculated using 238,791,482 genome‐wide SNPs. The f _b‐ratios are summarized using the topologies from: (a) whole mitochondrial genomes; (b) a network inferred from phased BUSCO gene regions.

FIGURE 4

Networks reconstructed from gene trees using the InferNetwork_MP method in PhyloNet: (a) constructed assuming Acanthoxyla prasina (AXP), A. prasina (AXG) and A. inermis are hybrid lineages with a maximum of three reticulations; (b) constructed assuming A. prasina (AXX) is a hybrid lineage with a maximum of one reticulation. Acanthoxyla lineages are coloured orange, Clitarchus species are coloured pale blue and hybridization events are in dark blue.

FIGURE 5

The distribution of phase block lengths (bp) from the 151 regions passing filtering.

FIGURE 6

PCA plots: (a) SNP variation from RNAseq data; (b) HEB signal towards the Clitarchus hookeri subgenome. Both analyses were performed using the C. hookeri reference assembly. The HEB PCA plot contains only those lineages for which HEB analysis was performed.

FIGURE 7

Plots showing the frequency distribution of Clitarchus hookeri alleles in the Acanthoxyla prasina (AXP) transcriptome, which has a history of hybridization with C. hookeri. The allele frequencies are measured using three reference genomes: (a) C. hookeri; (b) Acanthoxyla nov. sp. 1; (c) Argosarchus horridus (c). The red dotted lines represent the threshold for alleles to be inferred showing balanced expression (between the lines) or HEB towards the C. hookeri (right of the right‐hand line) or Acanthoxyla (left of the left‐hand line) subgenomes.

FIGURE 8

Heatmaps showing the direction and magnitude of HEB in genes from three KEGG pathways: (a) oxidative phosphorylation; (b) ribosome; (c) proteasome; (d) all genes showing a significant HEB signal. A positive value (blue) represents HEB towards the Clitarchus hookeri subgenome while a negative value (red) represents the HEB towards the Acanthoxyla subgenome.

FIGURE 9

UpSet plots showing the sharing of enriched GO categories among Acanthoxyla lineages: (a) GO categories enriched towards the Clitarchus hookeri subgenome; (b) towards the Acanthoxyla subgenome. The set size refers to number of enriched GO categories for each lineage and the intersection size is the number of enriched GO categories shared among lineages.