The study includes 131 ingroup taxa (Blattoidea), 19 near outgroup taxa (other Dictyoptera), and 9 far outgroup taxa (other Polyneoptera + Odonata), for a total of 159 taxa. The ingroup includes multiple representatives for all major blattoidean lineages (Blattidae, Lamproblattidae, Tryonicidae, Anaplectidae, and Cryptocercidae + Isoptera). The ingroup also includes multiple representatives for all currently recognised blattid subfamilies (Archiblattinae, Blattinae, Macrocercinae, and Polyzosteriinae) and a representative of the genus Duchailluia (Blattinae). Several putatively conspecific, but genetically distinct, individuals are included in the data set (see Table S1). This is the case for e.g. Lamproblatta albipalpus Hebard, for which two specimens, L. albipalpus MD-2014 and L. albipalpus TB-2018 exhibit at least 8% sequence divergence in all overlapping sequences; percent sequence divergence based on BLAST results from GenBank (www.ncbi.nlm.nih.gov/genbank). In general, sequences from different specimens were not combined as a single terminal taxon if the sequence divergence exceeded 3% in any overlapping sequences. The only exception to this rule was Eurycotis floridana (Walker) in which the complete mitochondrial sequence (GenBank # MG882177) was very similar to sequences from other included conspecifics based on 12S and COI+II (< 1% divergence), but highly divergent based on 16S (> 7% divergence). The 16S sequence used for E. floridana (GenBank # KP986295) is very similar (< 1% divergence) to the third available E. floridana 16S sequence (GenBank # JN615296). Named species were combined with congeneric ‘sp.’s or ‘cf.’s if the sequence divergence of all overlapping sequences was less than 1%. These cases are indicated in Table S1 as e.g. Shelfordella lateralis (Walker)/sp. MNHN BL113, but otherwise just referred to by the species name, here S. lateralis. We generally follow the taxonomy of Cockroach Species File (Beccaloni 2014), but within Blaberoidea we follow Djernæs et al. (2020) and Evangelista et al. (2020).

The data set consists of sequences from 9 genes that are widely used for resolving cockroach relationships (e.g. Inward et al. 2007; Djernæs et al. 2012, 2015, 2020; Legendre et al. 2015, 2017; Wang et al. 2017). Six mitochondrial genes: the ribosomal genes 12S (c. 375 nucleotides [nt] fragment) and 16S (c. 415 nt), the protein coding genes COI (c. 1490 nt fragment, but many sequences only a 658 nt fragment) and COII (c. 690 nt), and the tRNA genes tRNA-leu (c. 75 nt) and tRNA-lys (c. 40 nt fragment). Three nuclear genes: the ribosomal genes 18S (c. 1825 nt) and 28S (c. 520 nt fragment), and the protein coding gene H3 (c. 330 nt fragment). The total length of the aligned data set is 7147 nt. Our data set provides sequence data for 37 species of Blattidae and a species of Tryonicidae not previously sampled, as well as adding newly sequenced gene coverage to previously sampled species. Sequencing was done using standard methods: Samples contributed by JM (see Table S1) were extracted using the DNeasy Tissue Kit (Qiagen, Valencia, CA) and amplified using Ready-ToGo polymerase chain reaction Beads (GE Healthcare, Piscataway, NJ). The mitochondrial 12S and 16S fragments were amplified using the 12Sai ⁄12Sbi and 16SA⁄16SB primers respectively (Kambhampati 1995). The barcode fragment of the mitochondrial COI was amplified using the LCO1490/HCO2198 primers (Folmer et al. 1994). The nuclear 18S was amplified using the primer pairs 18S1F/18S5R, 18S3F/18Sbi and 18SA2.0/18S9R (Giribet et al. 1996; Whiting et al. 1997). The nuclear 28S was amplified using the 28Sa/28Sbout primers (Nunn et al. 1996). The nuclear H3 was amplified using the H3aF⁄H3aR primers (Svenson and Whiting 2004). The PCR protocol was 94°C for 5 min, followed by 35 cycles of 94°C (15 s), 50°C (15 s) and 72°C (15 s), and then a final extension of 72°C (7 min) on MJ Research Peltier Thermal Cyclers (MJ Research Inc., Waltham, MA). Each 25 µL reaction contained 1 µL of each 10 mm primer, 2 µL of template and 21 µL of water. PCR products were purified using AMPure magnetic bead purification (Agencourt Bioscience Corp., Beverly, MA) on a Biomek NX robot (Beckman Coulter, Fullerton, CA). Amplification products were then sequenced in both directions. Each reaction mixture contained 1 µL BigDye (Applied Biosystems), 1 µL of 3.2 mm primer, 1 µL BigDye Extender Buffer (Applied Biosystems) and 5 µL of DNA template. Sequencing reactions ran for 25 cycles of 96°C (15 s), 50°C (15 s) and 60°C (4 min). Sequences were purified using CleanSeq magnetic bead purification (Agencourt Bioscience Corporation) on a Biomek NX robot (Beckman Coulter) to remove unincorporated primers and dyes. Products were re-suspended in 40 µL of 0.5 mM EDTA and were electrophoresed in an ABI Prism 3730xl sequencer (Applied Biosystems). Samples contributed by MD (see Table S1) were produced following Djernæs et al. (2015). All new sequences were checked for contamination using unrestricted BLAST searches on GenBank. The new sequences were combined with sequences from GenBank and BOLD (http://www.boldsystems.org) to produce the data set.

The sequences were aligned in MAFFT 7.471 (Katoh et al. 2005, 2019; Katoh and Standley 2013) using the G-INS-1 algorithm. Distance trees (Neighbour Joining) were produced for each alignment in MAFFT to check for incorrectly identified GenBank or BOLD sequences. All alignments were checked visually and manual corrections were made in Mesquite v. 3.03 (build 702; Maddison and Maddison 2015).

Our taxon sampling approach resulted in a data set with missing data that could lead to lack of resolution and/or low support values. To alleviate this potential issue, we ran analyses both on the complete data set and on a reduced (trimmed) data set. To produce the trimmed data set, we excluded taxa that did not have coverage for at least three of the genes 12S, 16S, COI, COII, 18S, 28S or H3 (tRNA-leu and tRNA-lys not considered due to their short length). This resulted in the exclusion of 22 taxa from the trimmed data set: Catara rugosicollis (Brunner von Wattenwyl), Cartoblatta scorteccii Princis 471A, C. scorteccii 474A, Celatoblatta vulgaris (Johns), Hebardina concinna (Haan), Maoriblatta novaseelandia (Brunner von Wattenwyl), Pseudoderopeltis bimaculata (Walker), Macrocerca sp. 1 FL-2016, Anamesia maculosa Mackerras, Anamesia lambii Tepper, Desmozosteria scripta Mackerras, Desmozosteria cincta Shelford, Eppertia furcate (Tepper), Eppertia sp. ANIC, Eurycotis bahamensis Rehn, Euzosteria nobilis (Brunner von Wattenwyl), Euzosteria sordida Shaw, Megazosteria patula (Walker) 000186, Pallidionicus pandanorum Grandcolas, Polyzosteria limbata Burmeister, Zonioploca pallida Shelford, and Zonioploca sp. This exclusion did not affect our taxonomic coverage at the family or subfamily level, only at the genus and species level.

We partitioned the data according to origin (mitochondrial vs. nuclear) and gene type (protein coding vs. rRNA and tRNA), resulting in four partitions: 1) mitochondrial rRNAs and tRNAs (12S, 16S, tRNA-leu and tRNA-lys), 2) mitochondrial protein coding genes (COI and COII), 3) nuclear rRNAs (18S and 28S), and 4) nuclear protein coding gene (H3). We then added further partitioning by gene type (rRNA vs tRNA), and partitioning by level of variation (variable vs conserved regions of 18S and 28S; 12S and 16S in the present data set did not have the long conserved regions characteristic of 18S and 28S). Variable versus conserved regions was determined by visual evaluation of the alignment in Mesquite. This further partitioning resulted in six partitions: 1) mitochondrial rRNAs (12S and 16S), 2) mitochondrial protein coding genes (COI and COII), 3) mitochondrial tRNAs (tRNA-leu and tRNA-lys), 4) nuclear rRNAs, variable regions (variable regions of 18S and 28S), 5) nuclear rRNAs, conserved regions (conserved regions of 18S and 28S), and 6) nuclear protein coding gene (H3). We then added partitioning by gene, resulting in ten partitions: 1) 12S, 2) 16S, 3) COI, 4) COI, 5) tRNA-leu and tRNA-lys, 6) variable regions of 18S, 7) conserved regions of 18S, 8) variable regions of 28S, 9) conserved regions of 28S, and 10) H3. We tested additional partitioning by codon position (resulting in 16 partitions), but this led to a lack of convergence in the Bayesian Inference analyses and clear artefacts in the resulting trees (e.g. taxa included based on just the COI barcode fragment forming a clade apart from closely related taxa with more complete data), problems likely caused by over-parameterization (Ronquist et al. 2020).

