Silencing of the DNA methyltransferase 1 associated protein 1 (DMAP1) gene in the invasive ladybird Harmonia axyridis implies a role of the DNA methyltransferase 1‐DMAP1 complex in female fecundity

The invasive harlequin ladybird Harmonia axyridis is a textbook example of polymorphism and polyphenism as the temperature during egg development determines the frequency of melanic morphs and the number and size of black spots in nonmelanic morphs. Recent concepts in evolutionary biology suggest that epigenetic mechanisms can translate environmental stimuli into heritable phenotypic changes. To investigate whether epigenetic mechanisms influence the penetrance and expressivity of colour morphs in H. axyridis, we used RNA interference to silence key enzymes required for DNA methylation and histone modification. We found that neither of these epigenetic mechanisms affected the frequency of different morphs, but there was a significant impact on life‐history traits such as longevity and fecundity. Strikingly, we found that silencing the gene encoding for DNA methyltransferase 1 associated protein 1 (DMAP1) severely reduced female fecundity, which correlated with an abundance of degenerated ovaries in DMAP1‐knockdown female beetles. Finally, we observed significant differences in DMAP1 expression when we compared native and invasive H. axyridis populations with a biocontrol strain differing in egg‐laying capacity, suggesting that the DNA methyltransferase 1‐DMAP1 complex may influence the invasive performance of this ladybird.


