Genome-wide characterization of miRNA and siRNA pathways in the parasitoid wasp Pteromalus puparum

XIAO Shan,FANG Qi,LIU Ming-ming,ZHANG Jiao,WANG Bei-bei,YAN Zhi-chao,WANG Fang,David W.STANLEY,YE Gong-yin

1 State Key Laboratory of Rice Biology/Ministry of Agriculture and Rural Affairs Key Laboratory of Molecular Biology of Crop Pathogens and Insects,Institute of Insect Sciences,Zhejiang University,Hangzhou 310058,P.R.China

2 Biological Control of Insects Research Laboratory,USDA Agricultural Research Service,Columbia,Missouri 65203,USA

Abstract microRNAs (miRNAs) and small interfering RNAs (siRNAs) are small non-coding RNAs (ncRNAs) that trigger RNA interference(RNAi) in eukaryotic organisms.The biogenesis pathways for these ncRNAs are well established in Drosophila melanogaster,Aedes aegypti,Bombyx mori and other insects,but lacking in hymenopteran species,particularly in parasitoid wasps.Pteromalus puparum is a parasitoid of pupal butterflies.This study identified and analyzed two pathways by interrogating the P.puparum genome.All core genes of the two pathways are present in the genome as a single copy,except for two genes in the siRNA pathway,R2D2 (two copies) and Argonaute-2 (three).Conserved domain analyses showed the protein structures in P.puparum were similar to cognate proteins in other insect species.Phylogenetic analyses of hymenopteran Dicer and Argonaute genes suggested that the siRNA pathway-related genes evolved faster than those in the miRNA pathway.The study found a decelerated evolution rate of P.puparum Dicer-2 with respect to Dicer-1,which was contrary to other hymenopterans.Expression analyses revealed high mRNA levels for all miRNA pathway genes in P.puparum adults and the siRNA related genes were expressed in different patterns.The findings add valuable new knowledge of the miRNA and siRNA pathways and their regulatory actions in parasitoid wasps.

Keywords:miRNA pathway,siRNA pathway,annotation,Pteromalus puparum

1.Introduction

MicroRNAs (miRNAs) and small interfering RNAs (siRNAs)are non-coding RNAs (ncRNAs) that modulate gene expressionviainteractions with mRNA.The functions of these ncRNAs have been reported to be implicated in a wide range of biological processes (Carthew and Sontheimer 2009).Both ncRNAs are readily synthesized and are promising tools for a wide variety of medical and agricultural applications (Davuluriet al.2005;Baumet al.2007;Mendell 2008;Costa and Torchilin 2018).In insects,most genes encoding miRNA are transcribed by RNA polymerase II and generate long primary miRNA (pri-miRNA) transcripts(Leeet al.2003,2004).The resulting pri-miRNAs are bound and excised into approximately 70 nucleotide miRNA precursors (pre-miRNAs) by a protein complex composed of a ribonuclease III (RNase III) called Drosha and its cofactor Pasha,a dsRNA-binding protein (dsRBP) that is an ortholog to DGCR8 in vertebrates (Denliet al.2004).The pre-miRNAs are exported from the nucleus by Exportin-5 (Yiet al.2003).In the cytoplasm,the dsRBP Loquacious interacts with another RNase III,Dicer-1,to facilitate processing of pre-miRNA into a miRNA duplex (Kettinget al.2001;Jianget al.2005).The protein Argonaute-1 (Ago-1) interacts with one strand of the duplex to form the active miRNA-induced silencing complex (miRISC) (Forstemannet al.2007).After imperfect pairing between miRNAs and mRNAs,miRISCs induce mRNA destabilization or translational inhibition,which effectively silences target genes (Leeet al.1993;Rehwinkelet al.2005;Kiriakidouet al.2007).