We analysed both the complete and trimmed data sets using the three above-mentioned partitioning schemes (4, 6 and 10 partitions) allowing us to explore the effect of different partitioning schemes which can affect tree topology (Kainer and Lanfear 2015). Bayesian Inference (BI) analyses were performed in MrBayes 3.2.7a (Ronquist et al. 2012) on Cipres XCEDE (Miller et al. 2010) using model jumping with a gamma model for variation across sites. Model jumping allows MrBayes to sample across the entire GTR model space, a more elegant alternative to a priori model testing (Ronquist et al. 2020). The Temp parameter was set to 0.05, changed from the default of 0.1 to facilitate swapping between chains. Analyses of the trimmed data set were run for 10 million generations, sampled every 1000 generations, with burninfrac set to 0.5 and burnin set to 5000 trees. Analyses of the complete data set were run for 50 million generations, sampled every 10000 generation, with burninfrac set to 0.5 and burnin set to 2500 trees. Convergence was evaluated by a visual inspection of the log likelihood plot and by the average standard deviation of split frequencies being < 0.01. Maximum Likelihood (ML) analyses were performed in Garli 2.01 (Zwickl 2006) on Cipres XCEDE. For each partition we used the GTR submodel with the highest posterior probability from the corresponding BI analysis. Models and the overall rate were unlinked across partitions. Enallagma cyathigerum (Charpentier) was specified as outgroup to facilitate production of consensus trees. The ML analyses were terminated after 20,000 generations without significant change of topology. Ten independent runs of each ML analysis were performed and the best trees were chosen based on likelihood values (T-ML-4: 3 trees; T-ML-6: 9 trees; T-ML-10: 9 trees; C-ML-4: 3 trees; C-ML-6: 8 trees; C-ML-10: 2 trees). From these, a majority rule consensus tree was produced in Mesquite. Bootstrap analyses were run with 1000 bootstrap repetitions and one tree search repetition per bootstrap repetition, settings otherwise as above. Bootstrap values were calculated in Mesquite.

