Evolution and multiple roles of the Pancrustacea specific transcription factor zelda in insects

Gene regulatory networks (GRNs) evolve as a result of the coevolutionary processes acting on transcription factors (TFs) and the cis-regulatory modules they bind. The zinc-finger TF zelda (zld) is essential for the maternal-to-zygotic transition (MZT) in Drosophila melanogaster, where it directly binds over thousand cis-regulatory modules to regulate chromatin accessibility. D. melanogaster displays a long germ type of embryonic development, where all segments are simultaneously generated along the whole egg. However, it remains unclear if zld is also involved in the MZT of short-germ insects (including those from basal lineages) or in other biological processes. Here we show that zld is an innovation of the Pancrustacea lineage, being absent in more distant arthropods (e.g. chelicerates) and other organisms. To better understand zld´s ancestral function, we thoroughly investigated its roles in a short-germ beetle, Tribolium castaneum, using molecular biology and computational approaches. Our results demonstrate roles for zld not only during the MZT, but also in posterior segmentation and patterning of imaginal disc derived structures. Further, we also demonstrate that zld is critical for posterior segmentation in the hemipteran Rhodnius prolixus, indicating this function predates the origin of holometabolous insects and was subsequently lost in long-germ insects. Our results unveil new roles of zld in different biological contexts and suggest that changes in expression of zld (and probably other major TFs) are critical in the evolution of insect GRNs.


Introduction
Gene regulatory networks (GRNs) depend on the coevolution of transcription factors (TFs) and their relationship with the cis-regulatory modules (CRMs) they bind [1,2]. In insects, the detailed role of a number of TFs and CRMs have been well-described, particularly during the embryogenesis of the fruit fly Drosophila melanogaster [3].
In metazoans, the period following fertilization is typically characterized by rapid and nearsynchronous mitotic divisions and cleavages that occur under conditions of minimal cellular differentiation. Cleavages typically depend on maternally supplied factors and zygotic genome transcription is constrained during this early period of development [4]. A conserved process of metazoan embryogenesis is the maternal-to-zygotic transition (MZT), which is characterized by two critical steps: 1) the elimination of a maternal set of mRNAs and proteins and; 2) the beginning of zygotic transcription, which leads to the zygotic genomic control of development [5]. In D. melanogaster, the majority of the first set of zygotic transcripts are regulated by zelda (zld) [6], a zinc finger TF with particular affinity for promoter regions containing TAGteam sites-heptamers constituted by CAGGTAG and its variants [7][8][9]. zld (Dm-zld, in D. melanogaster) binding sites have been identified in D. melanogaster embryos (cycles 8 to 14) by chromatin immunoprecipitation coupled with high-throughput sequencing (ChIP-Seq) [8]. Dm-zld regulates a large set of genes involved in important processes such as cytoskeleton organization, cellularization, germ band development, pattern formation, sex determination and miRNA biogenesis [6]. Dm-zld was also suggested to participate in larval wing disc development [10] and its overexpression during wing imaginal disc formation led to wing blisters in adults, an indicative of improper adherence of ventral and dorsal wing epithelia [11]. Nevertheless, while zld´s functions have been thoroughly investigated in D. melanogaster MZT, its roles in other organisms and biological processes remain elusive.
D. melanogaster displays a long germ type of embryonic development, during which all segments are simultaneously generated along the whole egg. In contrast, short germ insects generate anterior (e.g. head) segments early in development, while the remaining segments are patterned from the posterior region, the growth zone (GZ). Since short germ development is considered to be the ancestral mode of development, short-germ insects have been established as developmental model systems [12,13]. The short-germ red flour beetle Tribolium castaneum (Tc) was the first beetle species to have its genome completely sequenced [14]. T. castaneum displays a short life-cycle and is amenable to gene silencing via RNAi [15], gene overexpression [16], specific tissue expression [17] and fluorescence labeling during early development [18]. Several developmental processes that have been investigated in T. castaneum were lost or extensively modified in the D. melanogaster lineage, such as GZ formation [19], extensive extra-embryonic morphogenesis [20] and the formation of a morphologically complex head during embryogenesis [21]. Early development of T. castaneum is similar to most other insect groups, in which synchronous rounds of division are followed by nuclear migration to the egg cortex and cellularization, forming the so called uniform blastoderm [18,22]. Taken together, all the genetic and morphological information on T. castaneum early development, along with the established techniques mentioned above, make this beetle an ideal model to understand the evolution of zld´s function during insect development.
In the present work, we provide the first comprehensive analysis of zld orthologs across a wide range of species. Further, we provide functional analysis of a zld ortholog in a non-Drosophillid insect, the short germ beetle T. castaneum (thenceforth referred as Tc-zld). Among our main results are: 1) The identification of some previously overlooked conserved domains in Zld; 2) An inference of the evolutionary origin of Zld, based on phyletic analysis in various hexapods and crustaceans; 3) The identification of a conserved set of 141 putative Tc-zld targets (i.e. genes with upstream TAGteam sequences), enriched in TFs, whose homologs that have been demonstrated to be zld targets in D. melanogaster; 4) Identification of key roles played by Tc-zld during the MZT; 5) Identification and experimental validation of two new biological roles of zld in T. castaneum: segment generation from the posterior GZ during embryogenesis and postembryonic imaginal disc development; 6) Demonstration that zld is also involved in GZ patterning in the hemipteran Rhodnius prolixus, supporting a conserved role in the GZ of a distant short-germ species. Altogether, our results unveil new roles of zld as a pleiotropic TF acting in various developmental processes across distantly-related insects.