The siRNA pathway is launched by dsRNAs derived from either exogenous virus or endogenous overlapping transcripts,structured loci and transposons (Ghildiyalet al.2008;Okamuraet al.2008).InDrosophilamelanogaster,the RNase III,Dicer-2,operates with a different Loquacious protein isoform to cleave long dsRNAs into duplex siRNAs(Bernsteinet al.2001;Marqueset al.2010).After cleavage,the protein R2D2,a dsRBP which also interacts with Dicer-2,subsequently recruits duplex siRNA to Argonaute-2 (Ago-2)(Liuet al.2003;Marqueset al.2010).Ago-2 interacts with one selected strand and forms a siRNA-induced silencing complex(siRISC),which induces RNA silencing by perfectly binding and cleaving the target mRNA (Matrangaet al.2005).

Studies on miRNAs and siRNAs have implicated both RNAi pathways in insect host-pathogen interactions and immunity (Hussain and Asgari 2014;Leggewie and Schnettler 2018).miRNAs are also associated with major biological processes such as development,metamorphosis and reproduction (Lucas and Raikhel 2013).The formation and biological significance of ncRNAs have been established in detail inD.melanogasterand a few other insect species,but not in hymenopterans,particularly not in parasitoid wasps.Such information on parasitoids may be valuable because these very small insects are applied in biological control programs globally to manage serious crop pest insects.Understanding and applying ncRNA mechanisms may lead to considerable enhancements in biological control programs.Pteromaluspuparum(Hymenoptera:Pteromalidae) is a parasitoid of some pupal butterflies such as the global cabbage white butterflyPierisrapae(Caiet al.2004).As seen in a related parasitoid,Nasoniavitripennis,the advantages of short life history,high fecundity,and diverse omics data makeP.puparuma model of choice for investigations into the co-evolution of host/endoparasitoid systems and into the possibilities of applying tools of modern biotechnology for the improvement of biological control operations (Werren and Loehlin 2009).Drawing on a completeP.puparumgenome,the study launched investigations into parasitoid ncRNA biology by posing the hypothesis that genes necessary for formation and operation of the ncRNA pathways are present inP.puparum.

2.Materials and methods

2.1.Insect rearing

Laboratory cultures ofP.puparumwere fed on a 10%sugar solution in distilled water,and reared at 25°C under a 14 h L:10 h D photoperiod as described by Caiet al.(2004).The hostP.rapaelarvae were collected from experimental farmlands around Zhejiang University (Hangzhou,China) and fed on cabbage leaves in the laboratory until pupation.Each mated female wasp parasitized aP.rapaepupa for 24 h.

2.2.Annotation of genes operating in the miRNA and siRNA pathways