We chose to use Garli rather than e.g. IQ-TREE or RAxML as the latter programs were developed to analyse phylogenomic data with thousands of taxa and hundreds of thousands to millions of nucleotides (Stamakis 2008; Minh et al. 2021). Our complete data set contains 159 taxa and 7147 nucleotides, and is thus not the type or size of data set that IQ-TREE and RaXML were developed for.

Data on geographical distribution are mainly from Princis’ catalogue (Princis 1965a, 1966, 1971). Distribution data for species described after 1970 are from the original descriptions. Data from Princis (1965a, 1966, 1971) and original descriptions were supplemented with field collecting data when available. Data for various unnamed species (sp.) are from GenBank, from Bourguignon et al. (2018: table S1), and collecting data (new specimens). When several conspecifics are included as terminal taxa (e.g. three nominal Tryonicus parvus (Tepper); see criteria for inclusion in section 2.1.), the distribution data are based on specimen data, not data for the nominal species. The exception is Cartoblatta scorteccii, in which both included specimens were intercepted at port facilities; in this case the data from Princis (1966) were used. We scored distribution of widespread pest species as inapplicable.

We followed the definitions of biogeographic realms of Olson et al. (2001), with the modifications used by Djernæs et al. (2020). Thus, we included all of Mexico in the Neotropics, included Oceania in Australasia, divided the Palearctic into East and West (along the Ural Mountains) and included all of China (unless more specific locality information was available) in East Palearctic. Additionally, as more than half of the included Blattoidea are from Australasia, we used a finer scale of geographic distribution within this region. This resulted in 10 biogeographic areas: 1) Nearctic, 2) Neotropics, 3) Afrotropics, 4) East Palearctic, 5) Indo-Malaya, 6) Australia, 7) New Guinea, 8) New Caledonia, 9) New Zealand, and 10) remaining Australasia, the latter five subregions of Australasia (no Blattoidea from West Palearctic included).