Identification of conserved domains and tracing the origin of zld
While previous studies reported that zld is involved in the MZT of D. melanogaster [6], its evolutionary history remains unclear. We investigated the phyletic distribution of zld and found a single ortholog in all the inspected insect genomes, including the beetle T. castaneum (Fig 1; S1 Table), indicating that insects are sensitive to increased copy number of this important regulator, which is interesting considering that different TF families are particularly prone to expansions across insect lineages [23]. Further, zld homologs were also found in some (but not all) collembolan and crustacean genomes ( Fig 1B; S1 Table). The canonical domain architecture of Zld has been reported as comprising a JAZ zinc finger (Pfam: zf-C2H2_JAZ) domain and a Cterminal cluster of four DNA binding C2H2 zinc finger domains (zf-C2H2) [6,11,24]. However, we performed a sensitive and detailed analysis of Zld proteins from multiple species and found other notable conserved protein domains and structural features. Firstly, we found two additional zf-C2H2 domains, N-terminal to the JAZ domain ( Fig 1A). The most N-terminal zf-C2H2 is absent or partially eroded (i.e. without the conserved cysteines and histidines) in some species (Fig 1B), including D. melanogaster. We define this N-terminal zf as ZF-Novel, since it has not been reported in previous studies. This observation was confirmed by inspecting Dm-zld alternative splicing isoforms. On the other hand, the second N-terminal zf-C2H2 domain is conserved in virtually all extant insects, but absent or degenerate in the other Pancrustaceans (e.g. Daphnia magna; partially conserved, with one lost cysteine) (Fig 1). Further, between this second zf-C2H2 and the JAZ domain, there is a strikingly conserved acidic patch that is characterized by an absolutely conserved motif of the form [DE]I[LW]DLD (Fig 1), which is predicted to adopt an extended conformation (using the JPRED software [25]) amidst surrounding disordered regions. A related conserved acidic motif was also found in the chordate protein CECR2, C-terminal to the DDT and WHIM motifs, which constitute a helical Asterisks and x marks represent presence and absence, respectively. Question marks denote that the feature is either partially preserved or could be flagged as absent due to sequencing or assembly errors (e.g. wrong start codons). The cladogram was organized according to a previously reported phylogenetic study [79]. Two outgroups without zld orthologs were also included.
https://doi.org/10.1371/journal.pgen.1006868.g001 zld in the beetle Tribolium domain involved in setting the nucleosome spacing in conjunction with the ISWI ATPase during chromatin remodeling [26]. An analogous conserved acidic patch is also seen in the HUN domain which functions as a histones chaperone [27]. Taken together, these observations raise the possibility that this Zld acidic region interacts with positively charged chromatin proteins such as histones. A largely disordered region, between the JAZ finger and the cluster of 4 widely-conserved zf-C2H2 domains, has been shown to be important for Dm-zld transactivation in vitro [24]. Our analysis of evolutionary constraints on the protein sequence also revealed a motif of the form hP [IVM]SxHHHPxRD, which appears to be under selection for retention despite the strong divergence in this region. Hence, it is possible that it specifically plays a role in transactivation. Further, between the two N-terminal zf-C2H2 domains, there is a highly conserved RYHPY motif, which could be involved in nuclear localization (Fig 1A), as predicted for other TFs [28]. Given the conservation of these additional domains, we hypothesize that they also play important roles in Zld functions that were previously attributed exclusively to the C-terminal domains.
Aiming to elucidate the origins of zld, we performed extensive sequence searches and were unable to find homologs with the conserved canonical domain architecture outside of Pancrustacea, indicating that zld is an innovation of this lineage. We detected clear zld homologs across insects, including the termite Zootermopsis nevadensis (Order Isoptera), the scarce chaser Ladona fulva (Order Odonata), the mayfly Ephemera danica (Order Ephemeroptera) and Machilis hrabei (Order Archaeognatha). zld homologs were also found in crustaceans belonging to different classes (Daphnia magna, Hyalella azteca and Eurytemora affinis), as well as in the collembola Folsomia candida. Curiously, we found no zld in the genome of Orchesella cincta (collembola) and Daphnia pulex (crustacean), suggesting either that it is not absolutely conserved outside of insects or missing due to incompleteness of the deposited genomes. Specifically, we carefully searched the genomes of other non-insect arthropods, including chelicerates (e.g. the tick Ixodes scapularis and the spider Parasteatoda tepidariorum), and found no proteins with the canonical Zld domain organization ( Fig 1B). General searches on Genbank against non-Pancrustacea arthropod proteins also returned no Zld orthologs. Although BLAST searches with Zld proteins against the nr and refseq databases recovered a number of significant hits in several distant eukaryotes, the similarity is almost always restricted to the Cterminal cluster of four Zf-C2H2 domains, which are very common across several TF families (e.g. glass, earmuff/fez, senseless/gfi-1 and jim). Taken together, our results support the early emergence of zld in the Pancrustacea lineage, with subsequent losses in particular species. Importantly, all the insect genomes we inspected have exactly one zld gene, indicating that this gene became essential in hexapods.

