The genomic and functional landscapes of developmental plasticity in the American cockroach

Abstract

Many cockroach species have adapted to urban environments, and some have been serious pests of public health in the tropics and subtropics. Here, we present the 3.38-Gb genome and a consensus gene set of the American cockroach, Periplaneta americana. We report insights from both genomic and functional investigations into the underlying basis of its adaptation to urban environments and developmental plasticity. In comparison with other insects, expansions of gene families in P. americana exist for most core gene families likely associated with environmental adaptation, such as chemoreception and detoxification. Multiple pathways regulating metamorphic development are well conserved, and RNAi experiments inform on key roles of 20-hydroxyecdysone, juvenile hormone, insulin, and decapentaplegic signals in regulating plasticity. Our analyses reveal a high level of sequence identity in genes between the American cockroach and two termite species, advancing it as a valuable model to study the evolutionary relationships between cockroaches and termites.

Fig. 1

Orthology and genome evolution of Blattodea. a Orthology assignment of four blattodean species and eight other representative insect species (Supplementary Table 6). Bars are subdivided to represent different types of orthology clusters as indicated. Universal groups represent common gene families across all species analyzed, but absence in at most one genome is tolerated; of them, “single-copy” tolerates absence or duplication in a single genome, while “multiple-copy” indicates other universal genes. “Blattodea only” indicates unique presence to the blattodean lineage and in at least three genomes. “Specific” groups indicate specific presence or duplication in only one species. “Homology” indicates partial homology detected with E < 10⁻⁵ but no orthology grouped. Remaining orthologs were assigned to “Present at half of species” or “Patchy”, depending on whether presence is in at least six genomes or not. The phylogeny was calculated using maximum-likelihood analyses of a concatenated alignment of 538 exactly single-copy proteins. The tree was rooted using the sister clades of Blattodea and Orthoptera species as described previously. Bootstrap values based on 100 replicates are equal to 100 for each node. b Distribution of sequence identities. The boxplots delineate the interval between the first and the third quartiles of the identity distribution between P. americana and one indicated blattodean species. Notch indicates the median value. Amino acid sequence identity was calculated based on multiple alignment of each universal single-copy ortholog (5911 in total) in four blattodean species. The red dashed line indicates the median value within Periplaneta species (96%). c Comparison of sequence identities on each ortholog cluster. A total of 7640 1:1:1 orthologs were identified among P. americana, B. germanica, and Z. nevadensis. Each dot represents such an ortholog. Value on x-axis indicates the sequence identity between P. americana (Pame) and B. germanica (Bger), while value on y-axis indicates that between P. americana and Z. nevadensis (Znev). Red line indicates the smoother of a locally weighted regression, with the coefficient of determination as 0.99999. Dashed line indicates positions where the sequence identity between P. americana and B. germanica is identical to that between P. americana and Z. nevadensis

Fig. 2

Gene families involved in chemoreception and detoxification in P. americana and other blattodean species. a Counts of chemosensory- and detoxification-related gene families in the genomes of three blattodean species and Drosophila melanogaster. OR olfactory receptor, GR gustatory receptor, IR ionotropic receptor, OBP odorant-binding protein, P450 cytochromes P450 (CYP), CCE carboxyl/choline esterases, GST glutathione S-transferase, ABC ATP-binding cassette. Of these, three gene families with massive expansions in P. americana were selected for maximum-likelihood phylogenetic analysis, as shown in b–d, representing GRs, IRs, and P450, respectively. Phylogenetic relationships of other gene families are shown in Supplementary Figs. 3–7. iGluRs ionotropic glutamate receptors, NMDA N-methyl-d-aspartate receptors, Mito the mitochondrial clan

Fig. 3