We mapped geographic distribution in Mesquite on our preferred tree, the tree resulting from the Maximum Likelihood analyses of the complete data set with 6 partitions (C-ML-6). We used parsimony reconstruction of ancestral states and treated the characters as unordered. When a terminal taxon occurred in more than one biogeographical area, we scored it as present in all relevant areas.

The phylogenetic analyses generally gave consistent results with all analyses finding the same relationships between the major lineages of Blattoidea (Blattidae, Lamproblattidae, Tryonicidae, Anaplectidae, and Cryptocercidae + Isoptera). Lamproblattidae was sister to the remaining Blattoidea and Tryonicidae was sister to Cryptocercidae + Isoptera (Fig. 1). All the major blattoidean lineages were monophyletic, except Anaplectidae: Anaplecta sp. FL-2015 was consistently placed outside Blattoidea as sister to Pseudomops oblongatus (Linnaeus) (Blattellidae). See Table S2 and Figures S1–S12 for support values. Within Blattidae, none of the currently accepted subfamilies were monophyletic and some relationships between clades differed between analyses, although most analyses found the same relationships between clades as depicted in Fig. 1. Differences from this pattern were found in some analyses of the complete data set, see Table S2. We chose the tree from the Maximum Likelihood analysis of the complete data set (C-ML-6) as our ‘preferred tree’. This tree has the same relationships between clades as found in the majority of analyses, and it was based on the complete data set, thus showing the phylogenetic position of all the included taxa.

Figure 1. 

Tree from Maximum Likelihood analysis of the complete data set using 6 partitions (C-ML-6). Subfamily, tribe and genus names and assignments reflect the taxonomic changes made in the present paper; ‘old’ genus names are given in parentheses (in grey). Bootstrap support values for clades of interest are shown, # indicate a bootstrap support < 50. Trees from other analyses essentially agree with this tree, with some minor differences within Blattidae (see Table S2).

The trees based on the trimmed data set generally had higher support values, both posterior probabilities and bootstrap values (Table S2). In the Bayesian analyses, the trees based on the trimmed data set were more resolved, while all Maximum Likelihood trees (majority rule consensus trees) were completely resolved. Bayesian analyses of the trimmed data set reached convergence (average standard deviation of split frequencies < 0.01) much faster than analyses of the complete data set (trimmed data set 1.8–8.7 million generations, 6–25 h; complete data set 30.3–42.3 million generations, 110–125 h).

Distribution mapping reveals clear geographic structuring with many clades restricted to one or two (neighbouring) geographic areas (Fig. 2). In general, the geographic structuring is clearer in the more densely sampled parts of the tree (compare e.g. Isoptera and Blattidae). Within Blattidae, the majority of Blattinae and Archiblattinae occurs in the Afrotropics and Indo-Malaya, while the majority of Polyzosteriinae is restricted to Australasia. Within Polyzosteriinae, several geographically restricted clades occur (restricted to Australia + New Guinea, to New Caledonia and to New Zealand). It is especially worth noting the clear geographic structuring in Blattidae compared to the lack of monophyly of the currently accepted subfamilies.

Figure 2. 

