The present study focuses on true cockroaches (Blattaria), the major component of Dictyoptera. Taxa included in my analyses also cover other Dictyoptera, namely, termites (Isoptera) and mantises (Mantodea), and the living sister of Dictyoptera, namely Eukinolabia + Xenonomia, as suggested by transcriptome data (Misof et al. 2014, Evangelista et al. 2019, Wipfler et al. 2019). In order to reduce the difficulty and inaccuracy in alignment, more distant insect groups were not included in this study. All DNA sequences were collected from whole mitogenome sequencing data deposited in GenBank. Isolated fragments do not conform to the main idea of the present study because they leave missing cells in the alignments. A certain number of missing cells may be harmless to phylogenetic inferences (Wiens 2006, Wiens and Moen 2008), but there is no universal standard for different datasets and methods to avoid pitfalls. Using fully sequenced mitogenome (thus avoiding isolated fragments) can also prevent concatenation of sequences from different specimens, of which the genetic distance is unknown. This practice also automatically rules out specimens that are misidentified as conspecific. Fundamentally, concatenation of sequences from different specimens is an artefact and does not represent a natural organism. This artefact may cause unpredictable errors (pers. observ.); however, the influence of this issue seemingly has not yet been addressed in the literature.

The initial data pool comprises 169 mitogenomes, including all available cockroaches and selected other insects (Suppl. material 1). GenBank files were imported to PHYLOSUITE v1.2.1 (Zhang D et al. 2020), in which duplicates were removed. To save time and computation resources, only one species was kept in each genus, with the exception of genetically diverse and speciose genera (e.g., Cryptocercus, Allacta and Ischnoptera). The final taxon set comprises 95 species (Suppl. material 2: Table S1). The DNA data in the GenBank files were extracted using PHYLOSUITE. Most of the mitogenomes (68.4%) are from Bourguignon et al. (2018), and the remaining from Yamauchi et al. (2004), Cameron et al. (2005, 2012), Zhang YY et al. (2009), Kômoto et al. (2010, 2012), Chen (2012), Wang et al. (2014), Jeon and Park (2015), Tian et al. (2015), Cheng et al. (2016), Ye et al. (2016), Chen et al. (2017), Dumans et al. (2017), Ma et al. (2017), Gong et al. (2018), and Zhang LP et al. (2018).

The character set includes 13 protein coding genes while all RNA genes were excluded. Aligning RNA gene sequences is dependent on the prediction of secondary structure (Buckley et al. 2000, Stocsits et al. 2009), the accuracy of which clearly influences alignment and tree reconstruction (Letsch et al. 2010). This approach is unfeasible in the present study, because even predicting the structure of a small fragment of closely related species is hard and exhausting (e.g., Li et al. 2017). To avoid errors introduced by inaccurate alignment, I do not use them. Besides, RNA gene sequences only account for a minor proportion in the mitochondrial genome, therefore excluding them does not virtually reduce the size of dataset.

Sequences of protein coding genes were aligned using MUSCLE in MEGA7 with default settings of codon mode (Edgar 2004, Kumar et al. 2016). Uneven ends were manually trimmed. 43 sequences were spotted containing missing or poorly-sequenced portions, which were deleted or question-marked (Suppl. material 3: Table S2). Then, all sequences were aligned again. The aligned sequences were translated into amino acids to check accuracy. Alignments of 13 genes were concatenated using PHYLOSUITE. The final alignment is 10932 bases long.

ALIGROOVE 1.05 (Kück et al. 2014) was used to assess alignment ambiguity for each gene, and DAMBE 7.2.7 (Xia 2018) to calculate the substitution saturation per codon position per gene (i.e. 3 bases × 13 genes). Amino acid alignments (translated from the nucleotide alignments) instead of nucleotide alignments were imported to ALIGROOVE, because the alignments were based on codon model. Saturation was calculated under GTR model by default, or under F84 model when an overflow error occurs (the error itself is a sign of saturation).

ALIGROOVE suggests high ambiguity in ATP8 alignment, followed by ND6 (Suppl. material 4). Other genes exhibit low general ambiguity, but a few taxa with too many missing cells in alignments show high ambiguity. Even if several incomplete sequences can be excluded in the following analyses, ATP8, in general, is still too ambiguous. On the other hand, the dataset as a whole is not significantly affected by scattered ambiguous alignments (Suppl. material 4).