Introduction
Individuals of the same species can often be assigned to one of several distinct phenotypes, a phenomenon known as polymorphism when the variation depends on genetic differences and polyphenism when an identical genome can give rise to different phenotypes according to the environment (Simpson et al., 2011). Such phenotypic plasticity is thought to be advantageous in dynamic environments (Levins, 1968;DeWitt et al., 1998), particularly if the adaptations prepare the organism for future environmental changes. The harlequin ladybird Harmonia axyridis (also known as the multicoloured Asian ladybird) is a textbook example of polymorphism because it occurs as a number of distinct morphs differing, for example, in overall colour and the number of the spots on the elytra (Dobzhansky, 1924(Dobzhansky, , 1933. These polymorphic colour patterns are regulated by a single locus, which comprises the four major alleles axyridis, conspicua, spectabilis and succinea, as well as 11 rare alleles (Hosino, 1940;Komai, 1956;Tan, 1946;Tan and Li, 1934). The melanic forms (axyridis, conspicua and spectabilis) may be advantageous in cold climates when dark surfaces absorb heat more quickly during exposure to sunlight. The frequency of different morphs therefore varies amongst populations with different geographical ranges, but it also exhibits seasonal changes (Jiang et al., 2008;Jing and Zhang, 2001;Yuan et al., 1994). The latter occurs because the temperature during the larval development stages influences the frequency of melanic morphs (Knapp and Nedvěd, 2013). H. axyridis First published online 10 September 2019. therefore combines the classic features of both polymorphism and polyphenism (Michie et al., 2010).
It is currently unclear how environmental stimuli translate into distinct phenotypes encoded by the same genome (Vilcinskas and Vogel, 2016). Recent concepts predict that epigenetic mechanisms control transcriptional reprogramming resulting in phenotypic plasticity (Flores et al., 2013). To test this hypothesis, we investigated whether the frequency of H. axyridis colour morphs in response to temperature is mediated by epigenetic mechanisms. To provide an easy readout system, we counted the black spots on the elytra of succinea beetles as previously described (Michie et al., 2010) and assigned them to three groups: low (0-7), middle (8-15) and high (> 16) spot number. The F1eggs of the first filial generation (F1) laid by the F0 generation were kept at three different temperatures (15, 21 or 28 C) and the emerging beetles were tested to determine the frequency of each group. We then used RNA interference (RNAi) to silence genes encoding enzymes responsible for DNA methylation and histone acetylation/deacetylation, two epigenetic mechanisms that are known to regulate the initiation of transcription (Glastad et al., 2011;Vilcinskas, 2017;Gegner et al., 2019).
DNA methylation is mediated by enzymes known as DNA methyltransferases (DNMTs), which catalyse the addition of methyl groups to cytidine residues and favour the formation of compact and inaccessible chromatin, generally blocking the interaction between DNA and transcription factors (Jurkowska et al., 2011). In this process the DNA methyltransferase 1 associated binding protein 1 (DMAP1) is the key activator of DNMT1, as the loss of this protein results in hypomethylation (Lee et al., 2010;Rountree et al., 2010).
In contrast, histone acetyltransferases (HATs) add acetyl groups to the histones, resulting in an open chromatin structure that encourages the binding of transcription factors and promotes gene expression. The activity of HATs is countered by histone deacetylases (HDACs) that remove the acetyl groups and restore the compact and inaccessible chromatin state (Marks et al., 2003).
Using the comprehensive H. axyridis transcriptome sequence (Vilcinskas et al., 2013) and gene predictions on the recently published genome data (Gautier et al., 2018), we selected HAT, HDAC and DMAP1 genes as targets for RNAi-mediated silencing. We found that the phenotypic plasticity of H. axyridis was not affected by the double-stranded RNA (dsRNA) treatment of any of these target genes, but we observed a striking and unexpected impact on life-history parameters, which we investigated in the context of differences in gene expression between invasive, non-invasive and laboratory-bred populations.

Temperature dependency of colour polymorphism
In agreement with the results reported by Michie et al. (2010) we found that beetles raised at lower temperatures from the egg stage produced more spots. For the lowtemperature group (raised at 15 C) all the beetles developed more than 16 black spots (high category). In contrast, in the mid-temperature group (raised at 21 C) two thirds of the beetles belonged to the middle and low spot number categories, and in the high temperature group (raised at 28 C) 92% of the beetles developed fewer than eight spots and were therefore assigned to the low spot number category (Table 1).
These values were used as controls to determine the impact of silencing the epigenetic regulators HAT, HDAC and DMAP1 in succinea colour morphs. Stage four (L4) larvae and F0 beetles were injected with dsRNA corresponding to each gene, but neither the adult beetles developing from the L4 larvae nor the F1 offspring of the injected F0 parents showed any changes in the distribution of spot numbers compared to water-injected mock controls. However, these experiments provided intriguing evidence that RNAi did affect the life-history traits of the F0 and F1 beetles, and we therefore investigated this aspect in greater detail.
Impact of HAT, HDAC and DMAP1 dsRNA treatment on the development, survival and phenotype of L4 larvae and the F1 generation reared at 15 C In L4 larvae, the development time from pupae to adults was shortened following the injection of HAT dsRNA (Fig. 1, Tables S1-S3). Furthermore, significantly fewer animals developed from larvae to pupae and from pupae to adult beetles after RNAi treatment of HAT and HDAC, with HDAC dsRNA having the most severe effect on the larvato-pupa transition. None of the DMAP1 silenced L4 larvae developed into pupae. They either died or remained arrested at the larval stage until the end of the experiment (Tables S1-S4). Following the injection of F0 beetles with HAT or HDAC dsRNA, the F1 generation took longer to develop from eggs to beetles than mock controls, whereas F0 beetles injected with DMAP1 dsRNA failed to produce any offspring at all (Fig. 2, Tables S5-S7). Compared to mock controls, we found that fewer offspring of the F0 beetles injected with HDAC dsRNA developed from larvae into pupae and from pupae into beetles, and the latter was also true for the offspring of beetles injected with HAT dsRNA (Tables S4 and S5). For HAT and HDAC, the injection of dsRNA did not affect the survival of the L4 larvae or the offspring of the F0 beetles. In contrast, the injection of DMAP1 dsRNA resulted in the mortality of 91% of the L4 larvae and, as stated above, there were no offspring from the injected F0 parents (Figs 1 and 2, Tables S1-S7).
Impact of HAT, HDAC and DMAP1 dsRNA treatment on life-history parameters of beetles reared at 21 C, and impact of DMAP1 knockdown on female reproductive organs Having investigated the effect of RNAi treatment of epigenetic regulators on L4 larvae and the F1 generation, we next analysed the impact on the survival, fertility and fecundity of the dsRNA injected beetles. The silencing of DMAP1 was lethal for 100% of the beetles after 24 days, but injection of HDAC dsRNA had only a slight, nonsignificant effect  (Tables S1-S4). DMAP1, DNA methyltransferase 1 associated protein 1; HAT, histone acetyltransferase; HDAC, histone deacetylase; d, days; dpi, days postinjection.
on survival, and HAT-RNAi had no effect at all. The injection of DMAP1 dsRNA caused the first beetles to die after 11 days and the median survival time was 15 days (Fig. 3, Table S8). When mock control females were mated with males injected with DMAP1 dsRNA, there was no difference in the total number of eggs, the hatching rate or offspring development compared to mock control females mated with mock control males, but the number of eggs per batch was lower (Fig. 4, Tables S9 and S10). In contrast, when both parents were injected with DMAP1 dsRNA, or when mock control males were mated with DMAP1 silenced females, the total number of eggs was significantly lower as well as the number of eggs per batch, and none of the eggs produced viable larvae (Fig. 4). The injection of HDAC dsRNA only marginally reduced the total number of eggs and number of eggs per batch compared to controls ( Fig. 4A, B). The hatching rate for the offspring of beetles injected with HDAC or HAT dsRNA was lower than the hatching rate of controls, but there was no difference in the developmental time of the larvae in each cohort (Fig. 4C, D).
Given that female beetles with silenced DMAP1 showed a significant loss of fertility and fecundity, we dissected their ovaries and recorded the number of ovarioles, follicles, mature eggs, and degenerated eggs per ovary. 10 days after injection, the number of ovarioles per ovary was similar in the beetles injected with DMAP1 dsRNA and the controls injected with water. However, the ovaries in the DMAP1 dsRNA treatment group contained significantly fewer follicles and mature eggs, and a much larger number of degenerate eggs, compared to the control group (Fig. 5, Tables S11 and S12). There were no morphological differences between the testes of  (Tables S5-S7). HAT, histone acetyltransferase; HDAC, histone deacetylase; d, days. males injected with DMAP1 dsRNA and those of mock control males (Fig. S4).
To verify that our findings on female fecundity and ovary development result from DMAP1 knockdown, we compared relative gene expression levels of DMAP1 dsRNAinjected beetles with those of mock-injected controls by quantitative real-time PCR (qPCR). And indeed our data showed that DMAP1 gene expression was suppressed about 3.23-fold (AE 1.45 SD) by the dsRNA at 3 days postinjection (dpi) and that the knockdown remained almost stable with 2.87-fold (AE 1.44 SD) at 10 dpi ( Fig. 6A, Tables S13 and S14).

DMAP1 gene expression level in five different H. axyridis populations
Given the importance of DMAP1 in traits such as fecundity, which are relevant for the invasive performance of H. axyridis, we investigated whether DMAP1 gene expression differs between native populations from China and Korea, invasive populations from Germany and France, and a laboratory-reared biocontrol strain with a particularly high reproductive capacity. The qPCR data showed that DMAP1 was expressed at a lower level in the native and invasive populations compared to the biocontrol strain, but there was no difference between the native and invasive populations (Fig. 6B, Tables S15 and S16).

Discussion
A very recent study combined whole-genome sequencing, population-based genome-wide association studies, gene expression, and functional analyses to determine that the transcription factor Pannier regulates melanic pattern polymorphism in H. axyridis and that highly variable discrete colour forms in natural populations result from cisregulatory allelic variation of a single gene (Gautier et al., 2018). In contrast, polyphenism in insects involves the environmental modulation of phenotypes that may be specified by an invariant genome or may already exist as a number of different allelic variants, thus contributing to the penetrance and expressivity of different alleles. The mechanism that underlies polyphenism is unclear, but epigenetic mechanisms may explain some of the features of this phenomenon, such as the ability of environmental pressure in one generation to affect the phenotype of later generations (Glastad et al., 2019;Simpson et al., 2011). In the invasive harlequin ladybird H. axyridis, we confirmed the already reported influence of temperature during the development of colour morphs (Michie et al., 2010). The exposure of eggs to a lower temperature increased the number and extent of black spots on the elytra of the resulting adults. Testing our hypothesis that epigenetic mechanisms acting at the level of transcriptional initiation could translate environmental stimuli into the frequencies of different morphs, we injected dsRNA corresponding to the genes encoding key regulators of DNA methylation (DMAP1) and histone modification (HATs and HDACs) into L4 larvae and F0 adults. However, no effects on the frequency of different morphs in the adults developing from injected larvae or the F1 offspring of the injected adults were observed, indicating that neither DNA methylation nor histone acetylation is solely responsible for the polyphenism observed in H. axyridis.
Although our initial hypothesis was rejected, we observed that RNAi treatment of HAT and HDAC genes influenced H. axyridis life-history parameters including survival and development. Injection of HAT dsRNA did not affect the survival of larvae under cold stress or beetles under normal conditions, but it did delay the developmental transition from larvae to pupae under cold stress compared to controls. In Tribolium castaneum, the inhibition of HAT by curcumin did not influence larval survival or development, but extended the life span of beetles under heat stress (Bingsohn et al., 2016). Similarly, the longevity of adult Drosophila melanogaster was increased after treating larvae with the same HAT inhibitor under normal conditions (Soh et al., 2013). However, the inhibition of HAT with curcumin did not affect D. melanogaster reproduction (Soh et al., 2013) but did reduce the fertility of T. castaneum (Bingsohn et al., 2016). In our experiments, we found that RNAi treatment of HAT reduced the fertility of H. axyridis but had no impact on its fecundity.
Interestingly, injection of HDAC dsRNA caused similar effects as shown for HAT, ie reduced fertility but no changes in fecundity or survival. However, one key . Kaplan-Meyer survival curves of Harmonia axyridis beetles injected with histone acetyltransferase (HAT)/histone deacetylase (HDAC)/ DNA methyltransferase 1 associated protein 1 (DMAP1) double-stranded RNA (dsRNA) or injected with water as a mock control. The survival curve of beetles injected with DMAP1 dsRNA differs significantly from those injected with HAT/HDAC dsRNA or with water as a mock control (P < 0.001), but no significant differences were observed for other treatments (mock vs. HAT: P = 1.000, mock vs. HDAC = HAT vs. HDAC: P = 0.120) (Table S8). d, days.
difference was that RNAi treatment of HDAC did not extend the transition times of any stages of development. The HDAC inhibitor valproic acid reduced the survival of T. castaneum larvae and heat-stressed beetles but had no effect on the longevity of D. melanogaster even though it reduced the fertility and fecundity of both species (Bingsohn et al., 2016;Gayathri and Harini, 2012). Interestingly, earlier studies concluded that the inhibition of HDAC by sodium 4-phenylbutyrate or trichostatin A and sodium butyrate increased the longevity of D. melanogaster but did not affect reproduction (Kang et al., 2002;Zhao et al., 2005). These contradicting findings highlight the need for more studies in insects to investigate the complex regulatory network of histone acetylation and deacetylation.
DNA methylation is also intimately linked to the regulation of genes involved in life-history traits, and a key component of the machinery is DMAP1 (Lee et al., 2010;Rountree et al., 2010). The RNAi-mediated silencing of DMAP1 resulted in severe effects, killing 50% of the treated beetles after 15 days and 100% after 24 days. To the best of our knowledge, this is the first report showing the lethal consequence of silencing the DNA methylation machinery in adult insects, although lethality has been observed in insect and mammalian embryos (Flores et al., 2013). In the mammalian embryos, demethylated DNA induced apoptosis, and we hypothesize that the same process is responsible for the lethality observed in H. axyridis.
Interestingly, we found that the reduced fertility and fecundity caused by DMAP1 silencing was only evident in females, corresponding to the presence of severely degenerated ovaries. To our knowledge this is the first time that  (Tables S9 and  S10). DMAP1(f), only females injected, mated with mock-injected males; DMAP1(m), only males injected, mated with mock-injected females; d, days. the essential role of the DNMT1 cofactor DMAP1 in insect reproduction has been shown. Until now, similar results were only reported for DNMT1 and DNMT3. For example, for the brown plant hopper, Nilaparvata lugens, it was recently reported that the knockdown of DNMT1 and DNMT3 reduced the number of offspring by suppressing ovary development (Zhang et al., 2015). Further, RNAi-mediated silencing was also used to demonstrate that DNMT1 is essential for egg production and embryo viability in the large milkweed bug, Oncopeltus fasciatus (Bewick et al., 2019). A recent study using T. castaneum elucidated that DNMT1 is expressed throughout the entire life cycle of this model beetle, but the highest levels of messenger RNA (mRNA) transcripts were found in eggs and Figure 5. Ovaries of Harmonia axyridis females 10 days after the injection of (A) mock or (B) DNA methyltransferase 1 associated protein 1 (DMAP1) doublestranded RNA, as well as the mean numbers of (C) ovarioles, (D) follicles, (E) mature eggs and (F) degenerated eggs per ovary: n(mock) = 12, n(DMAP1) = 11. Error bars show standard error of the mean and significant differences between groups are depicted by different letters (Tables S11 and S12). [Colour figure can be viewed at wileyonlinelibrary.com].
ovaries. Further, maternal silencing of this gene caused a developmental arrest in offspring embryos (Schulz et al., 2018). Taking together these studies and our findings, we postulate an evolutionarily conserved role of the DNMT1-DMAP1 complex beyond embryonic development including a role in ovary maturation and female fecundity.
In our study, DMAP1 silencing did not affect male fecundity or the integrity of the testes. However, we observed a slight decrease in the number of eggs per batch when mock control females (injected with water) were mated with males injected with DMAP1 dsRNA, suggesting that the DNMT1-DMAP1 complex plays a minor role in the development of male reproductive organs. This is supported by the fact that DNMT1 is strongly expressed not only in the ovaries of N. lugens, but also in the testes of the fire ant Solenopsis invicta (Kay et al., 2018;Lu et al., 2018).
In order to gain further insight into the role of the DNA methylation machinery in the fecundity of H. axyridis, we compared the expression of DMAP1 in two invasive populations, two native populations and a biocontrol strain, which reportedly differ in their reproductive capacity. The latter was specifically bred to achieve a higher daily reproduction rate for pest control applications (Tayeh et al., 2012;Tayeh et al., 2015). We found that DMAP1 mRNA was 27% more abundant in the biocontrol strain compared to native populations from Korea and China and 38% more abundant compared to invasive populations from France and Germany. These findings suggest that the DNMT1-DMAP1 complex helps to determine the fecundity and fertility of this invasive species and may therefore partly contribute to its invasive success.
To summarize, we found that the RNAi treatment of epigenetic regulators can affect the life-history traits of H. axyridis. The knockdown of DMAP1 in particular led to severe effects on survival and fecundity. We also showed that DMAP1 was expressed at much higher levels in a biocontrol strain bred for increased fecundity compared to native and invasive populations. Further research is needed to determine how epigenetic mechanisms, and in particular DNA methylation, promote the success of this highly invasive species.

Rearing and collection of H. axyridis
H. axyridis populations from Germany, France, Korea and China, as well as a laboratory-bred biocontrol strain, were collected and reared as previously described (Gegner et al., 2018). For mating, adult animals were paired in Petri dishes (94 × 16 mm) and fed ad libitum with pea aphids (Acyrthosiphon pisum). All experiments were carried out with beetles 2-3 weeks of age representing the colour morph succinea.

Analysis of colour polymorphism at different temperatures
The effect of temperature on colour polymorphism was analysed by separating F0 beetles into three breeding groups according to number of black spots on the elytra, namely low (0-7), middle (8-15) and high (> 16) spot number, as previously described (Michie et al., 2010). F1 eggs were maintained at 15, 21 or 28 C. The emerging F1 larvae were reared individually and the colour morphs and spot groups of the resulting F1 beetles were recorded separately for each breeding group.

Identification of epigenetic regulators and preparation of dsRNA
To identify putative H. axyridis HATs, HDACs and DMAP1, the corresponding amino acid sequences from Leptinotarsa decemlineata [National Center for Biotechnology Information (NCBI) accession numbers XP_023029651, XP_023012143 and Error bars indicate standard error of the mean. Significant gene expression knockdown is depicted by asterisks. Significance levels: P < 0.05 (*), P < 0.01 (**) (Tables S13 and S14). (B) Comparison of the DMAP1 gene expression in untreated Harmonia axyridis beetles of two invasive and two native populations relative to untreated individuals of a biocontrol strain. DMAP1 gene expression of the biocontrol strain is represented as baseline. Dots show mean values of 10 biological replicates (each comprising three males and three females). Error bars indicate standard error of the mean. Significant differences from biocontrol strain expressions are depicted by asterisks. Significance level: P < 0.001 (***) (Tables S15 and S16). XP_023020626] and T. castaneum (NCBI accession numbers XP_966673, XP_967425 and XP_008195930) were used as TBLASTN queries with default parameter settings against our inhouse H. axyridis transcriptome assemblies (Vilcinskas et al., 2013). Sequences with more than 55% amino acid sequence similarity to queries were collected for further analysis. All protein sequences were aligned in GENEIOUS vR11 (Biomatters Ltd, Auckland, New Zealand) using MUSCLE with default settings. The results were inspected for regions of high-quality alignment and were refined manually. During this step, candidates were also scrutinized for the presence of conserved amino acid patterns. To find possible gene homologues for the final candidates of HAT, HDAC and DMAP1, a BLASTN v. 2.6.0 search (Altschul et al., 1997) was performed against the H. axyridis genome HaxR_v1.0 (Gautier et al., 2018) with default settings and a word size of seven. We could not find any homologues for HAT and DMAP1 in the H. axyridis genome, as HAT mapped only to utg 63 pilon from position 2308311 to 2310853 and DMAP1 only to utg_158_pilon from position 488535 to 490820 (File S1). For HDAC we found two highly similar homologues, HDAC 4 and 7. One copy mapped completely to utg_635_pilon from position 102105 to 141466. The other copy mapped to utg_2690_pilon from position 160556 to 169854 and to utg_2046_pilon from position 273426 to 274154 (Files S1 and S2). For each possible homologous sequence a BLASTX (BLASTN v. 2.9.0) search was performed against the non-redundant protein sequences (nr) database using default settings. To exclude off-target effects, the sequence of dsRNA belonging to H. axyridis was split into all possible 21mers and a BLASTN search was performed against the H. axyridis genome HaxR_v1.0 as described above. Hits were only defined as potential off-targets if they had mismatches at position 1 and/or 21 and showed 100% identity in the remaining 19mer. In order to define the position of the potential off-targets 10 kb sequences upstream and downstream of these regions were extracted from the genome. Ab initio gene prediction was performed using the AUGUSTUS pipeline implemented in the OMICSBOX software (BioBam Bioinformatics S.L., Valencia, Spain) and the FGENESH HMM-based gene structure prediction and gene-finder tool with standard settings (Salamov and Solovyev, 2000). No off-targets could be found for our HAT or DMAP1 dsRNA constructs (File S2, Figs S1 and S2). The possible small-interfering RNAs (siR-NAs) resulting from the HDAC dsRNA (Fig. S3) were determined to address both on-targets described above, as well as two very short noncoding sections of the genome (utg_331_pilon from position 1696692 to 1696710 and utg_893_pilon from position 3749566 to 3749584; Files S2-S4). We can thus exclude that the selected siRNAs could potentially affect other, nontarget genes in Harmonia axyridis. Utilizing this approach, we selected mRNA sequences for HAT, HDAC and DMAP1 that were then used as templates for the preparation of dsRNAs (Table S17) using an Ambion MEGAscript T7 kit (Applied Biosystems, Foster City, CA, USA).

Monitoring life parameters after the RNAi treatment of epigenetic regulators
We injected 1 μg dsRNA into adult beetles as previously described (Gegner et al., 2018) or into L4 larvae laterally between the sixth and seventh tergites. The life parameters and colour phenotypes of 33 larvae per gene and 35 mock-larvae were monitored for 2 months. To investigate the impact of parental RNAi in offspring beetles, we injected one pair of F0 beetles from the low-spotnumber group (development temperature 21 C) three times with HAT, HDAC or DMAP1 dsRNA, or with water as a mock-control. The F1 generation was then reared at 15 C and the colour phenotype and life-history parameters were recorded. The experiment was terminated 3 months after injection. The survival of F0 beetles was monitored until 35 dpi. Fecundity (number of egg batches and number of eggs per batch) was recorded until the last DMAP1 dsRNA-injected female beetle died (24 dpi). We also determined the hatching rate and the developmental time of F1 larvae.

Influence of DMAP1 knockdown on reproductive organs
The influence of DMAP1 knockdown on reproductive organs was investigated 10 dpi by dissecting testes and ovaries in phosphate-buffered saline for analysis by light microscopy. In females, we also recorded the numbers of ovarioles, follicles, mature eggs and degenerate eggs.

Gene expression analysis
RNA isolation, cDNA synthesis and qPCR for relative gene expression analysis were carried out as previously described, using the ribosomal protein S3 gene (RPS3) for normalization (forward primer 5'-GGCTACCAGAACCGACAGAG-3 0 and reverse primer 5'-GTGCTATGGCGCATAATCCT-3 0 ; Gegner et al., 2018). DMAP1 gene knockdown mediated by dsRNA was verified with a total of three biological replicates per treatment and time point, each consisting of 10 pooled individuals with equal sex ratio, in a qPCR assay using forward primer 5'-TACCTGCAAATGTGG GTC-3 0 and reverse primer 5'-GACCATGTCGCTTCTAAGTTCG-3 0 . For each replicate 20 beetles (10 males and 10 females) were injected with either DMAP1 dsRNA or water as control. After 3 days five males and five females per treatment were flash frozen in liquid nitrogen and pooled for RNA isolation. 10 days after injection the procedure was repeated for the remaining beetles. Fold changes in gene expression relative to controls (2 ΔΔCт ) (Pfaffl, 2001;Schmittgen and Livak, 2008) were calculated and the 3 and 10 dpi means were tested for statistical significance.
Population-specific differences in DMAP1 gene expression were determined by qPCR amongst five untreated biological replicates representing two invasive populations from Germany and France, two native populations from China and Korea, and one biocontrol strain. Six beetles (three males and three females) for each treatment (RNAi and water controls) were pooled 10 dpi and the experiment was performed twice. For statistical analysis, the results for each invasive population (Germany and France) and each native population (China and Korea) were calculated and pooled. The population means were then tested against each other and against the biocontrol strain.

Statistical analysis
Statistical analysis was carried out using R v. 3.3.3 (R Core Team, 2015). The dunnTest function in package FSA v. 0.8.17 (Ogle, 2016) was used for Kruskal-Wallis multiple comparisons with Bonferroni adjusted P-values to analyse the development and survival rate of the F1 progeny of dsRNA-injected F0 beetles and of dsRNA-injected L4 larvae, as well as the development, fecundity and fertility of dsRNA-injected beetles, and the ovary parameters of DMAP1-RNAi females. The survival rate and the numbers of animals per developmental stage of the F1 progeny of dsRNA-injected F0 beetles and dsRNA-injected L4 larvae that developed at 15 C were analysed by pairwise comparisons using a two-sided Fisher's exact test with Bonferroni correction provided by the functions fisher.test and fisher. multcomp in package RVAideMemoire v. 0.9-68 (Hervé, 2017). The survival of F0 beetles reared at 21 C was investigated using the R package survival v. 2.40-1 (Therneau, 2015) to plot Kaplan-Meier survival curves and to calculate differences in survival rates between treatments with the log-rank test and Holm-corrected P-values. Statistical analysis of qPCR data was carried out according to Gegner et al. (2018), based on the ΔΔC T method established by Pfaffl (2001) and Schmittgen and Livak (2008).

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.
Table S13. Number of replicates (n), mean, standard deviation (sd), and standard error of the mean (se) for fold changes in gene expression of DNA methyltransferase 1-associated protein 1 after knockdown relative to mock control at 3 and 10 days postinjection. Table S14. Results of the multiple comparison of means analysis, utilizing a one-way analysis of variance and simultaneous tests for general linear hypotheses, for DNA methyltransferase 1-associated protein 1 (DMAP1) gene expression of DMAP1-injected beetles and mock controls 3 and 10 days postinjection. Significance levels: P < 0.05 (*), P < 0.01 (**), P < 0.001 (***).
Table S15. Number of replicates (n), mean, standard deviation (sd) and standard error of the mean (se) for fold change in gene expression of DNA methyltransferase 1 associated protein 1 relative to biocontrol strain. Table S16. Results of the multiple comparison of means analysis, utilizing a one-way analysis of variance and simultaneous tests for general linear hypotheses, for DNA methyltransferase 1-associated protein 1 gene expression of native, invasive and biocontrol strains. Significance levels: P < 0.05 (*), P < 0.01 (**), P < 0.001 (***). Table S17. Gene-specific primers without the T7 polymerase promoter sequence at the 5 0 end (5 0 -TAA TAC GAC TCA CTA TAG GGA G-3 0 ) generating double-stranded RNA of epigenetic regulators and the length of resulting PCR products: histone acetyltransferase, histone deacetylase and DNA methyltransferase 1-associated protein 1. Figure S1. Sense strand sequence (5 0 to 3 0 ) of the histone acetyltransferase messenger RNA, including forward (dark green) and reverse (light green) primer for double-stranded RNA (dark blue) synthesis. Figure S2. Sense strand sequence (5 0 to 3 0 ) of the DNA methyltransferase 1-associated protein 1 messenger RNA, including forward (dark green) and reverse (light green) primer for double-stranded RNA (dark blue) synthesis. Figure S3. Sense strand sequence (5 0 to 3 0 ) of the histone deacetylase messenger RNA, including forward (dark green) and reverse (light green) primer for double-stranded RNA (dark blue) synthesis. Figure S4. Testicles of (A) mock-and (B) DNA methyltransferase 1-associated protein 1-injected Harmonia axyridis males at 10 days postinjection.