Tree showing geographical distribution of Blattoidea and classification of Blattidae into subfamilies. Within Blattidae, geographic distribution shows greater congruence with phylogeny than the ‘old’ classification does. See section 2.4. for definition of regions and Table S1 for more detailed information on distribution. The tree is from analysis C-ML-6 (Fig. 1). Both the ‘old’ classification of Blattidae and the revised classification introduced in the present paper are shown.

The relationships between the major lineages of Blattoidea were consistent between all analyses, generally with a pp > 90, but with bs < 50 in all cases (see Table S2). We found Lamproblattidae as sister to the remaining Blattoidea (pp 100, bs 66–82), and Tryonicidae as sister to Cryptocercidae + Isoptera (pp 75–95, bs < 50). Anaplectidae was sister to Tryonicidae + Cryptocercidae + Isoptera (pp 94–98, bs < 50), and this clade was in turn sister to Blattidae (pp 92–97, bs < 50). This agrees partially with Djernæs et al. (2015: fig. 3), but not with Wang et al. (2017), Bourguignon et al. (2018), Evangelista et al. (2018, 2019), Bläser et al. (2020) or Djernæs et al. (2020). Wang et al. (2017) found Anaplectidae as sister to Cryptocercidae + Isoptera and placed Tryonicidae and Lamproblattidae as sister groups. Bourguignon et al. (2018) found Blattidae + Tryonicidae as sister to Cryptocercidae + Isoptera and placed Anaplectidae and Lamproblattidae as sister groups. Evangelista et al. (2018, 2019) and Bläser et al. (2020) found Lamproblattidae as sister to Cryptocercidae + Isoptera and placed Blattidae and Tryonicidae as sister groups. Thus, no consensus exists regarding the relationships between the major blattoidean lineages. Furthermore, both Djernæs et al. (2015) and Evangelista et al. (2018) reported lack of consistency between analyses; the positions of Lamproblattidae and Anaplectidae respectively were very labile. Wang et al. (2017: fig. 1) also found that the relationships between the major lineages differed between analyses. Compared to these previous studies, our study presents an improved taxon sampling, especially for Lamproblattidae and Tryonicidae, which may be responsible for the increased consistency between analyses run under different conditions.

Morphological phylogenetic analyses of these taxa did not recover Blattoidea (Klass and Meier 2006; Djernæs et al. 2015: fig. S13), and so are of limited use regarding the reconstruction of relationships between the major lineages.

Thus, neither molecular nor morphological data offer any firm conclusions regarding these relationships. However, if the internal consistency between analyses in this study is due to increased taxon sampling, increasing the taxon sampling could lead to better consistency between studies. In Blattidae, multiple representatives from all subfamilies should be included if possible (see section 4.6. for revised classification of Blattidae). Additionally, an increased taxon sampling of Anaplectidae, Lamproblattidae and Tryonicidae is desirable, see sections 4.2., 4.3., and 4.4. for details. The taxon sampling of Cryptocercidae + Isoptera is good across most studies, leaving limited room for improvement, although gene coverage for Cryptocercus clevelandi Byers could be improved.

Anaplectidae (excl. Anaplecta sp. FL-2015, see above) was divided into two groups, a Neotropical group and a group found in East Palearctic (China) and Australasia (Australia) with high support (pp 99–100, bs 99–100). Within the latter group, the two Australian specimens form a clade; both are identified as Anaplecta calosoma Shelford, but show > 17% divergence in the mitochondrial DNA, so are unlikely to belong to the same species. For comparison, a recent study of ten Chinese Anaplecta species found the largest interspecific difference in COI to be 16.8% (Deng et al. 2020). Furthermore, the type locality of A. calosoma is New Guinea (Shelford 1912), so it is possible that neither of the Australian specimens identified as A. calosoma actually belong the species described by Shelford.