Tc-zld is a master regulator of signaling genes and other transcription factors
Previous studies in D. melanogaster using microarrays and ChIP-Seq revealed that Dm-zld regulates the transcription of hundreds of genes during early embryogenesis [6,9,29]. In D. melanogaster, enhancers bound by Dm-Zld are characterized by a consensus sequence CAGGTAG (i.e. TAGteam sequence), which is overrepresented in early zygotic activated genes, including TFs involved in AP and DV patterning [7,8]. Since the TAGteam motif identified in D. melanogaster is conserved in A. aegypti [30], we investigated whether we could predict Tc-zld targets by detecting TAGteam motifs in the upstream regions of T. castaneum genes.
Firstly, an ab initio approach using DREME [31] was employed to analyze 2kb upstream regions of all T. castaneum protein-coding genes. This analysis uncovered a motif (i.e. GTAGGTAY) that is nearly identical to the TAGteam motif (Fig 2A). We used the D. melanogaster genome and experimental data [6,9,29] to validate our approach and found a significant overlap between experimental and predicted zld targets in D. melanogaster ( Fig 2B). The putative T. castaneum motif was then used to screen the T. castaneum genome, resulting in the identification of 3,250 putative zld targets, representing~19% of the T. castaneum genome ( Fig 2C, S2 Table). Comparison of the putative Tc-zld targets with 1,087 genes regulated by Zld during D. melanogaster MZT [8,29] allowed the identification of 141 D. melanogaster genes for which one-to-one orthologs figured among the putative Tc-zld targets (hypergeometric distribution, P < 4.5x10 -4 ; Fig 2C, S3 Table). Functional analysis of this gene set using DAVID [32] uncovered the enrichment of important categories, including a highly significant cluster of 26 homeobox TFs (S4 Table) and other significant clusters comprising genes involved in regionalization and segment specification, imaginal disc formation and metamorphosis (S4 Table). Interestingly, this gene set included multiple developmental regulators such as anteroposterior (AP), gap, pair-rule, homeotic and dorsoventral (DV) genes (S2 Table). Since several of these genes are involved in early developmental processes, we focused our initial analysis of Tc-zld´s function during embryogenesis.

Tc-zld is maternally provided and dynamically expressed during embryogenesis
Since zld is maternally expressed in the germ line in D. melanogaster, we compared its transcription in T. castaneum female ovaries and carcass by qRT-PCR. T. castaneum early development starts with synchronous divisions during the first three hours of embryogenesis (at 30˚C), followed by nuclear migration to the egg cortex and membrane segregation of nuclei into separate cells, approximately 7-8 hours after egg lay [18,22]. We found that Tc-zld is highly expressed in the ovaries (Fig 3A), supporting the transcription of Tc-zld in the germ line. The abundance of Tc-zld transcripts is also higher in the first three hours of development than in the next two 3-hour periods (i.e. 3-6 and 6-9 hours) (Fig 3B), suggesting that Tc-zld mRNA is maternally provided and degraded after the first 3 hours of development. An antibody against the transcriptionally active form of RNA polymerase II, previously used in other ecdysozoan species [33,34], showed that zygotic transcription in T. castaneum begins between three and six hours of development, shortly after the nuclei have reached the periphery (Fig 4).  [7] obtained by DREME and putative T. castaneum zld DREME motif. (B) Venn diagram of the D. melanogaster ChIP-Seq MZT regulated genes (green) from [8] and D. melanogaster genes predicted by FIMO analysis with the putative DREME motif (yellow). (C) Venn diagram of the D. melanogaster Dm-Zld MZT targets (green) [8] and D. melanogaster one-to-one orthologs of putative Tc-zld targets.
https://doi.org/10.1371/journal.pgen.1006868.g002 In situ hybridization confirmed maternal ubiquitous expression of Tc-zld in the first three hours of development (Fig 3C and 3D) and showed a progressive confinement of Tc-zld mRNA to the posterior region of the egg, where the germ rudiment will be formed (Fig 3E and  3F). Between 6 and 9 hours of development, Tc-zld expression is observed at the embryonic tissue ( Fig 3G), with higher levels at the posterior region (Fig 3H), where the GZ will generate new segments [19,35]. In addition, Tc-zld is also expressed in the ventral serosa, during the serosal window closure (Fig 3J). Later in development, Tc-zld expression is still detected at the GZ (Fig 3I), at the head lobes and a single gnathal segment ( Fig 3K) and, subsequently, in the nervous system ( Fig 3L). Although zld is maternally provided, and zygotically expressed at the neural progenitors of both D. melanogaster and T. castaneum [6,11], biased posterior expression is a feature so far described only in short-germ insects like the latter.