According to the 39 saturation plots (Suppl. material 5), ATP8, ND4L, ND6 and, as expected, the third codon positions are saturated. These three genes are among the less informative ones considered in Talavera and Vila (2011). Excluding ATP8, ND4L and ND6, data incompleteness was calculated by counting Ns and question-marks, both of which are regarded as missing cells. Eight taxa with missing cells greater than 1% were grouped into ‘BadSeq’: Anallacta methanoides (9.17% missing), Aposthonia borneensis (1.80%), Beybienkoa kurandanensis (1.18%), Eublaberus distanti (2.87%), Galiblatta cribrosa (3.64%), Megaloblatta sp. (2.92%), Metallyticus sp. (8.84%), and Platyzosteria sp. (2.19%).

Phylogenetic inferences were performed in MRBAYES 3.2.7 (Ronquist et al. 2012). Data were divided into two partitions: the first and the second bases of the codon. The third position of codon was excluded from all analyses. I did not use programs to select a ‘best-fit’ model, not only because this is unnecessary (Nascimento et al. 2017, Abadi et al. 2019), but also because the ‘best-fit’ is not necessarily the best or accurate (Gatesy 2007, Kelchner and Thomas 2007, Luo et al. 2010). Instead, I used the empirically universal model, GTR, with Gamma rates (+G). The manual and tutorials of MRBAYES (among others) recommend GTR+I+G as the universal model, but the invariable-sites model (+I) generates a strong correlation between the proportion of invariable sites and the gamma shape parameter, and becomes undesirable (Sullivan et al. 1999, Yang 2014). Using PHYLOBAYES with the CAT-GTR model implemented therein accounting for more exhaustive heterogeneities may improve the resolution and accuracy of phylogenetic inference (Lartillot and Philippe 2004, Lartillot et al. 2009, Moran et al. 2015), but this approach is currently impossible given the computation resources available to this study. Instead, the alignment ambiguity of the final dataset (see below) was reported, because the alignment ambiguity assessed by ALIGROOVE is also a measure of heterogeneous sequence divergence (Kück et al. 2014). I did not perform maximum likelihood analysis because the interpretation of bootstrapping (the assessment of the uncertainty of maximum likelihood estimates) is vague (Yang 2014), in contrast to posterior probabilities of Bayesian estimates (Huelsenbeck and Rannala 2004). Instead, the maximum-likelihood tree sampled from MCMC is reported. Each MRBAYES analysis involved two runs, each of which comprises four chains, running 1.5–2.5 million iterations depending on the difficulty of converging. Samples were taken once every 500 iterations and the initial 1%–5% (depending on the difficulty of converging) of samples were discarded (burn-in). I used TRACER 1.7.1 (Rambaut et al. 2018) to ensure sufficient effective samples (200 at least, 300–1000 in general).

The first analysis utilized all 13 genes. The results are only used for comparison with the second step analyses, to observe the influence of ATP8, ND4L and ND6.

The second step is to compare the trees inferred from three taxon sets. All analyses excluded ATP8, ND4L and ND6. (1) All-species analysis using complete taxon set. (2) Good-species analysis, excluding ‘BadSeq’. (3) Short-species analysis, excluding long-branched taxa detected from the all-species analysis. This step aims to detect the impact of incomplete data and long branch.

The third step analysis used only the ‘safe’ taxa. In this step, all taxa within ‘BadSeq’ were excluded even if they do not virtually affect the topology of other species. It is learnt from experience that more missing or poorly-sequenced bases imply more potential errors in the superficially intact data. Potential pitfalls of incompleteness (e.g. erroneous positions of these taxa per se) violate the main idea of this study. Long-branched taxa with low support are also to be excluded. In the present study, they are Aposthonia borneensis, Aposthonia japonica and Nocticola sp.. ‘Safe’ taxa comprise 85 species (Suppl. material 2: Table S1).