As currently defined, the family Anaplectidae contains two genera [Anaplecta (102 species) and Maraca (one species)] and occurs across five biogeographic regions [Neotropic, Afrotropic, East Palearctic, Indo-Malaya and Australasia (Beccaloni 2014; Deng et al. 2020)]. Anaplecta species from the Neotropics, East Palearctic, Indo-Malaya and Australasia have been included as terminal taxa in this and other phylogenetic studies (several Chinese species, including A. omei Bey-Bienko, occurs both in the East Palearctic and in Indo-Malaya, Deng et al. 2020). African Anaplecta have never been included in modern phylogenetic studies (unless Anaplecta sp. FL-2015 is a representative of the African Anaplectidae; no locality information was available from GenBank or the original manuscript). McKittrick (1964) studied three species of Anaplecta, two Neotropical and one African, and noted that they form a very distinct group, despite differences in the female genitalia indicating that they should not all belong to the same genus. Thus, the genus Anaplecta as currently defined is likely a coherent unit, but might need to be split into two or more genera as it consists of a number of morphologically and genetically distinct groups (McKittrick 1964; Djernæs 2015; Wang et al. 2017; this paper).

Future phylogenetic studies of Anaplectidae should include the Neotropical Maraca fossata Hebard to confirm whether it belongs in Anaplectidae, and if so, if it should be a separate genus or belongs among the Neotropical members of Anaplecta. The latter might be quite likely as Hebard (1926) considered M. fossata particularly close to specific Neotropical Anaplecta species, although the antenna and pronotum are different.

Apart from the question of the correct placement of M. fossata, some subdivision of the many species placed in Anaplecta would be appropriate as also suggested by McKittrick (1964). The geographical divisions within Anaplecta found in this study are consistent with results of other studies (Djernæs et al. 2015, 2020; Wang et al. 2017; Evangelista et al. 2018). If a denser sampling of Anaplecta with better geographic coverage confirms well-supported geographical groups, this might form the basis for a division of Anaplecta into several genera. For this, an increased sampling of Australian and Asian taxa is desirable, but inclusion of African Anaplectidae could be even more important as these have never been included in phylogenetic studies.

As currently defined, the family Lamproblattidae contains three Neotropical genera, Lamproblatta, Lamproglandifera and Eurycanthablatta (Beccaloni 2014). So far, only Lamproblatta species have been included in phylogenetic studies, the genus often represented by a single specimen (e.g. Djernæs et al. 2012, 2015, 2020; Legendre et al. 2015; Wang et al. 2017; Evangelista et al. 2018, 2019), although Bourguignon et al. (2018) included two specimens. We included four representatives of Lamproblatta, and found that the deepest split was between a Costa Rican Lamproblatta, nominally L. albipalpus, and the three Lamproblatta specimens from French Guiana, one of which is also nominally L. albipalpus. The two specimens identified as L. albipalpus show > 8% divergence in the mitochondrial DNA (12S 9%, 16S 11%, COI+II 12%), so are unlikely to belong to the same species. Furthermore, both specimens originate in localities rather distant from the type locality in Colombia (Hebard 1919); Lamproblatta is flightless, so is not expected to be highly mobile. Thus, neither of the two nominal L. albipalpus specimens might actually be L. albipalpus, underlining the need for a revision of this genus, which was last done by Rehn (1930).

The inclusion of Lamproglandifera and Eurycanthablatta in future studies is desirable both to test their placement in Lamproblattidae, and to improve the taxon sampling of this family. The inclusion of Eurycanthablatta is especially important as Fritzsche et al. (2008) noted similarities in the male genitalia to both Lamproblatta and Cryptocercus; thus, Eurycanthablatta might be a key taxon when trying to find the sister group of Cryptocercus + Isoptera.

The family Tryonicidae as currently defined contains two genera, Tryonicus and Lauraesilpha, occurring in Australia and New Caledonia (Beccaloni 2014), with the two genera considered sister taxa (Murienne et al. 2008; Legendre and Grandcolas 2018), a placement supported by Murienne et al.’s (2008) and Murienne’s (2009) phylogenetic analyses. However, results from Djernæs et al. (2015, 2020) and Wang et al. (2017) suggests that Tryonicus is paraphyletic with respect to Lauraesilpha. The increased sampling of Tryonicus in this study (seven terminal taxa compared to at most two in any previous studies) allowed us to investigate the possible paraphyly of Tryonicus more thoroughly. Our results place Lauraesilpha deep within Tryonicus, consistent with Djernæs et al. (2015, 2020) and Wang et al. (2017), but as sister to the New Caledonian species of Tryonicus, consistent with Murienne et al. (2008) and Murienne (2009).

