Introduction
Insects are the only invertebrates with the ability to fly. Their wings are not only the core structure of flight but also play an important role in orientation, protection, communication, and courtship (Shyy et al., Reference Shyy, Kang, Chirarattananon, Ravi and Liu2016). For a long time, the molecular mechanism of insect wing development has been an area of particular interest. Holometabolous insects need to undergo four developmental stages to become adults (Rolff et al., Reference Rolff, Johnston and Reynolds2019). Previous studies have shown that some genes are involved in the early development of Drosophila wings. These genes include transcription factors, such as engrailed (en), serum response factor (SRF) (Simmonds et al., Reference Simmonds, Liu, Soanes, Krause, Irvine and Bell1998), apterous (ap) (Bourgouin et al., Reference Bourgouin, Lundgren and Thomas1992), scalloped (sd) (Celniker and Rubin, Reference Celniker and Rubin2003), and vestigial (vg) (Williams et al., Reference Williams, Paddock, Vorwerk and Carroll1994). Studies have shown that BGIBMGA000319 gene encodes Cys2-His2 (C2H2) zinc finger protein, which is homologous to B lymphocyte maturation induction protein-1 (BLIMP-1) gene (Wu et al., Reference Wu, Tong, Li, Ding, Zhang, Fang, Tan, Hu, Liu and Dai2019). BMBLIMP-1 of silkworm has been specifically up-regulated during metamorphosis and has been highly expressed in the wings of the pupa on the third day of the pupal stage, and RNA interference (RNAi) has led to wing abnormalities (Wu et al., Reference Wu, Tong, Li, Ding, Zhang, Fang, Tan, Hu, Liu and Dai2019). Because BMBLIMP-1 is a C2H2 zinc finger gene, it may be related to the development of silkworm pupa wings.
Zinc finger protein is a type of transcription factors with a finger domain, playing an important role in gene expression regulation, cell differentiation, embryonic development, apoptosis, and other important life processes. In 1985, the zinc finger structure was first discovered in the transcription factor TFIIIA of Xenopus laevis (Miller et al., Reference Miller, McLachlan and Klug1985). The structure of zinc finger protein is defined as the existence of several Cys residues in an amino acid sequence of the regulatory protein, which is a zinc-dependent DNA-binding domain forming a short and stable finger-like structure by binding Zn2+ for self-folding (Frankel and Pabo, Reference Frankel and Pabo1988). Since then, researchers have found a large number of proteins containing zinc finger structure, collectively known as zinc finger proteins. According to the conserved motif characteristics of the zinc finger, Berg and Shi (Reference Berg and Shi1996) divided zinc finger proteins into nine categories: C2H2 (Cys2His2), C8 (Cys8), C6 (Cys6), ring type (Cys3HisCys4), Cys2HisCys, Lim type (Cys2HisCys5), C4 (Cys4), and Cys3His. Multiple types of zinc finger proteins may also be expressed in the same tissue (Miller et al., Reference Miller, McLachlan and Klug1985). Zinc fingers protein can bind not only DNA and RNA, but also DNA–RNA hybrid double-strand molecules and other zinc finger proteins or themselves to regulate gene expression at the transcriptional and translational levels (Zhao and Zhao, Reference Zhao and Zhao2009).
The different expression patterns of BgZelda in Blattella germanica seem to be related to metamorphosis (Ventos-Alfonso et al., Reference Ventos-Alfonso, Ylla and Belles2019). Holometabolous species, such as Drosophila melanogaster, produce larvae different from those of adults during embryogenesis. Therefore, in addition to MZT, the expression of Zelda in D. melanogaster may activate the continuous gene set needed to establish the morphology of vermicular larvae during embryogenesis (Ventos-Alfonso et al., Reference Ventos-Alfonso, Ylla and Belles2019). The Zelda gene is also expressed in later stage embryonic development of Apis mellifera, which is related to the development of precursor cells of the central nervous system (Pires et al., Reference Pires, Freitas, Cristino, Dearden and Simões2016). Other studies have found that Zelda is continuously expressed during embryonic development, including the developing central nervous system and trachea (Pearson et al., Reference Pearson, Watson and Crews2012). Additionally, it was reported that 55 and 58 CCCH zinc finger protein genes were identified in humans and mice, respectively. CCCH genes are expressed at different levels in all tissues. They are highly expressed in tissues with high energy consumption, such as the heart, muscle, blood vessels, and brown adipose tissue (Liang et al., Reference Liang, Song, Tromp, Kolattukudy and Fu2008b).
Zinc finger transcription factor Zelda (also known as vielfaltig) plays an important role in the transition from the mother cell to the zygotic cell (MZT) in D. melanogaster. In D. melanogaster, mutant embryos lacking Zelda have defects in the cell germ formation and cannot activate genes necessary for cell formation, sex determination, and pattern formation (Liang et al., Reference Liang, Nien, Liu, Metzstein, Kirov and Rushlow2008a; Schulz et al., Reference Schulz, Bondra, Moshe, Villalta, Lieb, Kaplan, McKay and Harrison2015). Moreover, Zelda studied in this paper is an important factor regulating early gene expression and plays a key role in the activation of the early zygotic genome of D. melanogaster (Liang et al., Reference Liang, Nien, Liu, Metzstein, Kirov and Rushlow2008a). T. castaneum is a holometabolous insect with short embryonic development, and Zelda also plays a key role in the posterior division and configuration of T. castaneum adult disk structure (Dönitz et al., Reference Dönitz, Schmitt-Engel, Grossmann, Gerischer, Tech, Schoppmeier, Klingler and Bucher2015; Ribeiro et al., Reference Ribeiro, Tobias-Santos, Santos, Antunes, Feltran, de Souza Menezes, Aravind, Venancio and Nunes da Fonseca2017; Ventos-Alfonso et al., Reference Ventos-Alfonso, Ylla and Belles2019). However, regulatory network about how TcZelda performs these functions in insects is still unclear. To further explore the signaling systems of TcZelda involved in insect physiology, RNA sequencing and RNA interference were combined to study the potential function of TcZelda in T. castaneum. The results will shed a new light into molecular regulatory mechanisms of other relative genes in T. castaneum.
Material and methods
Experimental samples
T. castaneum Georgia-1 (GA-1) strain was used for all experiments. It was placed in a flask with a mass ratio of flour to yeast powder at 19:1 and cultured in an artificial climate box with a set temperature of 30°C and a humidity of 40% (Begum et al., Reference Begum, Li, Beeman and Park2009).
Expression profiling of TcZelda gene
The tissue (head, antennae, ovary, testis, epidermis, midgut, and fat body) were dissected from 100 10-day-old adults and eight developmental stages were extracted by Trizol method. The eight stages of development included early eggs (1 day old, ~50 mg), late eggs (3 days old, ~50 mg), early larvae (1 day old, 3 individuals), late larvae (15 days old, 3 individuals), early pupae (1 day old, 3 individuals), late pupae (5 days old, 3 individuals), early adults (1 day old, 3 individuals), and late adults (10 days old, 3 individuals). Finally, the extracted RNA was reversed into circular DNA (cDNA) by reverse transcriptional reagents 4× gDNA wiper Mix and 5× HiScript Q RT Supermix II. Quantitative reverse transcription polymerase chain reaction (qRT-PCR) was used to analyze cDNA in different developmental stages and tissues. The reaction system consisted of 0.25 μl forward and backward primers (10 μM), 2 × 5 μl AceQ Universal SYBR qPCR Master Mix, 3.5 μl ddH2O, and 1 μl cDNA. Three technical and biological repeats were performed for each group of samples, and the melting curve was analyzed to ensure that only one PCR product was amplified.
Synthesis and injection of dsRNA
dsRNA of TcZelda and VER (T. castaneum vermillion, GenBank AY052390) was synthesized with primers designed by Primer Premier 5.0. The gene primers are shown in table S1. The promoter of T7 polymerase was added to the 5’ end of two forward primers, and dsRNA was then synthesized by TranscriptAi T7 high yield transcription kit. Using Nanoject III instruments (Drummond Scientific, Philadelphia, USA), the obtained dsRNA (200 ng/150 nl) was injected into the body cavity of each larva (15 days old), and the same volume of dsVER or injection buffer (IB) was microinjected as negative control and blank control, respectively. In each experiment, 30–40 larvae were injected, and the biological repeats were performed three times. The larvae injected with dsRNA were cultured normally, and the survival rate was recorded on the fourth day. Three individuals of each group were collected on day 4 after dsRNA injection to examine the knockdown level of the targeted gene by qRT-PCR.
Images of adults were taken using a Zeiss SteREO Discovery V8 (Zeiss, Oberkochen, Germany) stereomicroscope. ZEN2.3 microscope software (Blue edition, Zeiss) was used to control the microscope, image acquisition, and exportation of TIFF files. The images were processed using Photoshop CS6 software (Adobe Systems Inc., San Jose, CA, USA) for converting RAW format to EPS format, which then were compatible with the Adobe Illustrator version 6.0 software (Adobe Systems Inc.) vector graphics image editor program.
Extraction and sequencing of total RNA
Six samples, including two control and two dsZelda-treated samples harvested on the fourth day after dsRNA injection (each sample contained six larvae), were selected for RNA-Seq analysis. The concentration and purity of RNA were detected by Nanodrop2000/2000c. Only RNA with A260/A280 > 1.80 was used for the follow-up study, and the integrity of RNA was detected by Agilent 2100 Bioanalyzer. The reverse transcription was then performed with 1000 ng of total RNA (A260/A280 ratio >1.8 and A260/A230 ratio >2.0) using HiScript II reverse transcriptase (Vazyme, Nanjing, China), following the manufacturer's manual. Samples of RNA with established purity and integrity were sequenced at Beijing Genomics Institute (Shenzhen, China) on the BGISEQ-500 platform. The sequences were deposited at the National Center for Biotechnology Information (NCBI) database under accession numbers SAMN28580502 and SAMN28580503 for IB samples and SAMN2857338 and SAMN2857338 for dsZelda samples.
Filtering and mapping of RNA-Seq data
The original sequencing data contained reads with low quality, joint contamination, and high unknown base N content. These reads needed to be removed before data analysis to ensure the reliability of the results. SOAPnuke, a BGI-developed filtering software, was used for statistical analysis, while Trimmomatic was used for filtering. Specific steps included the removal of the reads containing the joint (joint contamination), the removal of the reads with an unknown base N content greater than 5%, and the removal of low-quality reads (defined as the reads with a mass value <15 that represented >20% of the total base number of the reads). The resulting clean reads were aligned with the reference genome sequence using HISAT (Kim et al., Reference Kim, Langmead and Salzberg2015) and with the reference gene sequences using Bowtie2 (Li and Dewey, Reference Li and Dewey2011; Langmead and Salzberg, Reference Langmead and Salzberg2012).
Differential gene expression analysis
To authenticate differentially expressed genes (DEGs) in the treated and control samples, RPKM values (reads per kilobase per million mapped reads) were used to calculate gene expression levels, and DEGs were defined as genes with |log2 − (T/C)| ≥ 0.5 and adjusted P-values of <0.05 based on DEseq2. DEGs were analyzed by hierarchical clustering using the heatmap function in R software. Furthermore, functional analysis and annotation of DEGs in T. castaneum were performed using the Gene Ontology (GO) database, which included molecular functions, cellular components, and biological processes, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.
Validating RNA-sequencing data by qRT-PCR
On the fourth day after dsRNA injection, the total RNA was extracted from six larvae of dsZelda or control groups, and two biological replicates were carried out for each sample. Total RNA of larval samples in the control and dsZelda groups was extracted by RNAisoTMPlus (TaKaRa), and then 1 μg of total RNA was reverse transcribed to synthesize the first strand of cDNA. RNA-sequencing data were verified by qRT-PCR according to the instructions of FastStart SYBR Green Master (Roche). Tcrps3 (GenBank accession number was XM_965494) was used as the internal reference gene (Horn and Panfilio, Reference Horn and Panfilio2016), and the gene expression was analyzed by the 2−ΔΔCT method (Livak and Schmittgen, Reference Livak and Schmittgen2001). The experiment was biologically repeated at least three times. The gene primers are shown in table S1.
Statistical analysis
Gene expression was calculated by the 2−ΔΔCT method (Livak and Schmittgen, Reference Livak and Schmittgen2001), and the gene expression data, the mean values obtained for the RNAi-treated larvae, and the mean values obtained for the control larvae were compared by one-way analysis of variance in combination with Fisher's least significant difference multiple comparison test using SPSS version 13.0 (SPSS Inc., Chicago, IL, USA). All data were presented as the mean ± standard error (SE). The error bars represented SEs among three biological replications. P-values of <0.05 were regarded as statistically significant, and P-values of <0.001 were considered extremely significant.
Results
Developmental stage and tissue-specific expression profile of TcZelda
The qRT-PCR results revealed that transcripts of TcZelda were detected at all stages from embryo to adult and strongly expressed at the early egg stage (fig. 1a). The expression of TcZelda in different tissues of larvae and adults was further studied. In the larval tissue, compared with the whole larva, the expression of TcZelda was the most significant in the head, followed by the epidermis, fat body, hemolymph, and the intestine, where it was the lowest (fig. 1b). In adult tissues, compared with whole adult beetles, the expression of TcZelda was the most significant in the ovary, followed by fat body, head and epidermis, and very low in the gut, antennae, testis, and accessory glands (fig. 1c).
Figure 1. Expression of TcZelda gene in the developmental stage (A), different tissues of larvae (B) and adults (C) of T. castaneum. Early egg (1 day old), late egg (3 days old), early larva (1 day old), late larva (20 days old), early pupa (1 day old), late pupa (5 days old), early adult (1 day old), and late adult (10 days old). Different letters on the column chart indicate that there are significant differences in different developmental stages. Late larvae (15 days old) and adults were dissected under the somatotype microscope, including the whole larva, head integument, fat body, gut, and hemolymph. Additionally, the whole adult, head, integument, fat body, gut, ovary, antennae, testis, and accessory gland were dissected. The transcriptional level of Tczelda in the whole larva and adult was used to calibrate the tissue-specific expression profile. The vertical line is the standard error of the mean (n = 3), and the vertical line letters indicate that the mean expression values at different developmental stages are significantly different at P < 0.05.
Observation of biological characterization after RNAi
TcZelda gene was silenced by RNAi to further determine the function of TcZelda in insects. After injection on the third day, the silence efficiency of TcZelda was measured by qRT-PCR, and the transcript levels of TcZelda were reduced by about 70% as compared to control (fig. 2a). We counted the mortality rate after RNAi, of which the mortality rate of insects injected with IB was 9%, that of insects injected with dsVER was 9.5%, and that of insects injected with dsZelda was 95% (fig. 2b). In addition, the death state of adults included the inability of wings to close, shortening of the antennae, and possible fusion of the ankles of the legs (fig. 2c).
Figure 2. The dsZelda-mediated suppression of TcZelda transcripts as determined by qRT-PCR (a). Mortality rate of the beetles after larval RNAi (b) phenotype (c). The insects in the control group were injected with the same amount of IB or dsVER. The error bar represents the standard error of three biological repeats. The letters at the top of the bar chart indicate that there are significant differences at a level of P of <0.05 between different treatment groups. Mortality statistics of 4 days after injection of IB, dsVER, or dsZelda (b). The data in the figure are the mean ± standard error (nasty 3). Different letters on the column chart indicate that there are significant differences among different treatments at P < 0.05 level (scale bars: 1000 μm for panel C).
Transcriptome profiling by RNA sequencing
The transcriptional group was sequenced, and 84.20 million clean reads were obtained from the control group and dsZelda-treated samples. Each sample has a clean reading of about 40 million, which is enough for quantitative analysis of gene expression. In the sequencing analysis, the control group produced 42.64 million raw reads, while the treatment group produced 43.82 million raw reads, of which about 96% passed the CleanReadsQ20 (table 1). The ratio of clean reads to raw reads of the four gene libraries was between 97.24 and 97.71%. The SOAPaligner/SOAP2 software was used to compare the sequenced readings from these sequences with the reference genome database of T. castaneum. In total, 85.24% of the control reads could match the genome, of which 47.37% had only one match with the genome, while 86.29% of the dsZelda groups readings could match the genome, and the only matching proportion in the genome was 48.39% (table 1). About 15% of the reads in the two gene libraries did not match the reference genome. The Pearson correlation coefficients of all gene expressions between every two samples were calculated and reflected in the form of a heat map to reflect the correlation of gene expression between samples. The correlation coefficient can reflect the similarity of overall gene expression among samples. The higher the correlation coefficient is, the more similar is the gene expression level (fig. S1). In sum, these data suggest that the two biological replicates were highly correlated for both the control and treated samples.
Table 1. Summary statistics of read numbers based on the T. castaneum transcriptome
Functional analysis of DEGs based on RNA-seq data
Fragments per kilobase million (FPKM) values of unigenes from the control and treatment groups with knockdown target gene were compared to further explore the function of TcZelda in T. castaneum. A total of 592 DEGs were identified, including 277 up-regulated genes and 315 down-regulated genes (fig. 3a, table S2). The cluster heat map of differential gene expression showed the similarity and difference of DEGs expression levels of two biological repeats from the control and dsZelda groups (fig. 3b). Ten genes (5 up-regulated genes and 5 down-regulated genes) were randomly selected and evaluated by qRT-PCR to further verify the accuracy of transcriptome data and DEGs results (table S1). The randomly selected DEGs in qRT-PCR and transcriptome data (fig. 3c) had a similar adjustment trend, indicating that the data of RNA-seq was reliable.
Figure 3. Analysis of differential gene expression after TcZelda gene knockdown. The number of differential genes (a) and differential gene expression clustering heat map (b) between the dsZelda and control groups. Comparison of gene expression profiles between RNA sequencing and qRT-PCR detection (c). The row represents a single gene, and the column represents the comparison of genes between the dsZelda treatment and control groups. The FPKM of the sample obeys standard normal distribution. The column on the right clusters the samples (N = 6) based on the similarity of Log2 (FPKM + 1). Group C was the control group, and group T was the dsZelda treatment group.
To annotate these DEGs, we performed a GO functional analysis (fig. 4, table S3). Gene Ontology can be divided into three functional categories: molecular function, cellular component, and biological process. According to the functional classification of differential genes, 592 differential genes (DEGs) were divided into 9, 3, and 17 main functional groups, respectively. Cellular anatomical entity (GO:0110165) with 322 genes (54.39%) and binding (GO:0005488) with 235 genes (39.70%) are the major biological processes. Cellular (GO:0009987) and metabolic processes (GO:0008152), containing 201 genes (33.95%) and 130 genes (21.96%), were dominant in the main categories of biological process, respectively. Additionally, GO term enrichment analysis showed that DEGs were associated with 30 GO terms (table S3), and the most enriched GO terms were oxidoreductase activity (GO:0016491), developmental process (GO:0032502), anatomical structure development (GO:0048856), and multicellular organism development (GO:0007275) (fig. S2). Further, KEGG enrichment analysis showed that there were 20 pathways significantly enriched in DEGs compared with the whole transcriptome background with a Q value of ≤0.05, and the most enriched pathways were lysosome (ko04142), protein digestion and absorption (ko04974), and focal adhesion (ko04510) (fig. 5, tables S4, S5)
Figure 4. GO analyses of DEGs in dsZelda and control samples. Different colors represent GO categories. The x-axis shows the number of DEGs. The y-axis shows the GO terms.
Figure 5. KEGG pathway analyses of DEGs in dsZelda and control samples. Pathway significance is shown together with the Q-value (color), rich factor (vertical ordinate), and the number of genes (size of circles).
Genes related to embryonic development, wing development, and fatty acid metabolism
In these DEGs databases, there were six genes related to embryonic development, including altered disjunction (ALD), TGF-beta, propeptide (TGF), AGAP005368-PA (GAP), and genes encoding cilia- and flagella-associated protein 44, advillin, and condensing complex subunit 1 (table 2). Additionally, nine DEGs were involved in the development of insect wings (table 3), and 11 DEGs were involved in the fatty acid metabolism of insects (table 4). Ubx ultrabithorax (UBX), Cx cephalothorax (CX), optomotor-blind-like (OMB), and five Endocuticle structural glycoprotein sgabd (ABD) genes were involved in the development of insect wings (table 3). The proteins involved in fatty acid metabolism included delta-12 desaturase, elongation of very long-chain fatty acids protein AAEL008004, fatty acid synthase, and palmitoyl-protein thioesterase 1 (table 4).
Table 2. Genes involved in reproductive process
C, control; T, dsZelda.
a Metamorphosis and survival after RNAi were analyzed (https://ibeetle-base.uni-goettingen.de/).
b Lethalities 11 days after larval injection (includes death as larva, prepupa, pupa).
Table 3. Genes involved in the development of insect wings
C, control; T, dsZelda.
a Metamorphosis and survival after RNAi were analyzed (https://ibeetle-base.uni-goettingen.de/).
b Lethalities 11 days after larval injection (includes death as larva, prepupa, pupa).
c 13 days after female pupal injection.
d 16 days after larval injection.
Table 4. Genes involved in fatty acid metabolism
C, control; T, dsZelda.
Discussion
Previous studies have shown that Zelda is expressed during MZT in D. melanogaster and T. castaneum. In these species, Zelda even plays a role in the post-embryonic development. This Zelda function beyond the range of MZT might be related to the evolution and innovation of holometabolous insects. The expression of Zelda related to blastoderm and gastrula formation has been observed in the early embryonic development of western honeybees (Pires et al., Reference Pires, Freitas, Cristino, Dearden and Simões2016), indicating that Zelda plays the same role as the activator of the early zygotic genome described in D. melanogaster (Liang et al., Reference Liang, Nien, Liu, Metzstein, Kirov and Rushlow2008a; Ventos-Alfonso et al., Reference Ventos-Alfonso, Ylla and Belles2019). Our results showed that TcZelda was expressed in eight key developmental stages (fig. 1a), indicating that TcZelda may be involved in various physiological processes. TcZelda highly expressed in embryos and larvae, then decreased expression in pupae, and slightly increased expression in late adults. It is speculated that TcZelda gene plays a key role in the activation of the early zygotic genome and late embryonic development. Through the analysis of our experimental results, TcZelda was mainly highly expressed in the head of larvae (fig. 1b), and it is speculated that TcZelda may be involved in the development of precursor cells in the central nervous system. It was found that TcZelda was also expressed in the epidermis, fat body, and hemolymph (fig. 1c); thus, we speculated that TcZelda gene might be involved in the regulation of fatty acid metabolism. The results showed that the expression of TcZelda in the ovary was the highest in adult tissues (fig. 1c), and the gene was probably related to the activation of the zygotic genome during the transition from mother cell to zygotic cell in early embryonic development.
Keratin and chitin are the main components of the insect exoskeleton and play an important role in wing development (Iconomidou et al., Reference Iconomidou, Willis and Hamodrakas2005). By analyzing the experimental results, it was found that TcZelda was expressed in the exoskeleton of adults (fig. 1c). Thus, TcZelda may play an important role in providing the stability of the stratum corneum during wing development. Moreover, TcZelda was expressed in the head and fat body (fig. 1B). It is speculated that TcZelda may also be involved in the metabolism of fatty acids in larval tissue. Additionally, the analysis of the experimental results showed that the expression and function of TcZelda significantly exceeded that of embryos in all metamorphosis species (fig. 1A), which was proved by its role in the disk structure formation of adults (Ribeiro et al., Reference Ribeiro, Tobias-Santos, Santos, Antunes, Feltran, de Souza Menezes, Aravind, Venancio and Nunes da Fonseca2017) and neuroblasts in D. melanogaster larvae (Reichardt et al., Reference Reichardt, Bonnay, Steinmann, Loedige, Burkard, Meister and Knoblich2018). Studies have shown that Zelda gene expression in D. melanogaster and T. castaneum extends to the middle and late stages of embryonic development, which may be a factor leading to the total metamorphosis innovation of the semi-abnormal ancestors of endoptera insects (Ventos-Alfonso et al., Reference Ventos-Alfonso, Ylla and Belles2019).
This study suppressed the expression of TcZelda gene by using RNAi. Four days after the injection of dsZelda, qualitative analysis was performed with RT-PCR, and then interference efficiency was quantitatively analyzed with qPCR. Compared with the control group, the expression of TcZelda was significantly down-regulated, and the interference efficiency was about 70% (fig. 2A). After the expression of TcZelda gene was suppressed, the mortality rates at different developmental stages were statistically analyzed. The cumulative mortality rate of the treatment group injected with dsZelda was 97%, including 72% in the pupa stage and 23% in the adult stage (fig. 2B). Therefore, we speculate that TcZelda might play an important role in the physiological and biochemical response of T. castaneum pupal stage. After the expression of the gene was suppressed, through the observation and analysis of adult biological characteristics, it was found that the wings of the adults were not closed, the antennae were shortened, the reproductive leaves of females were reduced, and the tarsal joints of the two legs were fused (fig. 2c), indicating that TcZelda also plays a key role in the posterior segmentation and configuration of the adult disk structure of the whole metamorphosis species with short embryonic development (Ventos-Alfonso et al., Reference Ventos-Alfonso, Ylla and Belles2019). However, in the course of the experiment, due to the death of a large number of pupae after RNAi, a small number of adults with successful eclosion had a wide gap between male and female insects; hence, it was difficult to perform phenotypic and reproductive statistics of the next generation. Therefore, the experiment needs further exploration.
By inhibiting the expression of vielfaltig gene (also known as Zelda) by RNAi technology, it can inhibit the growth of the stratum corneum in the egg stage, decrease the fecundity of larvae after injection, and die in the pupal stage, which is consistent with the results of iBeetle-Base (Dönitz et al., Reference Dönitz, Schmitt-Engel, Grossmann, Gerischer, Tech, Schoppmeier, Klingler and Bucher2015). Thus, Zelda plays a key role in the embryonic development, wing development, and fatty acid metabolism of insects. We used an IlluminaHisSeq2000 for RNA sequencing to fully understand the role of Zelda gene in reproductive development and pupa lethality in T. castaneum. After TcZelda gene knockdown, according to the functional annotation of the gene, 592 DEGs were assigned to 29 major functional groups (fig. 4, table S3), and 20 pathways were significantly rich in DEGs (fig. 5, table S4), of which 277 and 355 DEGs were up-regulated and down-regulated (fig. 3c, table S2), respectively. These results provide a basic understanding for the study of the role of Zelda gene in the functional regulation system of T. castaneum.
Zinc finger protein Zelda plays a key role in the activation of the insect's early zygotic genome. Six genes, including GAP, TGF, and ALD, were found by comparing the control group with the dsZelda group (table 2). Altered disjunction (ALD) is a magnetic-specific meiotic mutation-related protein. Studies have shown that the gene product of ALD plays a role in preventing the exchange of X chromosome and fourth chromosome in non-homologous segregation events (O'Tousa, Reference O'Tousa1982). It can be concluded that ALD is a key gene in the insect reproductive system. GAP genes can limit the region of the somatic segment of the embryo, and the mutation of these genes can change the number or polarity of the somatic structure (Boos et al., Reference Boos, Distler, Rudolf, Klingler and El-Sherif2018). Gene products of GAP divide embryos into regions equivalent to three body segments. Additionally, products of spacer genes form concentration gradients in their respective expression regions, and the location information provided by these gradients is transmitted to downstream genes, which orderly regulate the formation of embryonic development patterns and finally determine the boundaries between body segments and the fate of cells in different parts of the embryo (Sánchez and Thieffry, Reference Sánchez and Thieffry2001). After TcZelda was knocked down, the expression of ALD and GAP was down-regulated, which indicating that Zelda and ALD and GAP might be upstream and downstream genes in the reproductive regulatory system, respectively.
The TGFβ superfamily is composed of a class of polypeptide growth factors related to the structure and function of different tissues, including TGFβ, activin, inhibin, and BMP (Li et al., Reference Li, Lin, Geng and Xu2018). Members of the TGFβ family exist widely in various tissues of many organisms and play an important role in the development of organisms and organs. Reproductive stem cells in Drosophila ovary are known to be regulated by the Dpp signaling. The gene related to the survival of early male germ cells is DPP, and the gene related to the proliferation and differentiation of early male germ cells is BAM. Studies have shown that TGFβ signal transduction may affect the expression of these genes, regulating the proliferation and differentiation of germ cells (Chen and McKearin, Reference Chen and McKearin2003; Song et al., Reference Song, Wong, Kawase, Xi, Ding, McCarthy and Xie2004). After TcZelda knockdown, the expression of TcZeldaβ was up-regulated, indicating that the gene may complement the function of TcZelda. The above results confirmed that TGFβ was involved in insect embryonic development, which was consistent with the results of the previous tissue expression profile.
There were some genes related to the development of insect wings. As the only and earliest winged group of invertebrates, the flying ability of insects allows them to have great advantages in migration, foraging, and courtship (Guo et al., Reference Guo, Dan, Li, Jie and Entomology2017). Some studies have shown that the expression of UBX and ABD-A is suppressed by RNAi. There are 20 sheath-wing discs in ten body segments from the middle chest to the 8th ventral segment (T2–T8), which is the same as that of Drosophila. UBX is used to distinguish the front and rear wings, except that in the hind wings, UBX forms membrane wings by inhibiting the specific genes of the sheath wing, such as stratum corneum thickening genes, and the hind wings are specialized into sheath wings after being knocked out (Tomoyasu et al., Reference Tomoyasu, Wheeler and Denell2005). As a copy of SCR, the loss of function of Cephalo-thorax (CX) can lead to the formation of additional tissue similar to the sheath wing structure in the forechest (T1) (Guo et al., Reference Guo, Dan, Li, Jie and Entomology2017). After the injection of dsRNA of TcZelda, the expression of UBX and ABD genes decreased, while among the five ABD genes, the expression of ABD-8, ABD-5, and ABD-4 was significantly up-regulated, and that of ABD-1 and ABD-2 was significantly down-regulated (table 3). Therefore, we speculate that the function of TcZelda gene is similar to that of these genes, this is supported by the fact that it is also similar to the results of the early development stage and tissue expression profile. TcZelda plays an important role in the development of insect wings.
As mentioned above, we found that TcZelda is also expressed in the epidermis, fat body, and hemolymph; hence, we speculate that TcZelda gene may be involved in the regulation of fatty acid metabolism. In the KEGG enrichment pathway, there were also 11 DEGs related to insect fatty acid metabolism, including ELO, FAD, and FAS (table 4). FAS was located in the cytoplasm and can catalyze the reaction of acetyl-CoA with malonic acyl-CoA, while palmitic acid can be obtained by the reddest reaction (Jakobsson et al., Reference Jakobsson, Westerberg and Jacobsson2006). The study of FAS in insects focused on spatiotemporal expression and function and its effect on insect diapause. Some studies have shown that RNAi with FAS can reduce the content of FAS and fecundity in the ovary of brown planthopper (Lei et al., Reference Lei, Jiang, Liu, You and Wu2016). After RNA interferes with FAS-1 and FAS-3 of female mosquitoes, female mosquitoes cannot store lipids needed for winter; therefore, these two FAS play an important role in fat storage in the early diapause stage (Sim and Denlinger, Reference Sim and Denlinger2009). ELO is the rate-limiting condensase obtained in the first step of the fatty acid extension reaction (Leonard et al., Reference Leonard, Pereira, Sprecher and Huang2004). The study of ELO in insects has focused on Drosophila melanogaster, including the effect of ELO on reproductive function and the role of ELO in pheromone synthesis and the function of the epidermis. Some studies have shown that directional RNA interferes with NOA (ELOCG3971) in the late developmental stage of Drosophila testis germ cells, which can lead to male Drosophila infertility. Moreover, NOA damages the motor function of D. melanogaster and significantly reduces its survivability (Jung et al., Reference Jung, Hollmann and Schäfer2007).
FAD is an enzyme that can catalyze the formation of double bonds on specific sites of fatty acyl-CoA. Acyl-CoA Delta-11 desaturase belong to intact membrane proteins and is a type of fatty acyl-CoA desaturase (FAD) (Haritos et al., Reference Haritos, Horne, Damcevski, Glover and Gibb2014). FAD of insect has been widely studied, mainly focusing on the role of FAD in pheromone synthesis. By feeding FAD inhibitor CAY10566 to Drosophila, Wang et al. reduced the unsaturated fatty acids and slowed down the growth of Drosophila. After the first molting, the expression of ecdysone regulatory protein was no longer induced, normal molting could not take place, and the insect died a few days later (Wang et al., Reference Wang, da Cruz, Pulfemuller, Grégoire, Ferveur and Moussian2016). Additionally, FAD at different sites can also affect the diversity of insect pheromones (Albre et al., Reference Albre, Liénard, Sirey, Schmidt, Tooman, Carraher, Greenwood, Löfstedt and Newcomb2012). At the same time, 56 genes were down-regulated in TcZedla knockdown beetles, and these genes could lead to impaired embryonic development of insects, suggesting that TcZedla gene may regulate embryonic development of TcZedla through these genes in T. castaneum (table S5). The up-regulation or down-regulation of DEGs expression also indicates that TcZelda may have similar functions, which further confirms the important role of TcZelda in embryonic development and reproduction and provides clues for the high expression of TcZelda in the epidermis, fat body, and hemolymph. Still, the specific mechanism needs to be further studied.
Conclusion
There were 592 DEGs between ds-TcZelda and control insects. These results suggest that TcZelda may not only regulate insect embryonic development through ALD, GAP, and TGF genes but also participate in insect wing development through UBX, CX, and ABD. Additionally, the knockdown of TcZelda gene led to significantly different expression of genes, such as FAS, ELO, and FAD, which are involved in fatty acid metabolism, might be related to the high expression of TcZelda in insect head, epidermis, fat body, and hemolymph, and might be involved in insect reproduction and pheromone synthesis. These results may help to find new pest target genes and more suitable methods for pest control.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485323000263.
Financial support
This study was sponsored by the Staring Foundation for the Doctor, Anyang Institute of Technology (grant number BSJ2019009, BSJ2019022 and BSJ2021040); the Staring Foundation of Innovation and Practice Base for Postdoctors, Anyang Institute of Technology (grant numbers BHJ2020008); the Scientific and Technological Project of Henan Province (grant number 212102110444).
Competing interest
None.