The fourth, also the final, step yields the phylogeny that is regarded as the formal result. Prior to MRBAYES, sequences of ‘safe’ taxa were re-aligned and concatenated. This new, 9912-base-long alignment, as final dataset, was also imported to ALIGROOVE to assess alignment ambiguity. This 85-species dataset is less ambiguous than the original one (Suppl. material 4). The resulted tree was used as the fixed topology in dating analyses.

As calibrations, only two fossils fulfill the criteria of Parham et al. (2011) and are suitable for the present study. The earliest known termite Valditermes brenanae calibrates the split between Cryptocercidae and Isoptera (minimum age 130.3 Ma, see Wolfe et al. 2016), as in Misof et al. (2014) and Evangelista et al. (2019). The earliest known blattelline cockroach Piniblattella yixianensis Gao et al., 2018 (Gao et al. 2018) calibrates the split between Blattellinae and Blaberidae. Piniblattella yixianensis is used for calibration for the first time, setting a minimum age as 120.9 Ma (see Discussion). To observe the effect of this new calibration, I also ran an analysis without this fossil (i.e. only calibrated by V. brenanae) to compare with the two-fossil analysis. The minimum bounds of calibrated nodes were set to the minimum age of fossils. The minimum root age was set to 130.3 Ma, the age of the older calibration fossil. All maximum bounds were set to 412 Ma, the oldest age of Rhynie Chert, as justified in Evangelista et al. (2019).

Some studies used the so-called ‘roachoids’ (Eoblattodea, see Li 2019) to calibrate Dictyoptera+sister (which is the root herein), based on the hypothesis that those ambiguous fossils are stem members of Dictyoptera (Legendre et al. 2015, Tong et al. 2015, Bourguignon et al. 2018, Evangelista et al. 2019). To detect the impact of such fossils, I performed another dating analysis with the earliest ‘roachoid’, namely Qilianiblatta namurensis Zhang et al., 2012 (Zhang ZJ et al. 2012, Guo et al. 2013), thus three fossil calibrations were used in this analysis. The radioisotopic age of the Q. namurensis-bearing stratum is unavailable, instead, a preliminary stratigraphic correlation gives latest Bashkirian to middle Moscovian (Trümper et al. 2020). Therefore, I used the top age of Moscovian (306.9 Ma).

I used the MCMCTREE program in PAML 4.9j (Yang 2007) to estimate divergence times. Dating analyses used autocorrelated relaxed clock model and GTR+G model. Rate prior was set to 1 substitution per site per 100 Ma by reference to the empirical estimations (Papadopoulou et al. 2010, Andújar et al. 2012). Estimation of divergence times used the approximate method implemented in MCMCTREE. The first 20000 iterations were discarded as burn-in. 5000 samples were gathered, once every 200 iterations. A replicating MCMC was performed to check for convergence in TRACER.