Parental injection of Tc-zld dsRNA reduces number of eggs laid and impairs embryogenesis
It has been shown that injection of dsRNA in T. castaneum females leads to reduced expression of a given gene in the females and their offspring, in a phenomenon called parental RNAi (pRNAi) effect [15]. We injected zld dsRNA in females and analyzed Tc-zld transcriptional levels during embryogenesis by qRT-PCR. After zld pRNAi, Tc-zld mRNA levels were reduced in the first two weeks of egg laying, severely impairing larval hatching (Fig 5A-5C). Importantly, identical knockdown phenotypes during embryogenesis were obtained by using a second, in F. (G) Shortly before posterior invagination, zld expression occurs along the whole germ rudiment, with higher levels at the posterior region where growth zone will form (6-9 hours). (H) Shortly before the beginning of germ band elongation zld expression is highly expressed at the GZ, although lower levels of mRNA can be observed at the along the whole embryonic region. (I) Approx. 13 hours after egg lay (AEL), during the beginning of germ band extension, Tc-zld is expressed at the posterior GZ (arrow). (J) An embryo slightly older than the one in H, highlighting Tc-zld expression at ventral serosal cells during serosal window closure. (K) Approx. 18-21hours AEL zld is still expressed at the posterior most region (arrow), head lobes (arrowheads) and a single gnathal segment (asterisk). (L) Approx. 39-48 hours AEL, during early dorsal closure, expression is observed at the nervous system (arrow).
https://doi.org/10.1371/journal.pgen.1006868.g003 Western-blots of embryonic extracts from 0-1, 0-3, 3-6 and 6-9 hours after oviposition using an antibody against the transcriptionally active form of RNA pol II, as previously described [33]. (D) or an antibody against the αtubulin protein as a loading control (H). (E,F,G) Immunostaining showing nuclear RNA pol II at 3-6 hours (F) and 6-9 hours (G), but not at 0-1 hour (E) after oviposition. Coupling RNA pol II staining with nuclear DAPI staining shows that T. castaneum zygotic transcription starts between 3-6 hours after egg laying, when the energids reaches the periphery (Fig 2A-C-DAPI and E-G-RNA pol II). Similar results were obtained by western-blots using the same antibody (D).  Further, morphological analyses showed that cellularization was severely disrupted in over 50% of the zld dsRNA embryos (Fig 5D), similarly to what was previously reported in D. melanogaster zld mutants [6,11]. The remaining zld pRNAi embryos were not severely affected during cellularization and developed beyond that stage. Finally, we also found that some putative conserved target genes were down-regulated in embryos after zld pRNAi, such as the early zygotic genes involved in AP patterning, the serosal gene Tc-zerknullt [36] (Fig 5E and 5F) and the gap gene milli-pattes [37] (Fig 5G and 5H).
Changes in the spatial distribution of transcripts from predicted dorsoventral target genes were also observed after Tc-zld RNAi. In wild-type (WT) embryos, the TF Dorsal forms a dynamic transient gradient, which activates Tc-cactus (Tc-cact) and Tc-short-gastrulation (Tcsog) at the ventral region [38,39]. After Tc-zld pRNAi, the expression of Tc-cact and Tc-sog is observed in two lateral domains, in contrast to the expression at the single ventral domain in dsneo RNAi embryos (Fig 5I, 5J, 5K and 5L). These results suggest that Tc-zld is required for proper activity of Dorsal at the ventral-most region of the embryo.
As discussed above, Tc-zld is highly expressed at the posterior region of the embryo (Fig 3) and likely associated with segmentation and regionalization (S4 Table). Further, several putative Tc-zld targets are involved in posterior segmentation, such as caudal (Cdx), even-skipped (Eve) and several Hox genes (e.g. Ultrabithorax, Abdominal-A and Abdominal-B). Tc-eve, for example, is essential for the establishment of a genetic circuit required for posterior segmentation [40]. Interestingly, Tc-zld pRNAi embryos showed a continuous Tc-eve expression domain instead of the typical stripe patterning required for WT segmentation (Fig 5O and  5P).
Elegant studies on T. castaneum GZ patterning showed that cell proliferation is not essential for segment generation, which rather occurs by coordinated cell movement and intercalation [19,41]. We then evaluated whether Tc-zld regulates genes involved in cell intercalation, such as Toll2, Toll6 and Toll8 [42,43]. Interestingly, Toll7 (TC004474) and Toll8 (Tollo:TC004898) are among the common zld targets conserved between D. melanogaster and T. castaneum (S3 Table). Since Tc-Toll7 is expressed during early segmentation [43], we compared its expression in dsneo and Tc-zld RNAi embryos. While anterior expression of Tc-Toll7 is apparently unaffected, the striped expression at the posterior region is lost when Tc-zld expression is reduced (Fig 5M and 5N). Further support for the loss of posterior segmentation after Tc-zld RNAi is also provided by the analysis of expression of the segment-polarity gene, Tc-gooseberry (Tcgsb) (Fig 5Q and 5R). In summary, Tc-zld regulates the expression of several genes that are critical for early AP (zen, mlpt) and DV (sog, cact) patterning and, in a second phase, genes required for posterior elongation (e.g. Toll7, gsb).
Tc-zld plays specific roles in the posterior GZ While pRNAi diminishes maternal and zygotic expression of Tc-zld in T. castaneum [15], embryonic dsRNA injections (eRNAi) may affect only the zygotic component, since eggs can gap gene Tc-mlpt [37] in control (G) and zld RNAi (H). Tc-mlpt loses its anterior expression domain in zld RNAi when compared to control. (I,J) The expression of the dorsoventral gene Tc-cactus [38]  be injected after the MZT [38,44]. To investigate if Tc-zld is specifically required for embryonic posterior patterning, we injected Tc-zld dsRNA in transgenic embryos expressing nuclear GFP (nGFP), as previously described [18,19]. Embryonic injections of Tc-zld dsRNA after the MZT (see Methods for details) impaired segment generation from the GZ, while dsneoinjected embryos developed like WT ones (Fig 5S and 5T, S1 and S2 Movies). In addition, expression of the predicted target Tc-eve, a key TF involved in GZ patterning [35], has been largely down-regulated upon Tc-zld eRNAi, as previously observed for zld pRNAi (Fig 5O and  5P). In summary, our results imply that Tc-zld is involved not only in the MZT, early patterning and nervous system formation, as described for D. melanogaster [6], but also play roles in segment generation from the GZ, a structure found only in embryos of short-germ insects like T. castaneum.