Gene repertoire of the innate immune system and functional analyses of the Toll pathway in P. americana. a Representation of three main innate immune signaling pathways (Imd, Toll, and JAK-STAT, proposed in Drosophila)^– in P. americana. Genes from expanded gene families in P. americana are highlighted in red. Since all genes absent in P. americana were also found to be absent in other insect species, they were defined as Drosophila-specific components and outlined by gray dashed lines. 1, fas-associated DD protein; 2, death-related Ced-3/Nedd2-like protein; 3, effete; 4, TAK1-binding protein 2; 5, TGF-β-activated kinase 1; 6, kenny; 7, immune response deficient 5; 8, spätzle-processing enzyme; 9, myeloid differentiation primary response 88; 10, tube; 11, hopscotch; 12, socs36E; AMP, antimicrobial peptide; DOME, domeless; Dnr1, defense repressor 1; GNBP, Gram-negative binding protein; Gprk2, G protein-coupled receptor kinase 2; Grass, Gram-positive-specific serine protease; Iap2, inhibitor of apoptosis 2; IMD, immune deficiency; ModSP, modular serine protease; PGN, peptidoglycan; PGRP, peptidoglycan recognition protein; Pirk, poor IMD response upon knock in; Spz, spätzle; Spirit, serine protease immune response integrator; Stat92E, signal transducer and activator of transcription protein at 92E. Detailed information is additionally shown in Supplementary Tables 11–15. b Functional verification of genes in the Toll pathway against a classic Gram-positive bacterium, Staphylococcus aureus. Four major genes in the Toll pathway were knocked down by RNAi, compared to injection of control dsRNA (CK, a 92 bp non-coding sequence from the pSTBlue-1 vector). Corresponding mortality is shown upon S. aureus infection. Injection of cockroach saline solution (CSS, as a control) or S. aureus was performed 24 h after dsRNA injection. Error bars indicate standard deviation of three replicates. Bars labeled with different lowercase letters indicate significant difference between the two samples, with p < 0.05, one-way analysis of variance (ANOVA). c Functional relationship between the Toll pathway and expression of AMPs. Dorsal, a key component in the Toll pathway whose depletion caused the greatest mortality (shown in b), was knocked down by RNAi. Correspondingly, the relative expression of 11 antimicrobial peptide genes were measured and shown in c. Error bars represent s.d. of three replicates. Two-tailed Student's t-test: *p < 0.05; **p < 0.01

Fig. 4

Functional studies of pathways regulating development and reproduction in the American cockroach. a Regulation of molting by ecdysone receptor (EcR) and retinoid X receptor (RXR) genes in the 20E signal pathway. Mortality rates were checked after RNAi treatment of indicated numbers of individuals (n). CK, control dsRNA corresponds to that in Fig. 3. b Regulation of metamorphosis by methoprene-tolerant (Met) and kruppel homolog 1 (Kr-h1) in the juvenile hormone (JH) pathway. Adult proportion index was checked after the JH singling was disrupted. c Regulation of growth by the insulin signaling pathway genes. Insulin-like receptor (InR), phosphoinositide 3-kinase (PI3K), and target of rapamycin (TOR) genes were repeatedly depleted by dsRNA injection during 3 weeks. Three injections were given and cumulative growth rate of body weight (%) was calculated. Student’s t-test: ***p < 0.001. d Morphology change of ovary maturation during the first reproductive cycle in virgin females. Vitellogenin (Vg), double-sex (Dsx), and nine genes involved in the pathways of insulin, JH, and 20E were knocked down by RNAi. Gonadosomatic index and primary oocyte length were used to evaluate the ovary maturation degree. All the data were calculated as the mean value of three replicates. Error bars represent standard deviation. Two-tailed Student’s t-test: *p < 0.05; ***p < 0.001. Details of all involved pathways are shown in Supplementary Tables 16–18

Fig. 5

Leg regeneration ability and its regulation by the decapentaplegic (Dpp) pathway in P. americana. a Limb regeneration confirmation under different trauma severity indices. The diagram in the box represents a typical leg with indicated parts in which we performed the amputation experiment. Regeneration ability is shown in the bottom, according to the wild-type size. 4+, wild-type size; from 3+ to +, incomplete sizes; −, null. b Function of Dpp signaling pathway in leg regeneration. No leg regenerated (−) when Dpp and mothers-against-dpp (Mad) genes were knocked down by RNAi. CK, control dsRNA corresponds to that in Fig. 3