Mathematica 10

Demultiplexed FASTQ files had adapters removed and low-quality bases trimmed (Krueger, 2017; Trimgalore v. 0.4.3 options: -q 20; -stringency 3; -e 0.1; -length 30). Reads were assembled with Trinity v.2.0.6 (Grabherr et al., 2011; Haas et al., 2013) . Orthologs from OrthoDB v. 7 (Waterhouse et al., 2013) were identified in the 90 sequenced libraries and 37 additional Blattodea and Polyneoptera sequences (Evangelista et al., 2019) using Orthograph v. 0.6 (Petersen et al., 2017) with Drosophila melanogaster, Pediculus humanus, Tribolium castaneum, and Zootermopsis nevadensis as reference genomes [default options except -anysymbol option called in MAFFT (Katoh and Standley, 2013) , and blast-max-hits = 50]. 265 targeted loci were extracted from the ortholog set. Quality filtered reads are available on NCBI GenBank (SRP155429).
Each locus was aligned in MAFFT v. 7.3 (Katoh and Standley, 2013; options: -retree 4maxiterate 10 -adjustdirection) and then trimmed from the edges to eliminate leading or trailing sites missing >80% of data. A second, alignment was conducted in MAFFT v. 7.3 (-localpair -maxiterate 1000), which was then adjusted to maintain consistent reading-frame. Alignments were finalized with manual adjustment in AliView v. 1.18 (Larsson, 2014) to remove poor quality reads and correct misaligned sections. Custom scripts in Mathematica 10 (Wolfram Research, 2012) available in the package Phyloinformatica v. 0.9 (Evangelista, 2019) were used to manage sequence files, translate sequences, trim and concatenate alignments. We refer to the final concatenated alignment (127 taxa) as the "265_Full" alignment.
We ran PartitionFinder 2 (Lanfear et al., 2016) with the 265_Full alignment with blocks defined as codon positions per locus, possible models as GTR and GTR+G, branch lengths considered as unlinked, best model chosen with AICc, and rcluster search scheme (percent = 10; max = 1000). Using the resulting codon partitioning scheme, we ran a preliminary tree search with 10 independent runs in RAxML v. 8.2 (Stamatakis, 2014) , implemented on the CIPRES portal (Miller et al., 2010) . Assessment of the best preliminary tree showed that a few taxa (Amazonina platystylata, Doradoblatta sp., Ischnoptera galibi, Lanxoblatta sp., Panchlora stolata, Pycnoscelus femapterus, Pycnoscelus striata) had exceptionally long branches. The same taxa were among those with the largest proportion of missing data (supplementary data). After reassessing the alignments in which these species were present, we removed Pycnoscelus striata, P. femapterous, Ischnoptera galibi, Amazonina platystylata, and Doradoblatta sp. from the analysis under the grounds that (i) their data was low quality (short reads ambiguously aligned and with many nucleotide differences) and (ii) the pattern of data presence would not allow for testing of their hypothesized taxonomic assignment (see supplementary data). When running the tree searches there were no exceedingly long branch lengths, and Blaberidae was monophyletic.
Trees were inferred for three alignments: 1) the full alignment ("265_Full"), 2) using only the 2 nd codon positions (correcting for noise; "265_2nd"), 3) low missing data alignment (correcting for relationships inferred from missing data patterns; "265_Reduced"). The latter alignment was created by only retaining nucleotide positions having data for 51 or more taxa (Phyloinformatica function trimAlign2 missingProportion = 0.60; Evangelista, 2019) ). The same partitioning, modelling and RAxML parameters as used above were applied to each analysis but with 100 independent tree runs (GTRGAMMA, -f d, -N 100). We inferred one more tree using the 265_Full alignment in IQTree (Nguyen et al., 2015) using partitions determined by PartionFinder2, models determined by IQTree (Kalyaanamoorthy et al., 2017) and the options:ninit 200 -nbest 10 -allnni -ntop 40 -wbt -wsl -wsr. These four trees (265_Full, 265_2nd, 265_Reduced, 265_Full_IQ) were later used as a baseline for establishing subsample convergence on a reasonable topology. We assessed support for the three RAxML trees by bootstrap resampling using the auto-MRE stopping criteria (60, 300 and 108 for the first three trees respectively) and calculating bootstrap frequencies and node certainty scores (Kobert et al., 2016) .