Our sampling of Tryonicus includes T. mackerrasae Roth (Australia), three specimens of T. parvus (Australia), T. vicina (Chopard) (New Caledonia), and two New Caledonian Tryonicus sp. Within T. parvus, specimen FL-2015 shows 9–13% difference in mitochondrial DNA to MD-2014 and TB-2018, although the three specimens do form a highly supported clade (pp 100, bs 100).

These results highlight the need for a revision of Tryonicidae. Roth (1987) revised Tryonicus (including Lauraesilpha angusta (Chopard), which was then Tryonicus angusta) and considered Tryonicus monteithi Roth, T. mackerrasae and a possible third species more closely related to each other than to T. parvus or L. angusta. Our results suggest that either Lauraesilpha should be subsumed in Tryonicus, or Tryonicus should be split into at least three genera: one containing T. mackerasae (and likely T. monteithi, Roth 1987), one containing T. parvus sensu lato, and one containing the New Caledonian Tryonicus species. As Lauraesilpha forms a well-supported monophyletic group (Murienne et al. 2008; Murienne 2009; this study) with a distinctive mode of life (xylophagous, tunnelling in rotten branches) compared to other Tryonicidae which live under stones or pieces of dead wood (Roth 1987; Grandcolas 1997), we favour splitting Tryonicus rather than subsuming Lauraesilpha. However, as there seems to be several undescribed species of Tryonicus both in Australia (Roth 1987; this study) and in New Caledonia (Murienne pers. obs.), any taxonomic changes should await a better taxon sampling and description of new species. With regard to the exact placement of Tryonicidae within Blattoidea, an increased sampling of Australian Tryonicidae might be beneficial, due to the subordinate position of the New Caledonian taxa.

Cryptocercidae + Isoptera formed a highly supported monophyletic group in all analyses (pp 100, bs 94–100), consistent with previous studies (e.g. Inward et al. 2007; Djernæs et al. 2015; Legendre et al. 2015; Wang et al. 2017; Bourguignon et al. 2018; Evangelista et al. 2018, 2019). The relationship between the included Cryptocercus species is consistent with that found by Che et al. (2016, 2020) and the relationship between the included Isoptera species is consistent with the results of e.g. Cameron et al. (2012), Bourguignon et al. (2014), Wu et al. (2018), and Bucek et al. (2019). However, the taxon sampling in our study does not provide any new information about the relationships within Cryptocercidae or Isoptera compared to these studies, which focused explicitly on these relationships (for Cryptocercidae, see also e.g. Lo et al. 2006; Bai et al. 2018).

Blattidae was monophyletic in all analyses, but none of the subfamilies as currently defined (Archiblattinae, Blattinae, Macrocercinae and Polyzosteriinae) were monophyletic. Additionally, several genera could not be placed in any of these subfamilies, emphasizing the need for a revised classification of Blattidae. The non-monophyly of the blattid subfamilies is generally consistent with other studies with a sizeable sampling of blattids (Macrocercinae have not previously been included in any phylogenetic studies). Legendre et al. (2015) and Evangelista et al. (2018) found both Blattinae and Polyzosteriinae to be non-monophyletic (only a single representative of Archiblattinae included in both studies). Bourguignon et al. (2018) found Blattinae, but not Polyzosteriinae to be non-monophyletic (only a single representative of Archiblattinae included) while Wang et al. (2017) and Liao et al. (2021) found Blattinae, Polyzosteriinae, and Archiblattinae to be non-monophyletic. Distribution mapping showed a clear geographic structuring within Blattidae (Fig. 2), and a greater congruence between phylogeny and distribution than between phylogeny and current classification, similar to the patterns found by e.g. Svenson and Whiting (2009, Mantodea), Simon et al. (2019, Phasmatodea) and Djernæs et al. (2020, within Blaberidae). Based on the results from this and other studies, we present and discuss a revised classification of Blattidae. The revised placement of all blattid genera included or discussed in this study is shown in Table 1 with the relevant references.