Functional analysis in the hemipteran Rhodnius prolixus shows a conserved role of zld in short-germ insect posterior region
Although we demonstrated the involvement of zld in T. castaneum GZ, the conservation of this regulatory mechanism in other species remained unclear. Thus, we sought to analyze the functions of the zld ortholog in the hemipteran R. prolixus (Rp; Rp-zld gene), which is a hemimetabolous insect and lacks the complete metamorphosis present in holometabolous species such as T. castaneum [45,46]. Rp-zld knockdown via pRNAi resulted in two types of embryonic phenotypes: 1) severe defects in gastrulation and lack of any appendage development; 2) embryos that developed only the anterior-most embryonic regions comprising the head, gnathal and thoracic segments (Fig 6). These results support a role for zld in posterior elongation that was likely present in the common ancestor of hemimetabolous true bugs and holometabolous insects, if not earlier. To our knowledge this is the first direct description of zld function in insects other than D. melanogaster. A recent report showed zld as maternally transcribed in the hymenoptera Apis mellifera (also a long-germ insect), while zygotic transcripts are concentrated at the central region of the embryo during blastodermal stages [47]. Further, during the preparation of this manuscript, zld has been also proposed to be required for MZT and cellularization in the wasp Nasonia vitripennis [48].
As the first step towards the characterization of the post-embryonic role of zld, we analyzed Tc-zld expression by qRT-PCR in larvae of third (L3), fifth (L5) and seventh (L7) stages, and first pupal stage (P1) (Fig 7A-7D). Interestingly, Tc-zld expression increases during successive larval stages and sharply decreases after pupal metamorphosis (Fig 7E). This suggests that Tczld might be required for late larval stages, which take place during larvae-pupae metamorphosis, such as growth and patterning of structures derived from imaginal discs in D. melanogaster (e.g. antennae, legs, fore-and hindwings) [49,50]. Further, we found that three out of five predicted Tc-zld targets, namely Dll, Wg and Lim-1, also displayed an increase in expression during late larval and pupal development (Fig 7E).
To investigate zld´s post-embryonic roles, we injected two non-overlapping Tc-zld dsRNA constructs into early (L3) and late larval (L6) stages, as previously described [50]. qRT-PCR confirmed that Tc-zld was down-regulated after dsRNA injection (Fig 8A). Tc-zld dsRNA injections at early larval stages (L3) led to over 50% lethality during pupal stages. Atypical adult pigmentation in the pupal head and reduction in the wing size were observed after early Tc-zld dsRNA injection, suggesting that proper Tc-zld expression is required for metamorphosis and wing growth (Fig 8B and 8C). Interestingly, Tc-zld dsRNA injections at late larval stage (L6) displayed a different phenotype when compared with early larval dsRNA injections (L3). L6 larvae injected with dszld and dsneo reached adulthood at comparable rates (Fig 8D). Specifically, Tc-zld dsRNA adults showed a series of morphological alterations in tissues undergoing extensive morphological changes during metamorphosis, such as fore-and hindwings, antennae and legs (Figs 8E, 8F, 9 and 10).
The most visible effect of Tc-zld dsRNA beetles was a failure of the forewings (elytra) to enclose the hindwings, leading to the exposure of the dorsal abdomen (Fig 8E and 8F). Elytra, which are highly modified beetle forewings, have been proposed to be an important beetle innovation, being required for protection against mechanical stresses, dehydration, and predation [51]. In line with this hypothesis, Tc-zld dsRNA adults with exposed abdomens started to die a few days after metamorphosis, probably due to dehydration.
Next, we performed a detailed morphological analysis to investigate if patterning defects resulting from Tc-zld dsRNA occurred in the sclerotized elytra (forewing) or in the hindwing (Fig 9A-9D). This analysis showed that the parallel vein pattern of the elytra (Fig 9A) is disrupted after Tc-zld dsRNA in comparison to the control (Fig 9A and 9B). Nevertheless, hindwings of Tc-zld dsRNA beetles showed no signs of abnormal venation (Fig 9C and 9D). To address if fore or hindwing shapes have changed upon Tc-zld dsRNA knockdown, we applied the recently developed Proximo-Distal (PD) index, a morphometric analysis consisting of the measurement of the wing length and width at two positions [52]. While the PD index values of forewings (elytra) of Tc-zld and controls were similar, a slight but significant increase in the hindwing PD index was observed (Fig 9E). Neither fore-nor hindwings length were altered upon Tc-zld knock down (Fig 9F). In conclusion, Tc-zld is required for proper venation pattern, but not shape, of the elytra. Previous analysis of zld expression in D. melanogaster showed expression in wing imaginal discs, particularly where mitotically active cells are located [11]. Moreover, Dm-zld overexpression during wing imaginal disc formation leads to adult wing blisters or tissue loss [10,11]. Nevertheless, to our knowledge, our study provides the first direct evidence of zld´s role in insect wing formation.
Interestingly, all Tc-zld dsRNA beetles that showed this 'opened wing' phenotypes ( Fig 8F) also displayed defects in the legs and antennae. Antennae and legs share similar developmental GRNs, involving the so called serial homologs. Distaless (Dll), one of the putative Tc-zld targets, is essential for appendage segmentation [53]. Tc-Dll is also expressed during late embryogenesis on the distal part of the leg and, as its name suggests, disruption of its function leads to the absence of distal leg and antennae segmentation [54]. Interestingly, we found that Tc-zld RNAi resulted in a significant decrease of Tc-Dll mRNA levels (Fig 10A), indicating that this gene is indeed downstream of Tc-zld.
Insect antennae possess three primary segments: scape, pedicel and flagellum. In T. castaneum, the adult antennae display eleven segments, out of which nine form the flagellum. The three most distal flagellar segments are enlarged and form the club, while the six intermediate flagella are called the funicle (Fig 10B) [55]. In mild phenotypes, dszld caused a joint malformation in distalmost flagellar segments resulting in a fusion of the club (Fig 10C). On the other hand, strong Tc-zld dsRNA phenotypes resulted in fusion of the scape and pedicel, leading to severe loss of flagellar joints and formation of a single truncated segment (Fig 10D). In contrast, no differences in scape and pedicel were observed.
T. castaneum legs originate during late embryogenesis and can be recognized as a small outgrowth of the body wall, the limb bud. In the adult stage, there are three pairs of segmented legs with six segments: coxa, femur, trochanter, tibia, tarsus and pretarsus [56]. In T. castaneum the tarsus is subdivided in smaller segments (i.e. the tarsomeres), five in prothoracic and mesothoracic legs and four in metathoracic legs. Tarsal segmentation occurs during beetle . The statistical analysis was carried out by unpaired t-test assuming unequal variances (asterisks refer to P<0.0001 while ns stands for "no significance"), indicating that Tc-zelda RNAi affects hindwing shape. (F) Tc-zld RNAi elytra and hindwing do not show statistically significant differences in length when compared to their respective dsneo RNAi controls.
https://doi.org/10.1371/journal.pgen.1006868.g009 metamorphosis and this subdivision of the tarsus evolved in the common ancestor of insects, since the tarsus is not subdivided in non-insect hexapods [57]. After Tc-zld RNAi, tarsal segments were absent or fused, resulting in leg shortening; in dsneo insects, legs were identical to that of WT animals (Fig 10E and 10F). This indicates that some of the Tc-zld targets are involved in segment development or joint formation. Interestingly, a large domain of Dll expression is observed in beetle leg tarsus [58], and Tc-Dll knockdown in beetles also generated legs with tarsomere deletion [58,59]. Recently, Dm-zld has been demonstrated to be a key factor in the establishment of the early chromatin architecture, particularly for the formation of topologically associating domain (TAD) boundaries [60]. During D. melanogaster MZT, zld resembles a pioneer TF marking the chromatin of earliest expressed genes [8,61]. Dm-zld also increases chromatin accessibility of the most important TFs involved in DV and AP patterning (Dorsal and Bicoid, respectively) [62][63][64][65]. Further, the addition or removal of Dm-zld binding sites influences the timing of activation of Dorsal early zygotic targets in D. melanogaster [62,63], suggesting that Zld acts as a developmental timer.
Our results in T. castaneum showed an extensive and ubiquitous maternal contribution of Tc-zld mRNAs, followed by zygotic Tc-zld expression along the embryonic rudiment, particularly at the posterior region (Fig 3). It is possible that Tc-zld mediates a progressive anteroposterior opening of the chromatin in T. castaneum, shortly before these posterior GZ cells undergo convergent extension movements required for germ band elongation [19,41]. Loss of Tc-zld expression might lead to lack of convergent extension due to loss of Toll7 and eve expression and, ultimately, segmentation failure (Fig 5). Since zld is also important for the development of the posterior region of the hemimetabolous insect R. prolixus (Fig 6), zld´s role in posterior region dates back at least to the last common ancestor of Paraneoptera.
Tc-zld expression was also observed at the ventral serosa ( Fig 3J) during embryonic stages. While we cannot rule out that Tc-zld plays a role during normal serosal development, lack of segmentation after Tc-zld eRNAi cannot be explained by loss of Tc-zld expression in this tissue, since embryos without serosa, e.g. Tc-zen1 RNAi, do not display phenotypic defects under normal developmental conditions [36,66]. Future studies are warranted to determine if Tc-zld plays a specific role in the beetle serosa, particularly by regulating some of its predicted target genes involved in immune responses (e.g. Toll, cactus).
Tc-zld might also play a similar role in the regulation of leg and antenna segmentation during metamorphosis (Fig 10). The number of tarsomere segments in the leg and intermediate funicle are reduced after Tc-zld dsRNA injection, suggesting that the post-embryonic role of Tc-zld in appendage segmentation might also involve the regulation of complex GRNs, as observed in the posterior embryonic region.
The evolutionary success of hexapods is attributed to a combination of features: their segmented body plan and jointed appendages, which were inherited from their arthropod ancestor and; wings and holometaboly, two features that arose later in insect evolution [67]. It is interesting to note that the specific Pancrustacea gene zld is required for most of these processes, such as embryonic segment formation, wing (elytra) patterning, and appendage (antennae and leg) formation during beetle development. Like several other zinc finger proteins that show tremendous lineage-specific diversity in eukaryotes, zld appears to have specifically arisen within Pancrustacea and risen to play an important role as a "master TF". Hence, future studies focusing on how this TF was integrated to the conserved backdrop of developmental TFs and existing GRNs of arthropods would be of great interest.