This study identifiedP.puparummiRNA and siRNA pathway protein sequences by obtaining appropriate sequences inD.melanogaster,TriboliumcastaneumandApismelliferafrom NCBI and used them as query sequences for BLASTP(e-value:1e-5) searching within theP.puparumofficial protein set (http://www.insect-genome.com/Pteromalus/).To identify cognate genes,the reference protein sequences together with the resultingP.puparumproteins were subjected to TBLASTN (e-value:1e-5) analysis within theP.puparumgenome.Resulting new loci were subjected to FGENESH (http://linux1.softberry.com/) for gene prediction and the predicted proteins were inspected by BLASTP analysis in the NCBI NR database (Solovyevet al.2006).The conserved domains of all identified proteins were analyzed against SMART (version 7.1) and Pfam (version 32.0) databases with InterProScan,and illustrated using IBS (version 1.0) (Joneset al.2014;Liuet al.2015;Letunic and Bork 2018).Transmembrane helices of Sid1 proteins were predicted with TMHMM (version 2.0) (Kroghet al.2001).Similar procedures were performed for other insect species.The species and protein identifications and newly identifiedprotein sequences included in this study are listed in Appendices A and B,respectively.

2.3.Phylogenetic analysis

All resulting protein sequences were first aligned with MAFFT (version 7.123b) and then trimmed by trimAl(version 1.4.rev22) (Capella-Gutierrezet al.2009;Katoh and Standley 2013).The best substitution model was determined by the ModelFinder implemented in IQ-TREE(version 1.6.7) according to Bayesian Information Criterion(BIC) (Nguyenet al.2015;Kalyaanamoorthyet al.2017).Models suggested by ModelFinder are:LG+F+R4 model for Drosha and Dicer;LG+R4 for Pasha;LG+F+G4 for Exportin-5 and Argonaute;JTTDCMut+G4 for Loquacious;LG+R3 for R2D2;and LG+I+G4 for Sid1.Phylogenetic trees were constructed by IQ-TREE with 1 000 bootstrap replicates.TheTetranychusurticae,PenaeusvannameiorCaenorhabditiselegansortholog was applied as the outgroup for all phylogenetic trees,with the exception of Argonaute proteins,in which Piwi sub-family proteins ofD.melanogasterandP.puparumwere used.

To compare the molecular evolution rates of related genes in the two ncRNA pathways,multiple sequence alignments at the nucleotide level were generated using PAL2NAL with amino acids alignments as guidance(Suyamaet al.2006).The free ratio model was applied to calculate the nonsynonymous substitutions (dN),synonymous substitutions (dS),and their ratio (ω) for each gene using CODEML implemented in PAML (version 4.7)(Yang 2007).The phylogenetic trees were used as input trees for the CODEML analyses.Branches with dN or dS＞3,suggesting saturation of substitutions,were not included in the analyses.

2.4.RNA isolation and qPCR

Pteromaluspuparumsamples were collected across developmental stages with three replicates:2nd-,4th-and 6th-day larvae,white pupae,yellow pupae,half black pupae,1st-and 3rd-day adults.qPCR analysis was conducted according to Quantitative Real-Time PCR Experiment (MIQE) guidelines (Bustinet al.2009).In brief,the samples were processed in TRIzol?Reagent(Invitrogen,Carlsbad,CA,USA) for total RNA extraction according to the manufacturer’s instructions.TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix Kit(TransGen Biotech,Beijing,China) were used for reverse transcriptions.qPCR was performed using the Bio-Rad CFX 96 Real-Time Detection System (Bio-Rad,Hercules,CA,USA) with ChamQTMSYBR qPCR Master Mix Kit (Vazyme,Nanjing,China).The thermocycler was programmed at 95°C for 3 min,40 cycles of 95°C for 10 s and 60°C for 30 s.To verify the specificity of the amplification,a dissociation curve was included from 60-95°C at the end of each qPCR run.The stably expressedP.puparum18srRNA gene was selected as reference gene,and the quantitative variation for each gene was calculated using a relative quantitative method (2-ΔΔCT) (Livak and Schmittgen 2001;Wanget al.2013).Primers used in this study are listed in Appendix C.

2.5.Statistical analyses

The comparisons of dN,dS,and ω values between paralog genes were evaluated using Student’st-test (P＜0.05).To verify the normality of the data,variable homogeneity test was conducted before Student’st-test.All statistical analyses were performed using the Data Processing System software (version 14.50) (Tang and Zhang 2013).

3.Results

3.1.Structural and phylogenetic similarity of P.puparum miRNA pathway genes

As previously reported for most insect species,P.puparumhas one copy of each miRNA gene includingDrosha,Pasha,Exportin-5,Loquacious,Dicer-1andAgo-1(Table 1).Conserved domain analysis exhibits similar structures to proteins in other insect species.Pteromalus puparumRNases III Drosha and Dicer-1 contain two tandem RNase III domains,necessary for cleaving pri-and pre-miRNAs(Fig.1-A).Drosha includes a C-terminal dsRNA binding domain (dsRBD);Dicer-1 has three other domains:an RNA helicase domain,a Dicer dimer domain and a Piwi-Argonaute-Zwille (PAZ) domain.Pteromalus puparumPasha has two dsRBDs.Loquacious proteins include three dsRBDs in all six species.Exportin-5 has an Importin-beta N-terminal domain (IBN_N) and an Exportin 1-like protein domain (Xpo1).Ago-1 contains six domains including PAZ and Piwi (P element-induced wimpy testis) which are shared by other Argonaute proteins.

To confirm that theP.puparumproteins were correctly annotated,this study constructed phylogenetic trees for the proteins along with those of other insect species.Based on amino acid sequences,the phylogenetic analysis placedP.puparumandN.vitripennisproteins in the same clade with high bootstrap values in all trees except Ago-1,in whichCeratosolensolmsiwas closer toN.vitripennisthanP.puparum(Appendix D;Figs.2 and 3).All hymenopteran proteins together formed a distinct clade in the trees.

3.2.Gene expansion and variable molecular evolution rate of siRNA pathway

One copy ofDicer-2andSid1were identified in theP.puparumgenome and expansions ofAgo-2andR2D2were observed with three and two copies,respectively(Table 1).Protein domain comparisons showed thatP.puparumDicer-2 included two more domains than Dicer-1:DEAD-like helicase domain (DEAD) and dsRBD (Fig.1-B).Compared to Ago-1,Ago-2 proteins lack a middle domain(ArgoMid).Two dsRBDs at the N-terminal appear in all R2D2 proteins.For protein Sid1,P.puparumSid1 has 12 transmembrane helices (TM),one more thanN.vitripennisandA.mellifera.

Fig.1 Schematic depiction of domain organization of miRNA (A) and siRNA (B) pathway proteins in six insect species.Ppup,Pteromalus puparum;Nvit,Nasonia vitripennis;Amel,Apis mellifera;Tcas,Tribolium castaneum;Dmel,Drosophila melanogaster;Bmor,Bombyx mori;RNase III,ribonuclease III domain;dsRBD,double-stranded RNA-binding domain;Helicase,helicase conserved C-terminal domain;Dicer dimer,Dicer dimerisation domain;PAZ,Piwi-Argonaute-Zwille domain;IBN_N,importin-beta N-terminal domain;Xpo1,exportin1-like protein domain;ArgoN,N-terminal domain of argonaute proteins;ArgoL1,argonaute linker 1 domain;ArgoL2,argonaute linker 2 domain;ArgoMid,Mid domain of argonaute proteins;Piwi,Piwi domain;DEAD,DEAD-like helicase domain;TM,transmembrane helices.

Table 1 Description of the miRNA and siRNA pathway genes identified in Pteromalus puparum

Fig.2 Maximum likelihood tree of multiple sequence alignment of Dicer proteins.Tree rooted at Caenorhabditis elegans.Pteromalus puparum proteins are pointed out by pentagram.Nodes with colors represent different bootstrap support values.

TheP.puparumAgo-2genes were tandemly arranged in the genome with more than 74% identities at the protein level,indicating tandem duplications ofAgo-2(Appendix E).The expansions were also seen inN.vitripennis,Copidosoma floridanumandT.castaneum(Fig.3).A similar arrangement ofAgo-2genes was present inN.vitripennis,but not in other two insect species,indicating different gene duplication mechanisms.TwoR2D2genes were separated in different scaffolds and alignment of the proteins indicated ＜30%identities among insects,even within a species (Appendix E).The study identified twoR2D2genes in five other insect species includingN.vitripennis,C.floridanum,Solenopsis invicta,T.castaneumandAgrilusplanipennis(Fig.4).Phylogenetic analysis showed that theP.puparumR2D2-2occurred in a clade with aN.vitripennishomolog and were placed as the sister group to several hymenopteran R2D2 genes includingP.puparumR2D2-1.

Fig.3 Maximum likelihood tree of multiple sequence alignment of Argonaute-1 (Ago-1),Argonaute-2 (Ago-2) and Piwi sub-family proteins.Tree rooted at Piwi subgroup proteins.Pteromalus puparum proteins are pointed out by pentagram.Nodes with colors represent different bootstrap support values.

Fig.4 Maximum likelihood tree of multiple sequence alignment of R2D2 proteins.Tree rooted at Tetranychus urticae.Pteromalus puparum proteins are pointed out by pentagram.Nodes with colors represent different bootstrap support values.

Phylogenetic trees revealed substantially longer branches ofDicer-2andAgo-2than their paralogs involved in miRNA pathway,which inferred higher protein evolution rates in siRNA pathway genes than in miRNA pathway genes.To test whether there are significant differences among hymenopteran insects,the study calculated dN,dS and their ratio ω for each species using CODEML.It recorded greater ω values forDicer-2(mean=0.231,n=6) compared toDicer-1(mean=0.050,n=6) (Fig.5-A).Likewise,the dN and ω values forAgo-2(mean=0.078 and 0.187,n=12) genes were significantly higher compared toAgo-1(mean=0.003 and 0.006,n=7).The dS values forAgo-1(mean=0.511,n=7) andAgo-2(mean=0.443,n=12)genes were not different (Fig.5-B).Among theP.puparumgenes,the ω values forAgo-2genes ranged from 0.151 to 0.302 and forAgo-1it is 0.012 (Appendix F).ForDicergenes,the ω value forDicer-2(0.059) was lower compared toDicer-1(0.091),contrary to other hymenopterans noted in this study.

Fig.5 Comparison of nonsynonymous substitutions (dN),synonymous substitutions (dS),and their ratio (ω) values between hymenopteran Dicer (A) and Argonaute (B) genes.For Dicer,six Dicer-1 and Dicer-2 genes respectively identified from Pteromalus puparum,Nasonia vitripennis,Ceratosolen solmsi,Apis mellifera,Bombus terrestris and Solenopsis invicta were used for comparison.For Argonaute,7 Ago-1 and 12 Ago-2 genes identified from P.puparum,N.vitripennis,C.solmsi,C.floridanum,A.mellifera,B.terrestris and S.invicta were used.The histogram bars indicate mean relative accumulation of the values and the error bars indicate SE (n=3).Significance of the comparisons were determined by Student’s t-test.The original estimated parameters are shown in Appendix F.

3.3.Expression patterns of P.puparum miRNA and siRNA pathway genes

We investigated gene expression profiles related to miRNA and siRNA pathways across developmental stages.mRNA accumulations were stable throughout larval and pupal stages with the exception ofAgo-1,which had a peak at the 6th-day larvae (L6) (Fig.6-A).Higher accumulations of mRNAs encoding all miRNA proteins occurred in adults.The expression patterns differed among siRNA pathway genes.Two expression peaks ofDicer-2appeared at L6 and the first day post eclosion (Fig.6-B).Compared to other periods,expressions ofR2D2-1andR2D2-2was higher in the 2nd-day larvae and adults.Higher expression ofAgo-2aappeared in yellow pupae and adults.Ago-2cpeaked at L6 andAgo-2bwas relatively stable throughout development.Accumulations of mRNAs encoding Sid1 peaked at the third day post eclosion.

Fig.6 Expression profiles of miRNA (A) and siRNA (B) pathway genes across development stages.L2,L4 and L6 indicate 2ndday,4th-day and 6th-day larvae,respectively;WP,white pupae;YP,yellow pupae;HBP,half black pupae;1 DPE and 3 DPE indicate 1st-day and 3rd-day post eclosion.Data represent the mean±SE (n=3).

4.Discussion

The data presented in this paper strongly support our hypothesis that genes necessary for formation and operation of the ncRNA pathways are present inP.puparum.Substantial datasets supported the argument.First,the study identified genes encoding all proteins operating in the miRNA and siRNA pathways.Second,the structures of each of the proteins were consistent with known structures of analogous proteins in five other insect species.Third,extensive phylogenetic analysis placedP.puparumproteins in a clade withN.vitripennisand documented substantial similarities between theP.puparumproteins and such proteins from other insect species.Fourth,the expressions of genes were recorded across development stages,documenting the expression of the genes throughout the life cycle.Taken together,these points amounted to a solid argument supporting the hypothesis.

In general,miRNAs regulate expression of genes associated with a wide range of biological processes such as development,immunity and behavior (Lucas and Raikhel 2013;Belles 2017).The results determined that all genes in the miRNA pathway were present inP.puparumgenome with one copy,consistent with most insects (Dowlinget al.2016;Mongelli and Saleh 2016).Expansions of miRNA pathway genes have been reported only in the house flyMuscadomesticaand several aphid species in the Aphidinae (Jaubert-Possamaiet al.2010;Ortiz-Rivaset al.2012;Mongelli and Saleh 2016).

siRNAs act in defending against viral infections and transposon mobility (Piatek and Werner 2014;Zhu and Palli 2020).The study found the numbers of siRNA pathway genes were variable.More expansions ofAgo-2have been recorded in five orders,Dipteran,Hemipteran,Orthopteran,Coleopteran and Hymenopteran (Guoet al.2015;Dowlinget al.2016;Mongelli and Saleh 2016;Tianet al.2019).R2D2is present only in winged insects and is viewed as a derived feature of Pterygota (Dowlinget al.2016).ThreeAgo-2and twoR2D2genes inP.puparumwere identified in this study.TheAgo-2genes with high identities were tandemly arranged in theP.puparumandN.vitripennisgenomes.Here,it was conjectured that the expansions originated from tandem duplication which happened in the ancestors of both wasp species.Conversely,theAgo-2genes were separate in the caste-forming polyembryonic wasp,C.floridanumand the coleopteran,T.castaneum.It was concluded that there were different gene duplication mechanisms.Two R2D2 proteins had low identity to each other,though both possessed two dsRBDs at the N-terminal,as seen also in the two R2D2 proteins ofT.castaneum(Tomoyasuet al.2008).InDrosophila,R2D2was regarded as one of the most rapidly evolving genes (Obbardet al.2006).The low identities among proteins could be the result of rapid evolution,noting that theR2D2duplications occurred in several other hymenopteran species in this study.

The miRNA and siRNA pathways are related because both ncRNAs are processed by Dicer/dsRBP complex for gene silencingviainteractions with Argonaute proteins(Carthew and Sontheimer 2009).The two pathways are distinguished because they recruit different genes.Dicer-1,LoquaciousandAgo-1operate in the miRNA pathway andDicer-2,R2D2andAgo-2in the siRNA pathway.Previous studies revealed substantially accelerated protein evolution rates ofDicer-2andAgo-2inDrosophilaspecies and inAedesaegypti,compared to their paralogs in the miRNA pathway (Obbardet al.2006;Bernhardtet al.2012).This might be the outcome of co-evolution between RNA viruses and the host siRNA pathway.The HymenopteranDicer-2andAgo-2genes have significantly higher ω values compared to their paralogs,and the ω values of the threeP.puparumAgo-2genes are approximately 12 times greater compared toAgo-1.The ω value for theP.puparumDicer-2is lower thanDicer-1,most likely due to its higher synonymous substitutions rate.

Generally,higher miRNA expression was recorded in adults than in juveniles,which suggested higher miRNA abundances in the adults.Since miRNAs act in many areas of insect biology according to earlier studies (Iovinoet al.2009;Lucas and Raikhel 2013;Zhanget al.2016;Linget al.2017),it is speculated that miRNAs play specific roles,such as influencing ovarian development.Similarly,a peak ofAgo-1expression appeared in late larvae,which may act in metamorphosis.Expression patterns differed among siRNA pathway genes,even among the duplicate genes.For example,R2D2-1was expressed more highly in early larvae whileR2D2-2expression was higher in adults compared to other stages.Likewise,Ago-2aexpression was higher in pupae and adults andAgo-2cexpression was higher in larvae.Alternatively,Ago-2genes may have different functions inP.puparum,since theArgonautegenes have expanding roles in RNA alternative splicing (Maet al.2018).The roles ofAgo-2genes inP.puparumneed further investigation.

5.Conclusion

This study identified and recorded expression patterns of the genes associated with ncRNAs inP.puparum.The major significance of this work is proposing a series of hypotheses that will guide future work designed to understand the actions of individual genes in the miRNA and siRNA pathways.

Acknowledgements

The study was funded by the Key Program of National Natural Science Foundation of China (31830074),the Program for Chinese Outstanding Talents in Agricultural Scientific Research of the Ministry of Agriculture and Rural Affairs of China,the Program for Chinese Innovation Team in Key Areas of Science and Technology (2016RA4008).

Declaration of competing interest

The authors declare that they have no conflict of interest.

Appendicesassociated with this paper are available on http://www.ChinaAgriSci.com/V2/En/appendix.htm