Trees are visualized using FigTree 1.4.3 (Andrew Rambaut, http://tree.bio.ed.ac.uk/software/figtree) and modified using Adobe Illustrator CC 2017.

For the reader’s convenience and to enable a comparison of studies, familial taxonomy of cockroaches in this paper follows recent studies that are compared (e.g., Djernæs et al. 2015, Legendre et al. 2015, Wang et al. 2017, Bourguignon et al. 2018, Evangelista et al. 2019). It is noteworthy that recent molecular studies focusing on Blaberoidea raised subfamilies of Ectobiidae to families (Djernaes et al. 2020, Evangelista et al. 2020). Consequently, their Ectobiidae is identical to the Ectobiinae herein.

The relationship of major clades (suborder, superfamily, family, and subfamily) recovered herein is not identical to any previous studies. At the superfamily level, the sistergroup relationship of Corydioidea (only represented by Corydiinae) to the rest of Blattodea is consistent with that in Wang et al. (2017) and Djernæs et al. (2015, in part), both of which used three mitochondrial and at least two nuclear gene fragments, whereas conflicting with other recent phylogenies (Djernæs et al. 2015 in part, Legendre et al. 2015, Bourguignon et al. 2018, Evangelista et al. 2019). The superfamilial relationship of cockroaches is in dispute. On the other hand, the monophyletic Blaberoidea and the monophyletic blattoid complex (Blattoidea and Isoptera) are always supported.

Corydioidea are always undersampled. Species of Nocticolidae, Latindiinae, and Corydiidae incertae sedis (e.g. Ctenoneura) are lacking. Although one mitogenome of Nocticola is available, it is hardly serviceable unless the long branch is broken up by increased sampling (Poe 2003). Nonetheless, the transcriptome data support a monophyletic Corydioidea that include Corydiidae and Nocticolidae (Evangelista et al. 2019).

In the blattoid complex, only the sistergroup relationship between Cryptocercidae and Isoptera is universally recognized. These taxa constitute Xylophagodea (Engel 2011). The new phylogeny recovers Xylophagodea as sister to the remaining blattoid complex, of which the internal relationship is (Blattidae + Tryonicidae) + (Anaplectidae + Lamproblattidae). This is significantly different from other studies. Three major groups of Blattoidea are still undersampled, namely, Anaplectidae, Lamproblattidae and Tryonicidae. In addition, some mysterious taxa of Blattoidea incertae sedis (e.g. Oulopteryx) have not yet been sampled.

The paraphyly of Ectobiidae with respect to Blaberidae is a consensus among studies; the present study is not an exception. However, the relationships among Blaberidae and ectobiid subfamilies are conflicting among studies, especially in the positions of Ectobiinae and Pseudophyllodromiinae. The Ectobiinae contributes a weak point in the new phylogeny (pp = 79%), i.e. the node of Ectobiinae + (Blattellinae + Blaberidae). Regardless of the Nyctiborinae, which is not included in the final phylogeny herein, the sistergroup relationship between Blattellinae and Blaberidae is also supported in Bourguignon et al. (2018) and Evangelista et al. (2019, 2020). Although the considerably diversified Blaberidae are typically densely sampled, the phylogeny of them recovered by various studies is inconsistent.

A challenge to all molecular phylogenies is the reconciliation with morphological, behavioral, and other evidence. For example, oothecal property and rotation behavior are various and the taxonomic distribution of them in Blaberoidea is comparatively well known (McKittrick 1964, Roth 1967, 1968a, Bell et al. 2007). Note that the rotation behavior assigned to Ectobiinae in Evangelista et al. (2019) is contrary to the literature. The parsimonious scenario is that Pseudophyllodromiinae are sister to, or paraphyletic with respect to, the remaining Blaberoidea, which constitute a clade (e.g. McKittrick 1964, Klass and Meier 2006). The present study supports this scenario but the posterior probabilities of the “remaining Blaberoidea” node is relatively low (pp = 79%). Other studies, in which the phylogenies do not recover Pseudophyllodromiinae as sister to the remaining Blaberoidea, imply either parallel evolution of oothecal rotation in Ectobiinae and in Blattellinae + Nyctoborinae + Blaberidae, or an ancestral state of oothecal rotation in Blaberoidea and a loss in Pseudophyllodromiinae. At least, the relationship of Ectobiinae and Pseudophyllodromiinae to other Blaberoidea is debatable.

The only appropriate fossil calibration in Blattaria in the present study is Piniblattella yixianensis Gao et al., 2018 (Gao et al. 2018), which is used as a calibration for the first time. In the following, this fossil calibration is justified according to the five criteria in Parham et al. (2011).

Criteria 1 and 4. Information about the fossil-bearing stratum and museum collection is provided in Gao et al. (2018).

Criterion 2. Regardless of the determination of genus (which is in dispute, see Hinkelman 2019), P. yixianensis belongs in Blattellinae as evidenced by the oothecal rotation, reproduction type oviparity B, and the wing venation pattern, as explained below. Rotation feature and physical property of ootheca are crucial to the family-level phylogeny of cockroaches: the “advanced rotation” (i.e. rotating the ootheca and containing the anterior eggs inside vestibulum) is considered as a significant apomorphy in cockroaches, and distributed in Blaberoidea other than Pseudophyllodromiinae (McKittrick 1964, Roth 1967, 1968a, Bell et al. 2007). The rotation of the oothecae of P. yixianensis is unlikely due to taphonomic process: all preserved oothecae are horizontally positioned, none is perpendicularly or randomly positioned (Gao et al. 2018, Hinkelman 2019). The “primitive rotation” of Corydiidae is also ruled out: in the primitive rotation, the anterior eggs are outside the abdomen and the ootheca is obliquely positioned (Roth 1967) – this is not the case with P. yixianensis. Another difference between primitive and advanced rotation is the presence and absence of the flange, but which cannot be clearly observed in those fossils. In the phylogenies of the present study and some previous studies (e.g., McKittrick 1964, Klass and Meier 2006), Pseudophyllodromiinae are sister to or paraphyletic with the remaining Blaberoidea, supporting that the advanced rotation is autapomorphic for Blaberoidea excluding Pseudophyllodromiinae, and therefore P. yixianensis can at least calibrate the node of crown Blaberoidea (regardless of other evidence discussed below). However, other studies support several origins of the advanced rotation in Blaberoidea or loss of the advanced rotation in Pseudophyllodromiinae (e.g., Wang et al. 2017, Bourguignon et al. 2018, Evangelista et al. 2019). Under this hypothesis and regardless of other evidence (discussed below), P. yixianensis may only calibrate the split between Blaberoidea and the sister group with caution.

The reproduction type of P. yixianensis is oviparity B: (1) during reproduction, female cockroaches have a period of carrying the ootheca (if present) outside, but only the oviparity B carries the ootheca externally until hatching; other types only carry shortly before oviposition (oviparity A) or before retraction (ovoviviparity and viviparity) (Roth 1967, 2003, Bell et al. 2007), and have much less chance to leave fossils. Oviparity B likely contributes a lot to the preservation of ootheca-bearing fossils like P. yixianensis. (2) Based on the author’s observation during collecting, the oviparity A ootheca is easily dropped when the cockroach is caught, and almost certainly detached in the end. In comparison, some of the oviparity B oothecae remain attached in the abdomen even when the cockroach is preserved. This implies that it is unlikely that oviparity A cockroaches preserve fossils carrying oothecae, but oviparity B cockroaches may. (3) The oothecal keel of P. yixianensis is relatively underdeveloped and unornamented (Gao et al. 2018, Hinkelman 2019), and so accords with the characteristics of oviparity type B (Roth 1968a, 1971, Bell et al. 2007). Among living cockroaches, only some species of Blattella, Chorisia and Onycholobus are known of both the oviparity type B and the advanced rotation (McKittrick 1964, Roth 1967, 1968a, 1971, 1983, 2003, Bell et al. 2007). Chorisia and Onycholobus are not included in analyses here, but the former is considered as closest to Blattella (Roth 1983), while the latter has the ootheca resembling that of Blattella (Roth 1971).

However, oviparity B is homoplastic among Blaberoidea. Roth (1968b) found Lophoblatta, a pseudophyllodromiine genus, carrying the ootheca until the eggs hatch, but not rotating the ootheca. This discovery demonstrated that the oviparity B originated independently more than once within Blaberoidea. According to the reasonable hypothesis of Roth (1968a), oviparity B is the intermediate form between the ovoviviparity and the plesiomorphic oviparity A, i.e., ovoviviparity derived from oviparity B, which derived from oviparity A. Ovoviviparity occurs in most Blaberidae (with advanced rotation) but also, homoplastically, in two genera of Blattellinae, which rotate the ootheca (Roth 1982, 1984: Stayella, Roth 1995: Pseudoanaplectinia), and two genera of Pseudophyllodromiinae, which do not rotate the ootheca (Roth 1989: Sliferia, Roth 1997: Pseudobalta) (see also a review by Djernæs et al. 2020). Provided that Roth’s hypothesis is true, the homoplasy of ovoviviparity further demonstrates that the oviparity B is highly homoplastic.

Accordingly, it appears that Blaberoidea are preadapted to the advanced rotation and oviparity B (consequently ovoviviparity), but as far as known, these two features only co-occur in Blaberidae and Blattellinae. Blaberidae and Blattellinae were recovered as sister groups (Bourguignon et al. 2018 and the present study), or form a clade together with Nyctiborinae (Klass and Meier 2006, Evangelista et al. 2019), which is not included in analyses herein. This implies that species of the clade Blaberidae + Blattellinae (or Blaberidae + Blattellinae + Nyctiborinae) are more preadapted to allow (if not achieved) the co-occurrence of advanced rotation and oviparity B (consequently ovoviviparity). Therefore, the combination of advanced rotation and oviparity B may tentatively place P. yixianensis into that clade but outside of Blaberidae.

Other characters preserved in the fossils of P. yixianensis are barely discernible except the wing venation. The forewing of P. yixianensis conforms to the general form of Blattellinae (see Rehn 1951, Li et al. 2018): ScP with few branches, R pectinate proximally and dichotomous or irregular distally, M and CuA both developed and not essentially pectinate, claval furrow with sharp apical turn, and claval veins diagonal. These traits are distinct from other subfamilies of Ectobiidae. The hindwing has a simple ScP, a pectinate RA with four branches or so, a non-pectinate RP, a simple and feeble M, and a nearly pectinate CuA (Gao et al. 2018). This combination of hindwing traits is not characteristic of any taxon, although these traits are more common in Blaberoidea, particularly Blaberidae (see Rehn 1951, Li et al. 2018). Unfortunately, the polarity of wing venation characters above is barely clear, so that it is premature to conclude a phylogenetic position through these similarities in venation. It is noteworthy that both the tegmen and hindwing of P. yixianensis exhibit a developing characteristic posterior branch of R, i.e., the apicoposterior part of R is a branch with terminal branching only. Most cockroaches do not have a characteristic posterior branch (cpb), and this specialization is homoplastically derived among cockroaches, principally Ectobiidae (see Rehn 1951, Li et al. 2018). Nonetheless, the cpb and developing cpb vary in morphology, whereas the branching pattern of P. yixianensis is found in Blattella and related genera such as Episymploce, but not seen in others (Li et al. 2018, and unpublished observation). This evidence reinforces the hypothesis that P. yixianensis belongs in Blattellinae, although it is premature to conclude that P. yixianensis is a close relative of Blattella or even sister to Blattella.

So far, the evolution of reproduction type, ootheca handling behaviour and wing venation of cockroaches is not well understood, and might be more complicated than currently known. In view of this, the phylogenetic position of P. yixianensis is not securely settled. Nevertheless, P. yixianensis can be tentatively considered as a member of Blattellinae, and thus calibrates the node of Blattellinae + sister (Blaberidae herein). In summary, P. yixianensis as a calibration should be used with caution, and comparative analyses with/without this fossil should be performed to accommodate its uncertainty.

Criterion 3. Reconciliation between molecular and morphological phylogenies is partially achieved. As mentioned above, regardless of the Nyctiborinae that is not included in the final data, the sistergroup relationship between Blattellinae and Blaberidae is supported herein and in recent big data analyses (Bourguignon et al. 2018, Evangelista et al. 2019, 2020), and Ectobiinae and Pseudophyllodromiinae are not nested in the clade of Blattellinae + Blaberidae + Nyctiborinae. In the most comprehensive ever morphological (and ethological) phylogeny of Dictyoptera (Klass and Meier 2006), the above relationships within Blaberoidea are also recovered, except the absence of Ectobiinae. According to the phylogenetic discussion in McKittrick (1964), the Ectobiinae is nested in the clade of Blattellinae + Blaberidae + Nyctiborinae (Fig. 1), but that study is somewhat outdated and not strictly phylogenetic. There is a lack of recent morphological phylogeny covering all major groups; therefore, it is impossible to thoroughly compare the morphological phylogeny with the molecular phylogeny.

Criterion 5. Piniblattella yixianensis is from Huangbanjigou, Beipiao, Liaoning, northeastern China (Gao et al. 2018). Isotopic age of the fossiliferous layers in Huangbanjigou ranges from 121.2 Ma to 129.8 Ma (Swisher et al. 1999, Yang et al. 2007). However, the horizontal correlation between this fossil and the radiometric samples is unknown, and the radiometric sampling is insufficient, therefore the age range above does not necessarily represent the age of fossils. I conservatively use the age of the overlying stratum (of top Yixian Formation), 120.9 Ma (Smith et al. 1995), as the minimum age of P. yixianensis.

The other fossil for calibration, V. brenanae, has been frequently used for Xylophagodea (e.g., Misof et al. 2014, Bourguignon et al. 2018, Evangelista et al. 2019). Its identity as a termite is secured by the presence of basal suture (Jarzembowski 1981), one of the defining characters (autapomorphies) of Isoptera (Ax 1999, Krishna et al. 2013). Its validity as a calibration was justified by Wolfe et al. (2016), and I have no comments on this fossil.

Fossil calibrations contribute considerably to the discrepancy in the age estimation among studies. The ‘roachoid’ fossils, remarkably, were frequently assigned as “stem Dictyoptera” (e.g., Legendre et al. 2015, Tong et al. 2015, Bourguignon et al. 2018, Evangelista et al. 2019). Although the ambiguity of them and alternative interpretations were considered (e.g., Kjer et al. 2015, Bourguignon et al. 2018, Li and Huang 2019), a formal report on the impact of them is lacking. Because of the same assignment of ‘roachoid’ fossils, the age estimates with three calibrations herein (which is only for comparison) are close to that in Bourguignon et al. (2018) and Evangelista et al. (2019) (Permian origin of crown Dictyoptera and Triassic origin of crown Blattodea), and only somewhat younger than that in Legendre et al. (2015). ‘Roachoid’ fossils often play a decisive role in the dating of cockroaches, pushing the age estimates older. In comparison, the formal age estimates herein (two fossil calibrations) are close to that in Misof et al. (2014), both studies do not use ‘roachoid’ fossils as calibrations, suggesting Jurassic origins of crown Dictyoptera and crown Blattodea. Without the ‘roachoids’, other fossils will take over them as decisive calibrations and result in various, and usually younger, estimations (e.g., Wang et al. 2017, Bourguignon et al. 2018: fig. S12). A comparison of age estimates among studies is shown in Fig. 4.

Figure 4. 

Comparison among the ages estimated in various studies. The fossils are: Valditermes brenanae Jarzembowski, 1981; Piniblattella yixianensis Gao et al., 2018; Nodosigalea burmanica Li & Huang, 2018; Cretaperiplaneta kaonashi Qiu et al., 2020; Stegoblatta irmgardgroehni Anisyutkin & Gröhn, 2012. Abbreviation: ALTB, Anaplectidae + Lamproblattidae + Tryonicidae + Blattidae.

Unfortunately, many of the fossil calibrations other than ‘roachoids’ are also unjustified. Subsequently, comparisons among the age estimates from various studies could be pointless. For example, the “stem Mantodea” Homocladus grandis (Djernæs et al. 2015, Bourguignon et al. 2018) is highly questionable (Evangelista et al. 2019 and references therein); the “oldest Mantoidea fossil” Prochaeradodis enigmaticus (Djernæs et al. 2015, Wang et al. 2017) may be a cockroach (Cui et al. 2018); the “blattid” Balatronis libanensis is unlikely a true cockroach (Blattaria), not to mention Blattidae (Evangelista et al. 2017, Qiu L et al. 2020a); the “Diploptera fossils” (Bourguignon et al. 2018) cannot be identified to Diploptera and the higher-rank placement of those fossils also remains undetermined (Evangelista et al. 2017, Li et al. 2017); the “first modern cockroach” Zhujiblatta anofissilis (Bourguignon et al. 2018) is phylogenetically unsettled and has to be redescribed, which I am preparing elsewhere. Even more surprisingly, an unnamed “Epilampra fossil” found in an extant cockroach database was used (Bourguignon et al. 2018).

A critical review of cockroach fossil calibrations was not achieved until Evangelista et al. (2017), who recommended four cockroach fossils for node calibration. Evangelista et al. (2019) discarded one of them and retained corydiid Cretaholocompsa montsecana Martinez-Delclos, 1993, blaberid “Gyna” obesa (Piton, 1940) and ectobiid Ectobius kohlsi Vršanský et al., 2014. However, these fossils are still debatable.

Cretaholocompsa montsecana was determined as a close relative of extant Holocompsa (Martinez-Delclos 1993), and used as a calibration for corydiid nodes (Legendre et al. 2015, Wang et al. 2017, Evangelista et al. 2019). However, Cretaholocompsa significantly differs from Holocompsa in the presence of large spines along the ventral margin of midfemora (other legs unknown) (Qiu L et al. 2020b). Besides, according to recent accounts of Corydiidae, such spines are absent from all femora in this family (e.g., Estrada-Alvarez and Guadarrama 2012, Hopkins 2014, Crespo et al. 2015, Qiu L 2017, Qiu L et al. 2017, 2019a, 2019b, 2020b).

“Gyna” obesa was used as a calibration for blaberid nodes (Bourguignon et al. 2017, Evangelista et al. 2019). Evangelista et al. (2017) identified this fossil to Blaberidae based on (1) the stout cerci, (2) approximately parallel edges of tegmina, (3) elongated CuP, (4) shape of the pronotum, (5) asymmetrical male subgenital plate and (6) large body size. Although this fossil appears to be blaberid in overall appearance, the evidence above is weak or invalid to place “Gyna” obesa in Blaberidae. First, as acknowledged by Evangelista et al. (2017), traits 1, 2 and 3 are homoplastic with Blattidae. Second, the shape of the pronotum is not entirely clear (Evangelista et al. 2017). Third, the subgenital plate is not clearly discernible: according to the figures in Evangelista et al. (2017), the “concave margin” is implausible in favour of poor preservation. Besides, the subgenital plate appears to be large and cover three segments as in extant female cockroaches, in comparison to the male subgenital plate that is as small as one normal segment; therefore, the specimen may be a female. Fourth, large-sized species are common in Blattidae and Corydiidae in addition to Blaberidae.

Ectobius kohlsi was identified based on a comparison with extant species (Vršanský et al. 2014). However, the preserved characters are not unique enough to indicate the genus; i.e., diagnostic characters of Ectobius or of Ectobiinae are not clearly observed, e.g., elongate male genital elements (Roth 2003) and a pinnate R+M+CuA system of tegmina (unpublished observation). Instead, spot-and-line macula patterns on the pronotum and forewings are common in Ectobiinae and Pseudophyllodromiinae. Vršanský et al. (2014) reported a female with valvate subgenital plate. If this is true, then this species likely belongs to Pseudophyllodromiinae because suchlike females do exist in Pseudophyllodromiinae (e.g. Euphyllodromia, Anisyutkin 2011). If this fossil is a male, then it is more likely to be a member of Pseudophyllodromiinae, some genera of which have a valvate subgenital plate in males (e.g. Balta, Qiu ZW et al. 2017).

Interestingly, the fossil discarded by Evangelista et al. (2019), Cariblattoides labandeirai Vršanský et al., 2011, is likely to be a genuine pseudophyllodromiine species, even though the genus is uncertain. According to Vršanský et al. (2011b), the forewing of C. labandeirai bears venational characters in common with most Pseudophyllodromiinae: ScP simple and short, R essentially pectinate, M pectinate and more developed than CuA, claval veins oblique or diagonal (see Rehn 1951, Li et al. 2018). The venation alone is a weak indicator of the taxonomic identity, whereas the combination of venation and macula pattern is stronger reasoning. Unfortunately, C. labandeirai is not suitable for calibrating Pseudophyllodromiinae + sister (= Blaberoidea herein) even if its inclusion in Pseudophyllodromiinae is proven, because it would be suppressed by the older fossil, P. yixianensis, which already calibrates an internal node.

Only one true-cockroach fossil is used as a calibration in the present study, but this does not imply that other fossils are substandard. Every informative fossil (with high phylogenetic resolution and ascertained geological context) has the potential to be a competent calibration, but the incorporation of them is hampered by the fact that relevant living species are under-sampled or have not yet been sequenced. Noteworthy examples of fossils include those of extant genus, e.g. Supella (Nemosupella) miocenica Vršanský et al., 2011 (see Vršanský et al. 2011a), and those of Corydioidea, e.g., Proholocompsa fossilis (Shelford, 1910) (see Gorokhov 2007), Paraeuthyrrhapha groehni Anisyutkin, 2008 (see Anisyutkin 2008), Crenocticola Li & Huang, 2019 (see Li and Huang 2019). These fossils could become powerful calibrations for smaller clades (younger nodes) if the data of related extant species are available, otherwise they can only calibrate larger clades (older nodes), and become ineffective when older fossils (e.g. those used herein) calibrate the same or internal nodes.