Beetle rearing
T. castaneum beetles were cultivated in whole-wheat flour. For sterilization, the flour was kept for 24h at -20˚C and another 24h at 50˚C. The beetles were maintained inside plastic boxes of approximately 15x15cm with humidity between 40-80%.

Bioinformatic analyses
Protein sequence analyses were performed using BLAST [68] and PFAM searches using HMMER3 [69] against proteins available on Genbank [70] and Vectorbase [71]. Dm-Zld (NP_608356) was used as initial query for sequence searches. Genomic data from the following genomes were obtained from the Baylor College of Medicine Human Genome Sequencing Center: Eurytemora affinis, Hyalella azteca, Blattella germanica, Catajapyx aquilonaris, Machilis hrabei, Libellula (Ladona) fulva and Ephemera danica. BLAST results and domain architectures were manually inspected. In silico searches in the T. castaneum genome for over-represented motifs were performed using the DREME software [72], part of the MEME suite software toolkit [73]. The motif with highest identity with Dm-Zld binding site (CAGGTAY) was compared to previously described motifs in the FlyFactorSurvey [74] database with TOMTOM. Further, this motif was used to scan 2000 bp upstream of all predicted T. castaneum genes with FIMO [14,75,76]. Upstream regions and orthology information were retrieved from Ensembl Metazoa Biomart (http://metazoa.ensembl.org/biomart/) [77].

dsRNA synthesis, parental RNAi and embryonic phenotypic analysis
Two non-overlapping PCR fragments containing T7 promoter initiation sites at both ends were used as templates for Ambion T7 Megascript Kit (Cat. No. AM1334) following the manufacturer instructions (for details see S1 Fig). The amount and integrity of the dsRNA samples were measured by spectrophotometry and agarose gel electrophoresis, respectively. For parental RNAi (pRNAi) analysis,~0.5μl of dsRNA were injected from a solution containing 1μg/μl of dsRNA into adult female beetles [36]. Eggs were collected for four egg lays (2 day each) and zld´s down regulation was estimated by quantitative Real Time PCR (see below).

T. castaneum embryonic dsRNA injections
Egg injections were performed as previously described [38,44]. Briefly, for the analysis of Tczld zygotic role, embryos containing nuclear-localized green fluorescent protein (GFP) were collected for one hour and let to develop for an additional three hours (30˚C) [19,35]. After this period, twenty embryos were dechorionated with bleach (2% solution), aligned onto a glass slide and covered with Halocarbon oil 700 (Sigma). Embryos were immediately microinjected at the anterior region with zld or neo dsRNA at 1 μg/μL concentration with the help of a Nanoinject II instrument (Drummond Scientific Company). After injection, a single nGFP embryo was photographed every five minutes during the following 16 hours (25˚C) in a Leica DMI4000 inverted microscope using a GFP filter. Single photographs were used to generate a movie using Windows Movie Maker (S1 and S2 Movies). Phenotypes of all injected embryos (neo or zld dsRNA) were scored at the end of the experiment.

T. castaneum larval dsRNA injections
Larvae were injected with zld or neo dsRNA as previously described [50]. Knockdown phenotypes in pupae and adult beetles were generated by injection of zld or neo dsRNA solutions at a concentration of 1 μg/μL in the dorsal abdomen of individuals on third and sixth larval instars (n = 40). Following injection, larvae were reared in flour at 30˚C and collected periodically for RNA extraction and phenotype annotation. Adult beetles were then fixed in ethanol 95% overnight for further morphological analysis. Immunostainings have been performed as previously described [38].

Morphological analysis of imaginal disc derived tissues
Antennae, legs, elytra and wings were dissected using forceps and placed in a petri dish for observation. Phenotypic analyses and documentation were performed under a Leica stereoscope model M205.

PD index
The methods for wing and elytra measurement and PD index were performed according to described by [52]. Leica AF Lite software was used for the wing measurements. Image properties were adjusted in Adobe Photoshop CS4.

Quantitative real-time PCR
For experiments using embryos, total RNA was isolated from 100 mg of eggs collected from specific development stages (0-3, 3-6 and 6-9 hours after egg laying), ovary and carcass (whole beetle without ovary) using Trizol (Invitrogen), according to the manufacturer's instructions. Three independent biological replicates were used for each assay. First strand complementary DNA (cDNA) was synthesized from 2 μg of RNA using Superscript III reverse transcriptase (Invitrogen) and oligo(dT). The cDNA was used as template for real time qRT-PCR analysis using SYBR green based detection. qRT-PCR reactions were carried out in triplicate, and melting curves were examined to ensure single products. Results were quantified using the ''delta-delta Ct'' method and normalized to rps3 transcript levels, as previously described [78]. Primer sequences used during the study are provided at the supplemental data (S5 Table).
zld pRNAi in the hemiptera Rhodnius prolixus (Rp) zld cDNA sequence was initially identified by BLAST and included in Fig 1. Parental RNAi against Rp-zld was performed as previously described [46]. Supporting information S1 Fig. Schematic drawings of Tc-zelda non-overlapping dsRNA constructs. Tc-zld gene corresponds to the Beetlebase ID: TC014798. The following primer pairs were used: C1-Forward ggccgcggAGCGCATCTTCTCCCTATCA and C1-Reverse cccggggcGCCGTTTTGT CGTTTCTCAT. This primer pair amplifies leads to the amplification of a fragment of 528 bp (349-857). A second primer pair C2-Forward ggccgcggACGACGAGTACCGCTTGACT and C2-Reverse-cccggggcCTTACCACAGGTGTCGCAGA was also used and lead to similar results. This primer pair amplifies a fragment of 476 bp (1823-2279) covering the four zincfinger domains. The lowercase letters contain the primer sequence used as a template for a second PCR using universal primers which adds a T7 promoter at both sides of the PCR template as previously described [81]. (TIF) S1