Dominic Evangelista

Manpreet France Thouzé

Kaur Kohli

Philippe Lopez

Frédéric Legendre

Manpreet Kaur Kohli

Topological support and data quality can only be assessed through multiple tests in reviewing Blattodea phylogeny

Introduction

Phylogenetic trees are increasingly used across biological disciplines and their interpretation can have deep implications (e.g. [START_REF] Beaulieu | Synthesizing phylogenetic knowledge for ecological research[END_REF][START_REF] Kozak | Phylogeny, ecology, and the origins of climate-richness relationships[END_REF][START_REF] Maganga | Ebola virus disease in the democratic Republic of Congo[END_REF]. Yet phylogenetic data are not straightforward and trees inferred from them can be conflicting. Traditional phylogenetic tests such as bootstrapping or tree comparisons are used to assess such conflicts (e.g. [START_REF] Blaimer | Phylogenomic methods outperform traditional multi-locus approaches in resolving deep evolutionary history: a case study of formicine ants[END_REF][START_REF] Garrison | Spider phylogenomics: untangling the Spider Tree of Life[END_REF][START_REF] Johnson | Phylogenomics resolves evolutionary relationships among ants, bees, and wasps[END_REF][START_REF] Kjer | Progress, pitfalls and parallel universes: a history of insect phylogenetics[END_REF][START_REF] Trautwein | Advances in insect phylogeny at the dawn of the postgenomic era[END_REF]. These tests, however, have their own limitations and biases and can sometimes give artefactual results [START_REF] Dell'ampio | Decisive data sets in phylogenomics: lessons from studies on the phylogenetic relationships of primarily wingless insects[END_REF]. Other tests and metrics have been used to complement traditional ones and provide a more complete picture of data signal and topological support (e.g. spectra analysis, saturation plots, likelihood mapping; [START_REF] Borowiec | Convergent evolution of the army ant syndrome and congruence in big-data phylogenetics[END_REF][START_REF] Dell'ampio | Decisive data sets in phylogenomics: lessons from studies on the phylogenetic relationships of primarily wingless insects[END_REF][START_REF] Kobert | Computing the internode certainty and related measures from partial gene trees[END_REF][START_REF] Wägele | Visualizing differences in phylogenetic information content of alignments and distinction of three classes of long-branch effects[END_REF]. In addition, they provide better guidance on crafting future datasets through locus choice, which is essential in the age of genomics because next-generation sequencing allows systematists to choose from a variety of loci whose quality varies (as in [START_REF] Borowiec | Convergent evolution of the army ant syndrome and congruence in big-data phylogenetics[END_REF][START_REF] Chen | Selecting question-specific genes to reduce incongruence in phylogenomics: a case study of jawed vertebrate backbone phylogeny[END_REF].

Pre-omics molecular phylogenetic approaches

Pre-omics molecular phylogenetic datasets typically comprise 2-7 markers (e.g. [START_REF] Bradler | A molecular phylogeny of Phasmatodea with emphasis on Necrosciinae, the most species-rich subfamily of stick insects[END_REF][START_REF] Chaboo | Origins and diversification of subsociality in leaf beetles (Coleoptera: Chrysomelidae: Cassidinae: Chrysomelinae)[END_REF][START_REF] Kambhampati | Evidence for sibling species in Cryptocercus punctulatus, the wood roach, from variation in mitochondrial DNA and karyotype[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Maekawa | The evolution of soil-burrowing cockroaches (Blattaria: Blaberidae) from wood-burrowing ancestors following an invasion of the latter from Asia into Australia[END_REF][START_REF] Marvaldi | Structural alignment of 18S and 28S rDNA sequences provides insights into phylogeny of Phytophaga (Coleoptera: Curculionoidea and Chrysomeloidea)[END_REF][START_REF] Murienne | Molecular data confirm family status for the Tryonicus-Lauraesilpha group (Insecta: Blattodea: Tryonicidae)[END_REF][START_REF] Song | 300 million years of diversification: elucidating the patterns of Orthopteran evolution based on comprehensive taxon and gene sampling[END_REF]. These markers usually come from both mitochondrial DNA (mtDNA) (e.g. [START_REF] Bradler | A molecular phylogeny of Phasmatodea with emphasis on Necrosciinae, the most species-rich subfamily of stick insects[END_REF][START_REF] Kambhampati | Phylogenetic relationship among cockroach families inferred from mitochondrial12S rRNA gene sequence[END_REF][START_REF] Kambhampati | Evidence for sibling species in Cryptocercus punctulatus, the wood roach, from variation in mitochondrial DNA and karyotype[END_REF][START_REF] Legendre | The phylogeny of termites (Dictyoptera: Isoptera) based on mitochondrial and nuclear markers: Implications for the evolution of the worker and pseudergate castes, and foraging behaviors[END_REF][START_REF] Mandal | Mitochondrial markers for identification and phylogenetic studies in insectsa review[END_REF][START_REF] Muraji | Phylogenetic analysis of water striders (Hemiptera: Gerroidea) based on partial sequences of mitochondrial and nuclear ribosomal RNA genes[END_REF][START_REF] Vogler | Exploring rate variation among and within sites in a densely sampled tree: species level phylogenetics of north american tiger beetles (genus cicindela)[END_REF][START_REF] Wiens | Discordant mitochondrial and nuclear gene phylogenies in emydid turtles: implications for speciation and conservation[END_REF] and nuclear DNA (nucDNA) (e.g. [START_REF] Bradler | A molecular phylogeny of Phasmatodea with emphasis on Necrosciinae, the most species-rich subfamily of stick insects[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Song | 300 million years of diversification: elucidating the patterns of Orthopteran evolution based on comprehensive taxon and gene sampling[END_REF] for animals, and T (2015) discussed that while the monophyly of Blaberoidea had strong morphological support, the molecular support was weak or lacking for one mtDNA marker (COII). Recent studies have demonstrated low support for evolutionary relationships among Blattoidea prior to the split between sub-social cockroaches (Cryptocercidae) and eusocial termites. This includes the position of Tryonicidae, Lamproblattidae, Blattidae [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF] and Anaplectidae [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF].

There is not just lack of support within studies, but lack of congruence among studies. For example, the competing molecular hypotheses about the position of Tryonicidae [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Murienne | Molecular data confirm family status for the Tryonicus-Lauraesilpha group (Insecta: Blattodea: Tryonicidae)[END_REF] confound arguments for [START_REF] Grandcolas | Reconstructing the past of Cryptocercus (Blattaria: Polyphagidae): phylogenetic histories and stories[END_REF][START_REF] Grandcolas | Systematique phylogenetique de la sous-famille des Tryonicinae (Dictyoptera, Blattaria, Blattidae)[END_REF] or against [START_REF] Klass | A phylogenetic Analysis of Dictyoptera (Insecta) based on morphological characters[END_REF] morphological similarity with Blattidae. Also, the existence of some xylophagous species of Tryonicidae [START_REF] Grandcolas | Systematique phylogenetique de la sous-famille des Tryonicinae (Dictyoptera, Blattaria, Blattidae)[END_REF] possibly suggests shared ancestry with Cryptocercidae and termites, which has sometimes been hypothesized with molecular data [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Murienne | Molecular data confirm family status for the Tryonicus-Lauraesilpha group (Insecta: Blattodea: Tryonicidae)[END_REF]. Similarly, putative morphological [START_REF] Klass | A phylogenetic Analysis of Dictyoptera (Insecta) based on morphological characters[END_REF] and behavioral [START_REF] Mckittrick | A contribution to the understanding of cockroach-termite affinities[END_REF] synapomorphies in Lamproblattidae are left ambiguous because of the unclear position of the family with respect to Cryptocercidae (e.g. [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF]. Sister group relationships at the first split in Blattodea are also controversial (e.g. [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF]. In Blaberoidea, nearly every combination of subfamilial relationships has been obtained in a molecular or combined analysis, the only constant being the monophyly of Blaberidae [START_REF] Bourguignon | Phylogeny of cockroaches (Insecta, Dictyoptera, Blattodea), with placement of aberrant taxa and exploration of out-group sampling[END_REF][START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Inward | Death of an order: a comprehensive molecular phylogenetic study confirms that termites are eusocial cockroaches[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF]. To summarize all these incongruences, [START_REF] Evangelista | Fossil calibrations for the cockroach phylogeny (Insecta, Dictyoptera, Blattodea), comments on the use of wings for their identification, and a redescription of the oldest Blaberidae[END_REF] reviewed eight phylogenetic studies of Blattodea and showed that the majority of the internal relationships were supported by less than half of the studies.

Using seven tests and metrics, we attempt to assess support of relationships across the Blattodea tree, determine strengths and weaknesses of the dataset contributing to various levels of support, and recommend future improvements based on these findings. We focus on several clades (full list in Appendix A) chosen for biological or methodological reasons. Thus, to better understand: (i) the evolutionary precedents to eusociality, we investigate the position of Tryonicidae Stolotermitidae. Note that we use "Anaplectidae" in a sense equal to that of Anaplectinae [START_REF] Roth | The cockroach genera Anaplecta, Anaplectella, Anaplectoidea, and Malaccina (Blattaria, Blattellidae: Anaplectinae and Blattellinae)[END_REF] and consider it a family-level lineage as did [START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF]. Conversely, we use "Nocticolinae" whereas other authors consider this a family, despite its position within Corydiidae s.l. [START_REF] Bourguignon | Phylogeny of cockroaches (Insecta, Dictyoptera, Blattodea), with placement of aberrant taxa and exploration of out-group sampling[END_REF][START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF].

chloroplast DNA for plants (e.g. [START_REF] Bremer | Molecular phylogenetic dating of asterid flowering plants shows early cretaceous diversification[END_REF][START_REF] Schulte | Phylogeny of Bromelioideae (Bromeliaceae) inferred from nuclear and plastid DNA loci reveals the evolution of the tank habit within the subfamily[END_REF]. Investigators often choose regions coding for ribosomal RNA (rDNA) because they contain a variety of fast and slow evolving sites and should be able to inform a phylogeny across multiple time-scales [START_REF] Hillis | Ribosomal DNA: molecular evolution and phylogenetic inference[END_REF][START_REF] Mandal | Mitochondrial markers for identification and phylogenetic studies in insectsa review[END_REF]. Indeed, all these types of data have allowed for large progresses in reconstructing the tree of life.

The small number of generic "toolkit" loci typical of pre-omics datasets necessitates consideration of their limitations. For instance, mtDNA markers are known to evolve quickly [START_REF] Simon | Evolution, weighting and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers[END_REF] and might be most appropriate to address recent evolutionary questions as opposed to deep ones [START_REF] Mandal | Mitochondrial markers for identification and phylogenetic studies in insectsa review[END_REF]. Also, mtDNA can have very different evolutionary histories from nuclear genes [START_REF] Fisher-Reid | What are the consequences of combining nuclear and mitochondrial data for phylogenetic analysis? Lessons from Plethodon salamanders and 13 other vertebrate clades[END_REF][START_REF] Kodandaramaiah | Deceptive single-locus taxonomy and phylogeography: Wolbachia-associated divergence in mitochondrial DNA is not reflected in morphology and nuclear markers in a butterfly species[END_REF][START_REF] Wiens | Discordant mitochondrial and nuclear gene phylogenies in emydid turtles: implications for speciation and conservation[END_REF] and phylogenies from nucDNA loci can be different from each other as well [START_REF] Lanier | Applying species-tree analyses to deep phylogenetic histories: challenges and potential suggested from a survey of empirical phylogenetic studies[END_REF]. Finally, both 18s + 28s and all mtDNA have limited biological independence. All mtDNA is transcribed at once [START_REF] Cameron | Insect mitochondrial genomics: implications for evolution and phylogeny[END_REF] and similarly the nuclear rDNA are encoded by the same transcription unit [START_REF] Lodish | Molecular Cell Biology[END_REF].

Since there is no ideal molecular marker, biases resulting from their analysis should be accounted for through phylogenetic tests and metrics. One common example in parsimony and maximum likelihood analyses is to calculate clade bootstrap support [START_REF] Pattengale | How many bootstrap replicates are necessary?[END_REF][START_REF] Soltis | Applying the bootstrap in phylogeny reconstruction[END_REF]. This test assesses internal alignment conflicts [START_REF] Soltis | Applying the bootstrap in phylogeny reconstruction[END_REF], which are dependent on signal conflicts among or within markers. Though, given the complicated nature of phylogenetic inference, a single test is far from universally informative. For instance, bootstrap values are highly dependent upon the number of taxa included [START_REF] Soltis | Applying the bootstrap in phylogeny reconstruction[END_REF] and can be biased by artefactual signal [START_REF] Dell'ampio | Decisive data sets in phylogenomics: lessons from studies on the phylogenetic relationships of primarily wingless insects[END_REF]. Complementary assays have been proposed to more thoroughly assess data quality and tree support. For example, loci information content and bias can be calculated prior to phylogenetic inference by examining: nucleotide compositional bias [START_REF] Song | When phylogenetic assumptions are violated: base compositional heterogeneity and among-site rate variation in beetle mitochondrial phylogenomics[END_REF], mutation saturation (Wenzel and Siddall, 1999), or information content [START_REF] Wägele | Visualizing differences in phylogenetic information content of alignments and distinction of three classes of long-branch effects[END_REF]. After tree reconstruction, molecular marker quality can be assessed through calculation of consistency indices [START_REF] Farris | The retention index and homoplasy excess[END_REF][START_REF] Klassen | Consistency indices and random data[END_REF], and locus saturation by branch length (e.g. [START_REF] Borowiec | Extracting phylogenetic signal and accounting for bias in whole-genome data sets supports the Ctenophora as sister to remaining Metazoa[END_REF]. Support for a given set of relationships can be estimated through statistical tests of alternative relationships [START_REF] Shimodaira | An approximately unbiased test of phylogenetic tree selection[END_REF] by comparing the effect of analysis permutations on tree reconstruction (e.g. [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Tang | Effects of phylogenetic reconstruction method on the robustness of species delimitation using single-locus data[END_REF], or measures of certainty and stability of trees and clades [START_REF] Kobert | Computing the internode certainty and related measures from partial gene trees[END_REF][START_REF] Legendre | One hundred years of instability in ensiferan relationships[END_REF]. Ideally, synthesis of all tests would allow one to conclude which relationships are strongly supported, which markers are mainly responsible for dark nodes (i.e. unresolved or unsupported relationships), and to suggest further improvements.

Molecular phylogenetics of Blattodea and dark areas of the tree

Using molecular data to decipher the relationships in Blattodea was revolutionary, as it was for other organisms [START_REF] Kjer | Progress, pitfalls and parallel universes: a history of insect phylogenetics[END_REF]. It helped shed light on the position of termites as sister to Cryptocercus and nested within the cockroach clade Blattoidea [START_REF] Inward | Death of an order: a comprehensive molecular phylogenetic study confirms that termites are eusocial cockroaches[END_REF][START_REF] Lo | Evidence for cocladogenesis between diverse dictyopteran lineages and their intracellular endosymbionts[END_REF], which had been controversial for decades (e.g. [START_REF] Grandcolas | The phylogeny of cockroach families: a cladistic appraisal of morpho-anatomical data[END_REF][START_REF] Grandcolas | Reconstructing the past of Cryptocercus (Blattaria: Polyphagidae): phylogenetic histories and stories[END_REF][START_REF] Klass | Morphological evidence on Blattarian phylogeny: "phylogenetic histories and stories" (Insecta, Dictyoptera)[END_REF][START_REF] Mckittrick | A contribution to the understanding of cockroach-termite affinities[END_REF][START_REF] Thorne | Phylogeny of Dictyoptera[END_REF]. It elucidated the co-evolutionary patterns of cockroaches and one of their most important endosymbionts [START_REF] Lo | Evidence for cocladogenesis between diverse dictyopteran lineages and their intracellular endosymbionts[END_REF]. It supported many other relationships and the monophyly of many groups (e.g. [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Legendre | Molecular phylogeny of Blaberidae (Dictyoptera, Blattodea) with implications for taxonomy and evolutionary studies[END_REF]).

Yet, current molecular datasets are unable to recover many blattodean relationships, despite strong taxon sampling (Fig. 1). For instance, one study highlights the difficulty in reconstructing the relationships within this clade using a set of six markers and 128 taxa [START_REF] Legendre | Molecular phylogeny of Blaberidae (Dictyoptera, Blattodea) with implications for taxonomy and evolutionary studies[END_REF]. In particular, some clades reconstructed with their full dataset were found only part of the time, or even rarely among individual gene trees [START_REF] Legendre | Molecular phylogeny of Blaberidae (Dictyoptera, Blattodea) with implications for taxonomy and evolutionary studies[END_REF]. Similarly, Djernaes et al. every partition.

A tree search was then implemented in RAxML [START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF]) on CIPRES [START_REF] Miller | Creating the CIPRES Science Gateway for Inference of[END_REF], running independent tree searches on 4 threads using the "-f o" option. Preliminary runs with 100, 20 and 10 random starting trees gave the same topology with only changes in branch lengths. Thus, we conducted the final tree search using 10 starting trees to minimize computation time (hereafter called "focal tree"; Appendix D).

Tests and metrics

Test 1. Testing the effect of removing taxa with low data completeness. Here, five trees were generated with different taxon reduction schemes (one full tree + four reduced trees). For each reduction scheme, taxa were eliminated based on gene sampling completeness. Starting with the full alignment (1geneMin) we removed taxa missing four or more genes (3geneMin), three or more genes (4geneMin), two or more genes (5geneMin), and one or more genes (6geneMin). The last and most complete of these alignments included 92 taxa, and still represented most major lineages of Blattodea.

Maximum likelihood trees were inferred as in Section 2.1, using each of the reduced alignments. To compare trees, we examine both the targeted clades (Appendix A; Section 1.2) and their support, which we calculated as a weighted average of the bootstrap clade support scores of the three clades comprising the relationship (i.e. the clade including the common ancestor, and the two daughter clades individually). The average weights the lowest supported clade by a factor of two. This is based on the following rationale. Given the tree (((A, B), (C, D)), E) the relationship of clade AB sister to clade CD relies on the support of three nodes: (A, B), (C, D) and ((A, B), (C, D)). If any of these nodes are weakly supported one could argue for collapsing them into a tritomy with the superior node [e.g. ((A), (B), (C, D)) or ((A, B), (C, D), (E))]. Collapsing any of these three nodes removes the sister relationship between AB and CD from the tree's topology. Therefore, the support for the relationship depends most heavily on which of the nodes has the lowest support (and is therefore the strongest rationale for collapsing the relationship).

Test 2. Testing the effect of reducing alignment completeness with constant taxon sampling. Six gene trees were generated starting from an alignment of full character sampling of 92 taxa (described in Section 2.2) and 92 alignment blocks were removed at random. Each alignment had subsequently more alignment blocks removed (i.e. a single random gene for a single random taxon; -0Blocks -92Blocks, -184Blocks, -276Blocks, -368Blocks, -460Blocks). Trees were inferred as described in Section 2.1 applied to each alignment. Results were recorded as for Test 1.

Test 3. Assessing node support through bootstrapping and internode certainty. Exhaustive tree searches for bootstrap pseudoreplicate trees using the autoMRE stopping criteria [START_REF] Pattengale | How many bootstrap replicates are necessary?[END_REF] were done for all the trees generated in the above tests. Additionally, relative tree certainty (RTC) and tree certainty (TC) scores were calculated for all Table 1 Molecular markers used in this analysis and descriptive information about their contribution to the total alignment. Information was taken from RAXML or manually calculated. Nucleotide composition is illustrated by proportion of each nucleotide, and relative bias is illustrated by base proportion standard deviation. and Lamproblattidae [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Murienne | Molecular data confirm family status for the Tryonicus-Lauraesilpha group (Insecta: Blattodea: Tryonicidae)[END_REF]; (ii) the early morphological evolution of Blattodea, we focus on Blattoidea, Corydioidea, and Blaberoidea [START_REF] Bourguignon | Phylogeny of cockroaches (Insecta, Dictyoptera, Blattodea), with placement of aberrant taxa and exploration of out-group sampling[END_REF][START_REF] Klass | A phylogenetic Analysis of Dictyoptera (Insecta) based on morphological characters[END_REF]. To assess potential artefactual reconstructions, we focus on taxa supported by: (iii) long branches (Nocticolinae and Pseudophyllodromiinae; [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF]; (iv) short branches (Blaberidae subfamilies; [START_REF] Legendre | Molecular phylogeny of Blaberidae (Dictyoptera, Blattodea) with implications for taxonomy and evolutionary studies[END_REF]; (v) rogue taxa (Diplopterinae; [START_REF] Legendre | Molecular phylogeny of Blaberidae (Dictyoptera, Blattodea) with implications for taxonomy and evolutionary studies[END_REF]; or (vi) taxa with incongruent position among studies (Blaberoidea and Blaberidae; [START_REF] Evangelista | Fossil calibrations for the cockroach phylogeny (Insecta, Dictyoptera, Blattodea), comments on the use of wings for their identification, and a redescription of the oldest Blaberidae[END_REF]. Finally, for the purpose of comparison, we focus on (vii) taxa with highly supported and congruent positions among studies (Mastotermitidae, Kalotermitidae, dampwood termites; Bourguignon et al., 2015;[START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF].

Methods

Phylogenetic methodology, tree reconstruction, and support

The alignment was composed of six genetic loci from 575 taxa. Sampling primarily included cockroaches and termites, but with a wide variety of outgroup taxa (Mantodea, Polyneoptera, and Palaeoptera). The six loci were: four mtDNA (12S rDNA, 16S rDNA, COI, COII), and two nucDNA (18S rDNA, 28S rDNA). We obtained the dataset from [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF] with 65 additional sequences from Genbank (Appendix B), and with many Mantodea removed such that we could focus on relationships in Blattodea.

Starting with the alignment from [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF] rDNA sequences were manually re-aligned to a structural model from [START_REF] Ware | Relationships among the major lineages of Dictyoptera: the effect of outgroup selection on dictyopteran tree topology[END_REF] following [START_REF] Kjer | Structural and evolutionary considerations for multiple sequence alignment of RNA, and the challenges for algorithms that ignore them[END_REF] (Appendix C). This eliminated unalignable hypervariable loop regions but retained hypervariable loop regions for some taxa (mantises and termites separately) that were indeed alignable. Forcing this pattern of missing data could enforce monophyly of the three groups. Yet, the relationships among the three "orders" of Dictyoptera have been settled for at least 10 years [START_REF] Deitz | Phylogeny of the Dictyoptera Re-examined (Insecta)[END_REF][START_REF] Inward | Death of an order: a comprehensive molecular phylogenetic study confirms that termites are eusocial cockroaches[END_REF][START_REF] Klass | Morphological evidence on Blattarian phylogeny: "phylogenetic histories and stories" (Insecta, Dictyoptera)[END_REF][START_REF] Lo | Evidence for cocladogenesis between diverse dictyopteran lineages and their intracellular endosymbionts[END_REF] and deviation from these relationships would be a clear sign of error in the analysis. The missing data patterns could still be inflating support for Mantodea or termites, which should be considered when interpreting the results. The retained hypervariable regions were aligned with gaps for all other taxa. The final alignment was 7676 nucleotides long (Appendix D; Table 1), whereas it was about 10 kb in [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF].

Data blocks were defined as the boundaries of loci, hypervariable regions and codon positions (for the protein coding genes: COI and COII). In PartitionFinder [START_REF] Lanfear | PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses[END_REF], using these block definitions, output models for RAxML, and the BIC optimality criterion, the following partitioning scheme was determined: 18S + 28S, 12S + 16S (including hypervariable regions), 28S "termite" hypervariable region, 28S "mantis" hypervariable region, COI + COII codon 1, COI + COII codon 2, and COI + COII codon 3. The model GTR + G was chosen for bootstrap sets. TC scores are raw sums of all node IC's and RTC are those values normalized in relation to the number of branches whose removal induces a non-trivial bipartition [START_REF] Kobert | Computing the internode certainty and related measures from partial gene trees[END_REF]. Certainty values are a measure of diversity of alternative topologies among a set of trees [START_REF] Kobert | Computing the internode certainty and related measures from partial gene trees[END_REF]. Thus certainty values differ from bootstrap frequencies, which only indicate abundance of the reported topology. Using certainty scores, one can differentiate between scenarios where the focal tree topology is the most common topology, or less common than another tree topology within the pseudoreplicate set. For example, a clade with 40% bootstrap support could be the most commonly recovered clade (IC > 0) or less commonly recovered than another clade (IC < 0). Certainty scores were computed using the RAxML "-L MR" options on the bootstrapped pseudoreplicate sets [START_REF] Kobert | Computing the internode certainty and related measures from partial gene trees[END_REF]. IC-all (ICA) scores, which take into account all recovered clades instead of just the two most abundant, were calculated for the fullysampled tree (i.e. the "focal tree" obtained as in Section 2.1) using the "-f i" options.

Test 4. Assessing statistical support for alternative topologies. The approximately unbiased (AU) test was used to determine whether the full dataset statistically rules out alternative topologies [START_REF] Shimodaira | An approximately unbiased test of phylogenetic tree selection[END_REF]. The AU test is a tree choice test that corrects for the selection and comparison biases of similar tests [START_REF] Shimodaira | An approximately unbiased test of phylogenetic tree selection[END_REF]. Yet, all the tests are based on the same underlying principle of assessing variance in the likelihoods of alternative trees inferred from bootstrapped tree sets and assigning p-values to topologies based on those variances [START_REF] Kishino | Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in Hominoidea[END_REF]). The "focal tree" described in Section 2.1 was used as a starting tree. The tree was manually edited in Mesquite (Maddison and Maddison, 2017) to create alternative input trees and to remove branch lengths. AU tests were done in the software IQ-Tree [START_REF] Nguyen | IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF] using 10,000 resampled estimated log-likelihood (RELL) bootstraps and the partitioning scheme define in Section 2.1.

Test 5. Visualizing split signal and conflict. The purpose of Test 5 was to evaluate signal and conflict within the data (as opposed to conflict in inferences from the data as in Test 3). Optimized phylogenetic trees inferred from datasets show relationships heavily favored by emergent signal in the dataset but do not show alternative relationships partly supported by the data, as phylogenetic networks can. However, phylogenetic networks cannot be entirely represented in a 2-dimensional space [START_REF] Wägele | Visualizing differences in phylogenetic information content of alignments and distinction of three classes of long-branch effects[END_REF] making them difficult to assess qualitatively. Lento plots, also known as spectral plots/analysis, are a simple way of representing the abundance of nucleotide characters supporting specific splits in a phylogenetic network, splits that may conflict with one another. Bars, whose height represent the number of characters supporting a given split, are plotted for any given number of splits in descending order of support. Using the software SAMS [START_REF] Wägele | Visualizing differences in phylogenetic information content of alignments and distinction of three classes of long-branch effects[END_REF] split noise and conflict are also indicated. Visually inspecting these plots in parallel with a tree can give evidence for longbranch attraction phenomena, signal/noise ratio and information content of the data [START_REF] Wägele | Visualizing differences in phylogenetic information content of alignments and distinction of three classes of long-branch effects[END_REF]. Since alignments with high proportions of missing data can be problematic in SAMS, the alignments were trimmed both vertically and horizontally using a custom script in Mathematica (Appendix E) before generating Lento plots. SAMS was used with default parameters (all max noise = 0.25, consensus threshold = 0.5, pairwise comparison, analyze 150 occurring splits) except gaps were treated as missing data and total (tree) split compatibility was visualized.

Test 6. Assessing locus saturation. Quantifying saturation in molecular markers informs on their phylogenetic utility and the time-scale of their optimal utility [START_REF] Salichos | Inferring ancient divergences requires genes with strong phylogenetic signals[END_REF]. Mutation saturation was assessed in two ways: plotted raw pairwise genetic distances against GTR model corrected pairwise genetic distances, and plotted raw pairwise genetic distances against pairwise tree branch-length distances using the "focal tree" with all taxa. Genetic distances were calculated in PAUP* [START_REF] Swofford | PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods)[END_REF] and plotted in Mathematica [START_REF] Research | Mathematica Edition: Version 9[END_REF]. Tree distances were calculated in R using the "cophenetic" function [START_REF] Sneath | Numerical Taxonomy: The Principles and Practice of Numerical Classification[END_REF]. Input data can be found in Appendix F.

Test 7. Assessing homoplasy in molecular markers. Here we distinguished markers evolving mostly parsimoniously with tree structure from more homoplastic markers using consistency index (CI) and retention index (RI) scores. Though homoplastic characters are not equivalent to noisy characters, extremely homoplastic characters contribute more to tree error than other characters [START_REF] Breinholt | Phylotranscriptomics: saturated third codon positions radically influence the estimation of trees based on next-gen data[END_REF]Wenzel and Siddall, 1999; but see [START_REF] Vogler | Exploring rate variation among and within sites in a densely sampled tree: species level phylogenetics of north american tiger beetles (genus cicindela)[END_REF]. Scores were calculated for alignment regions in R using the package "phangorn" [START_REF] Schliep | phangorn: phylogenetic analysis in R[END_REF] and the "focal tree".

Results

Alternative topologies and support (Tests 1-4)

Taxon reduced tree inferences (Test 1; Appendices G & H) were mostly congruent with the "focal tree" (Fig. 2; Appendix D). Among the relationships examined, the position of Anaplectinae and the relationships within Blaberidae differed the most. The results of the gene reduction schemes (Test 2; Appendix H) showed the positions of termite lineages and the sister relationship between Blaberidae and Pseudophyllodromiinae + Ectobiinae are largely robust to perturbation.

The relative tree certainty scores (Table 2; Appendices G & H) correlated with mean bootstrap values (Test 3). The gene reduction tree with 184 data block removals (-184Blocks) had comparatively less certainty relative to its bootstrap scores. In the "focal tree", the most uncertain relationships were the sister group to Blaberidae, and the position of Blattellinae + Nyctiborinae (Appendix H; Fig. 2). Both of these relationships had a negative weighted-average internode certainty, indicating that alternative topologies were more common than the relationships recovered in the "focal tree". The relationships with the highest internode certainty relative to bootstrap support were the positions of Lamproblatta and Tryonicus.

Test 4 (Appendices H & I) showed that even relationships recovered consistently in Test 1 (and with relatively high bootstrap and internode certainty support values) were not unambiguously supported by the data. Multiple alternative topologies, some of which were never recovered in any of the alignment permutation analyses (Tests 1 and 2), could not be rejected. The alternative topology with Ectobiinae sister to remaining Blaberoidea was never recovered in any of our analyses but actually produced higher log-likelihood (logL) score in the AU test (Appendix I).

Signal and conflicts (Tests 5-7)

Lento plots (Test 5; Fig. 3; Appendix E) showed that the majority of nucleotide support for splits (signal) are in 18S and 28S, while COI + COII had moderate signal, and 12S + 16S had low signal. Intramarker conflict was highest in the protein coding genes, and lowest among the nuclear rDNAs. The least conflict occurred in the combined alignment, with only 46% of the first 50 splits conflicting and none of the top 10 splits conflicting. In general, the more markers considered the less internal conflict in the alignment. However, the two nuclear markers together contain more signal and less internal conflict than all four mitochondrial markers (Appendix E).

The two data quality analyses (Tests 6 & 7; Figs. 4-6; Appendix F) both showed that COI, COII and 12S seem to have a more limited phylogenetic utility in Blattodea than 16S, 18S and 28S. 12S and COII in particular had both low consistency (Fig. 6) and high saturation at deep branches (Fig. 5). Consideration by codon position showed that first codon positions in COII and third codon positions in both COI and COII are highly saturated (Fig. 4).

Discussion

Support for the phylogeny of Blattodea

The phylogenetic hypothesis we recovered from the full dataset (1GeneMin tree; Fig. 2; Appendix D) mostly agreed with the topology of [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF]. The tree in [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF] has higher bootstrap support values comparatively (Appendix H). One notable difference though is in the support for Panchlora as sister to the remaining Blaberidae, which is strongly supported in [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF] but given nearly negligible support here. Most other differences between the two trees are relationships within Blaberidae. Also, since [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF] did not include most of the Anaplectidae included here, it should be noted that this clade was recovered as polyphyletic, with the majority of the family as highly volatile within Blattoidea. One individual ("Anaplecta sp.") fell in Pseudophyllodromiinae and it must be investigated further whether this placement results from a specimen misidentification or a poorly defined taxon that needs revision.

Alignment permutations tests (Appendix H) helped identify some unstable relationships in the tree. Other studies removing taxa and affecting alignment completeness have shown similar variations in tree topology [START_REF] Chen | Selecting question-specific genes to reduce incongruence in phylogenomics: a case study of jawed vertebrate backbone phylogeny[END_REF]. The position of Anaplectidae (except "Anaplecta sp.") and the lineages of Blaberidae were the most affected by alignment permutations. Most of the other relationships were more stable, which contradicts their volatility when comparing among [START_REF] Evangelista | Fossil calibrations for the cockroach phylogeny (Insecta, Dictyoptera, Blattodea), comments on the use of wings for their identification, and a redescription of the oldest Blaberidae[END_REF] or within [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF] past studies. Similarly, the phylogenetic relationships reconstructed here are more stable overall than those [START_REF] Ware | Relationships among the major lineages of Dictyoptera: the effect of outgroup selection on dictyopteran tree topology[END_REF] found when experimenting with taxon removal/addition in Dictyoptera phylogenetics. The higher stability here [and in [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF] because their trees were mostly congruent] is likely due to the larger overall sampling and use of more diverse outgroups than in [START_REF] Ware | Relationships among the major lineages of Dictyoptera: the effect of outgroup selection on dictyopteran tree topology[END_REF]. This is in agreement with theoretical and empirical studies underlining the importance of large taxon sampling [START_REF] Heath | Taxon sampling and the accuracy of phylogenetic analyses[END_REF][START_REF] Zwickl | Increased taxon sampling greatly reduces phylogenetic error[END_REF].

The approximately unbiased (AU) test results mostly corroborated topological support (both strong and weak) from previous tests, but gave some unique insight (Tests 4; Appendix I). The AU tests did not reject any of the topologies that were volatile in the alternative analyses (see Section 2.2, Test 1; Appendix H). Similarly, AU tests did not reject alternative positions for Corydioidea or Nyctiborinae + Blattellinae, yet their placement was not volatile according to Tests 1 and 2. This illustrates that even if the phylogenetic analysis of a dataset results in a best unique tree, this dataset is not necessarily statistically incompatible with alternative topologies [START_REF] Buddenhagen | Anchored Phylogenomics of Angiosperms I: Assessing the Robustness of Phylogenetic Estimates[END_REF]. Importantly, in

Table 2

Alignments used in various analyses. 10 alignments used in our main analysis (1geneMin) and subsequent tests of taxon (Test 1) and gene (Test 2) reduction. First 5 alignments are the full alignment with subsequent reductions of taxa by removal of taxa with increasing amounts of missing data (e.g. 3geneMin equals 1geneMin but with all taxa missing more than 3 markers removed). Last 5 alignments are the 6geneMin (alignment with all 92 taxa having all 6 markers present) with the specified number of alignment blocks removed (e.g. -92Blocks alignment has 92 random partitions from random taxa removed). Various alignment statistics are taken from RAXML or calculated in Mathematica 9. a dataset with heterogeneous partition quality and missing data patterns (Table 1) the AU test may lack the statistical power to favor one topology over another, which likely explain why this test was very conservative in comparison to the others. As such, we should consider topologies not uniquely supported by the AU test as given tentatively low support, and other test results should take precedence.

Some relationships we focused on (Appendix H) had negative ICA values. In these cases, the relationships in the 1GeneMin tree were less frequent than other relationships recovered among the pseudo-replicate trees inferred from the same alignment. This reveals conflicts, incongruent phylogenetic signal, among or within the phylogenetic markers used, which is common in datasets with a small number of markers Fig. 3. Splits support spectrograms for alignments. We focus on alignments for 18S and 28S together, 12S and 16S together, COI and COII together, all mtDNA together and a full-length alignment (4geneMin alignment). Horizontal axes are tree splits (bipartitions) supported by the data. Vertical axis shows how many nucleotides support each split. Orange and tan bars indicate noisy data supporting the splits. Grey bars indicate a split that is supported by the alignment but is conflicting with another split that has higher nucleotide support. Only support for the first 50 splits are shown. Pairwise genetic distance with GTR model corrected distance on the horizontal axis and uncorrected genetic distances on the vertical axis. A dashed line is included so that the genetic distances can be compared with parity. Parity would indicate that there is no locus saturation (no homoplasy, and high signal) and high deviation from parity indicates that there is high saturation (high homoplasy and noise). For consistency, all plot ranges are given as between 0 and 1. Plots indicate that COI, COII and 12S are the most saturated and 16S is moderately saturated. 18S and 28S are not saturated but have low genetic distances overall, indicating high conservation. For the protein coding genes, COII first codon positions and COI and COII 3rd codon positions are highly saturated. [START_REF] Narechania | Random Addition Concatenation Analysis: a novel approach to the exploration of phylogenomic signal reveals strong agreement between core and shell genomic partitions in the cyanobacteria[END_REF]. The reasons for this pattern are unclear, however we suspect it might have to do with internal non-independence of rRNA's (due to paired-mutations in self-bonded regions), patterns of missing data [START_REF] Dell'ampio | Decisive data sets in phylogenomics: lessons from studies on the phylogenetic relationships of primarily wingless insects[END_REF], or heterogeneous partition quality (Table 1). The low scores overall and the heterogeneity among pseudoreplicate trees indicates that: more markers are needed to stabilize the tree topology [START_REF] Narechania | Random Addition Concatenation Analysis: a novel approach to the exploration of phylogenomic signal reveals strong agreement between core and shell genomic partitions in the cyanobacteria[END_REF]) and dilute potential biases; better modelling strategies are needed (e.g. [START_REF] Letsch | Potential pitfalls of modelling ribosomal RNA data in phylogenetic tree reconstruction: evidence from case studies in the Metazoa[END_REF][START_REF] Zhang | Phylotranscriptomic analysis based on coalescence was less influenced by the evolving rates and the number of genes: a case study in ericales[END_REF]; or missing data patterns should be filled in through new sequencing. In parallel, signal analyses must be conducted to differentiate the markers supporting the relationships from those with highly conflicting signal.

Signal and conflicts in molecular markers

Lento-plots (Fig. 3; Appendix E) give insights into three classes of long-branch effects [START_REF] Wägele | Visualizing differences in phylogenetic information content of alignments and distinction of three classes of long-branch effects[END_REF]: symplesiomorphy driven (class 1), signal-erosion driven (class 2), and homoplasy driven (class 3). We did not find evidence of class 1 in the data because there are few splits where large clades are lumped with outgroups. When ingroup taxa matched with outgroups they also tended to group incorrectly with a few other ingroup taxa supported by long branches. These taxa were: Latindiinae spp., Nocticolidae sp., Anaplecta sp., Supella longipalpa, Dendroblatta sp. and Isoldaia sp. Also, the taxa clustering with outgroups were individual species or members of the same genus, suggesting that the states driving the long branch effect are derived. These taxa are probably clustering under a class 2 long branch effect derived from saturated nucleotide positions (Figs. 4 and5; Test 6; as in [START_REF] Omilian | Rate acceleration and long-branch attraction in a conserved gene of cryptic daphniid (Crustacea) species[END_REF]. Lento plots for 3rd codon positions clearly indicate class 2 long branch effects and all splits are polyphyletic when mapped onto the "focal tree". Class 3 effects would result from homoplastic sites contributing to conflicting signal (as in the data of [START_REF] Crandall | On the phylogenetic positioning of the South American freshwater crayfish genera (Decapoda: Parastacidae)[END_REF], but conflict was low in the combined data Lento plots (Fig. 3) so presence of a class 3 effect seems less likely here.

Interpretation of splits analysis methods (such as Lento plots; e.g. [START_REF] Wägele | Phylogenetic support values are not necessarily informative: the case of the Serialia hypothesis (a mollusk phylogeny)[END_REF], and integration with other analytics, can also give insight into the phylogenetic value of specific molecular markers. 18S and 28S provided the highest support to splits (Fig. 3) despite being conserved overall (Fig. 5) and having the lowest density of unique patterns (Table 1). Yet, including mtDNA markers, which had low phylogenetic signal on their own (Fig. 3), did improve split support and lessened internal conflict (Fig. 3). Given that many of the nodes on the final tree have low support (Appendix H) it is unclear how the addition of the mtDNA markers, which were comparatively more conflicting, affected them. In fact, given the much stronger nucleotide compositional bias in the mtDNA (Table 1) it is possible that nodes supported ambiguously by nucDNA markers would have been swayed by mtDNA compositional bias. Additional tests are needed to further assess these effects (i.e. site rate variation, compositional heterogeneity and violation of sequence evolution assumptions) from individual markers [START_REF] Song | When phylogenetic assumptions are violated: base compositional heterogeneity and among-site rate variation in beetle mitochondrial phylogenomics[END_REF].

Combining several tests for a better understanding of phylogenetic hypotheses

The assays employed here proved useful for elucidating some dark regions in the Blattodea tree and for understanding the discordance among previous studies (as illustrated in [START_REF] Evangelista | Fossil calibrations for the cockroach phylogeny (Insecta, Dictyoptera, Blattodea), comments on the use of wings for their identification, and a redescription of the oldest Blaberidae[END_REF]. However, no single test or metric was sufficient to identify either all areas of the tree that are in need of new data (Fig. 7) or the relative contribution of molecular markers to tree instability [START_REF] Chiapella | Fishing for significance in phylogenies: too many alternatives for the same outcome, or an appeal to journal editors[END_REF][START_REF] Vogler | Exploring rate variation among and within sites in a densely sampled tree: species level phylogenetics of north american tiger beetles (genus cicindela)[END_REF]. This is unsurprising given that phylogenetic uncertainty comes from a multitude of sources and that no single test or metric is able to identify them all, which underline the importance of confronting them.

A dark region identified using multiple tests is the position of Corydioidea (= Corydiidae s.l. = Polyphagoidea). On one hand, alignment permutation tests (Tests 1 and 2) consistently recovered Corydioidea as sister to Blattoidea (Appendix H). There were two exceptions (-92Blocks and -460Blocks trees) but this was still one of the more stable relationships and it had moderately high support (weighted bootstrap in the "focal tree" = 41). On the other hand, weighted ICA for this relationship was near 0 (Test 3), and AU tests (Appendices H and I) showed that the alternative of Corydioidea as sister to all of Blattodea is also a plausible topology (Test 4). [START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF] claimed to have recovered Corydioidea + all other Blattodea with strong support, but their assessment of support was limited to bootstrap frequency and Bayesian posterior probability. We have not reanalyzed their data so we do not know if other tests would agree with their topology. Given the discordance between our tests, we doubt the certainty of either proposed relationship. To summarize, we found the phylogenetic position of Corydioidea as stable and with moderate bootstrap support, but with very low ICA scores and no statistical support in the AU test. Using only one of these tests would have resulted in drastic conclusions (in one way or another), which would have been premature.

Similarly, the set of analyses conducted here brought mixed results for the position of Tryonicidae. This lineage is thought to be important in understanding the evolutionary past of termites because [START_REF] Murienne | Molecular data confirm family status for the Tryonicus-Lauraesilpha group (Insecta: Blattodea: Tryonicidae)[END_REF] and [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF] recovered Tryonicidae as sister to (Cryptocercidae + termites). This relationship was never recovered in any of the alignment permutations (Tests 1 and 2), and IC values were relatively high in the 1GeneMin tree (Test 3). This set of results is congruent with [START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF]. AU-tests, however, did reject some alternative topologies but not all (Test 4; Appendix I) and Lento plots gave no indication of a potential long branch effect with this taxon (Test 5; Appendix E). The different positions of Tryonicidae in past studies could be due to: missing Anaplectinae and Lamproblattidae and a relatively limited taxon sampling in [START_REF] Murienne | Molecular data confirm family status for the Tryonicus-Lauraesilpha group (Insecta: Blattodea: Tryonicidae)[END_REF]; inclusion of morphological data, different alignment partitioning schemes or tree reconstruction methods in [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF].

Suggestions and future directions

The tree of life will undoubtedly become more fully resolved as transcriptomic and genomic data are increasingly collected from diverse taxa (e.g. [START_REF] Blaimer | Phylogenomic methods outperform traditional multi-locus approaches in resolving deep evolutionary history: a case study of formicine ants[END_REF]. "Big data" are not yet available for every group of organisms and some progress could still be made with pre-omics datasets, such as the six molecular markers studied here. However, we must analyze how these molecular data contribute to phylogenetic inferences using tests that will help identify both current dark areas in phylogenetic trees and the target profile of future markers to decipher these unresolved areas. An effective strategy would be to mine genomic studies (e.g. Granados [START_REF] Granados Mendoza | A genome-scale mining strategy for recovering novel rapidly-evolving nuclear single-copy genes for addressing shallow-scale phylogenetics in Hydrangea[END_REF][START_REF] Misof | Selecting informative subsets of sparse supermatrices increases the chance to find correct trees[END_REF][START_REF] Misof | A priori assessment of data quality in molecular phylogenetics[END_REF] for new, more informative, and more independent markers (as in [START_REF] Chen | Selecting question-specific genes to reduce incongruence in phylogenomics: a case study of jawed vertebrate backbone phylogeny[END_REF]. Some of the tests discussed here should be used to tailor locus choice to fit the needs of a particular phylogenetic question. For instance, our evidence of class 2 or 3 long branch effects could be remedied by targeting characters more conserved than the saturated or low consistency ones found here. This aim is not easily achieved though, as the interaction of mutation rate, timescale and phylogenetic context is rarely straightforward. The phylogeny of Blattodea illustrates how a given set of molecular markers provide excellent resolution at a certain timescale for one clade (i.e. termites) but low resolution for similar time-scales on nearby clades (i.e. Blaberoidea; considering dates from [START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF]. Though, it is good to keep in mind that some nodes in phylogenies are inherently difficult to infer [START_REF] Narechania | Random Addition Concatenation Analysis: a novel approach to the exploration of phylogenomic signal reveals strong agreement between core and shell genomic partitions in the cyanobacteria[END_REF][START_REF] Salichos | Inferring ancient divergences requires genes with strong phylogenetic signals[END_REF] and identifying which nodes these are is key to our understanding of evolutionary histories [START_REF] Whitfield | Deciphering ancient rapid radiations[END_REF].

For numerous taxa, pre-omics datasets are still useful and, to benefit

Conclusions

Molecular datasets are ever-increasing in size so that interpreting the relationship between input and output is becoming more difficult. However, comprehensive tests can disentangle dataset weaknesses and strengths. Assessments of data quality can be done both before and after tree inference but combining the results gives the most meaningful assessment. Multiple lines of evidence are needed to fully identify tree support and data robustness.

In the phylogeny of Blattodea, only AU tests and ICA scores, and not bootstrap scores or alignment permutations, illustrated the lack of support for the placement of Corydioidea. Split support analysis showed that few deep relationships in the ingroup had visible character support. As such, tree inference and the calculation of support values are vulnerable to data biases, signal saturation, conflicting signal and other issues. Indeed, our contradictory results and alternative topologies could not be rejected using the information in the alignment. More optimistically, although there was evidence of long-branch attraction in the split support plots, long-branch taxa appear correctly placed in the even further from them, a number of steps could be considered prior to tree reconstruction (Misof et al., 2014b). With our dataset, structural alignment of rRNA had minimal effect on tree reconstruction. Still, not aligning with respect to structural homology (as in Bourguignon et al., 2015;[START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF] could introduce error associated with phenetic inference [START_REF] Kjer | Opinions on multiple sequence alignment, and an empirical comparison of repeatability and accuracy between POY and structural alignment[END_REF], which could reduce tree accuracy and artificially inflate support values. Further improvements could be made by down-weighting rRNA stems, and modelling stem and loop regions separately, but only if these regions do not violate model assumptions [START_REF] Letsch | Potential pitfalls of modelling ribosomal RNA data in phylogenetic tree reconstruction: evidence from case studies in the Metazoa[END_REF]. Finally, omitting saturated nucleotide regions could improve tree reconstruction [START_REF] Breinholt | Phylotranscriptomics: saturated third codon positions radically influence the estimation of trees based on next-gen data[END_REF]Wenzel and Siddall, 1999), although sometimes highly homoplastic sites provide "hidden support" to nodes, even deep ones [START_REF] Vogler | Exploring rate variation among and within sites in a densely sampled tree: species level phylogenetics of north american tiger beetles (genus cicindela)[END_REF]Wenzel and Siddall, 1999). Saturated sites could be removed by alignment masking software (e.g. [START_REF] Castresana | Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis[END_REF][START_REF] Wu | Accounting for alignment uncertainty in phylogenomics[END_REF] or in protein coding genes analyzed with a "R-Y" coding. Running multiple tests prior to tree reconstruction might help in making educated guesses as to the best strategy to follow for a given dataset. Relationships recovered consistently in taxon reduction tests but that also have low node support (such as the position of Corydioidea) could be better resolved by improving alignment completeness with the current molecular markers. Conversely, highly volatile taxa (like Blaberidae lineages) will probably not benefit from additional sequencing to improve completeness of these six markers, since multiple individuals representing major lineages have highly complete data but are still very volatile. Anaplectidae, although volatile, currently has such a low data completeness that further sequencing might provide resolution.

Our results show that assessing the strength of a phylogenetic hypothesis requires integration of tests and metrics. Tree independent methods of visualizing support, such as Lento Plots (via SAMS; [START_REF] Wägele | Visualizing differences in phylogenetic information content of alignments and distinction of three classes of long-branch effects[END_REF] or split-networks [START_REF] Bryant | Neighbor-net: an agglomerative method for the construction of phylogenetic networks[END_REF][START_REF] Wägele | Phylogenetic support values are not necessarily informative: the case of the Serialia hypothesis (a mollusk phylogeny)[END_REF], directly measure support for relationships and can be used to identify data features causing long branch attraction [START_REF] Wägele | Visualizing differences in phylogenetic information content of alignments and distinction of three classes of long-branch effects[END_REF]. Bootstrap resampling can give estimates of conflicting signal for clades [START_REF] Kobert | Computing the internode certainty and related measures from partial gene trees[END_REF][START_REF] Soltis | Applying the bootstrap in phylogeny reconstruction[END_REF] but bootstrap frequencies alone may not be fully indicate the extent of conflict and can be complemented with certainty scores [START_REF] Kobert | Computing the internode certainty and related measures from partial gene trees[END_REF]. Similarly, statistical tests of topology (i.e. AU test) can test alternatives not recovered in the most likely tree [START_REF] Shimodaira | An approximately unbiased test of phylogenetic tree selection[END_REF] but their power is limited by data homogeneity. These multifaceted approaches improve understanding of the dependence of the results on sampling [START_REF] Narechania | Random Addition Concatenation Analysis: a novel approach to the exploration of phylogenomic signal reveals strong agreement between core and shell genomic partitions in the cyanobacteria[END_REF][START_REF] Ware | Relationships among the major lineages of Dictyoptera: the effect of outgroup selection on dictyopteran tree topology[END_REF], modelling [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF], alignment bias [START_REF] Dell'ampio | Decisive data sets in phylogenomics: lessons from studies on the phylogenetic relationships of primarily wingless insects[END_REF], or inference method [START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Zhang | Phylotranscriptomic analysis based on coalescence was less influenced by the evolving rates and the number of genes: a case study in ericales[END_REF].

"focal tree". Yet, the long-branch signal may be contributing to low tree support.

Fig. 1 .

 1 Fig. 1. Phylogenetic context for the present study of Blattodea. Solid edges represent relationships considered well-supported among all previous molecular and morphological studies on Blattodea and Dictyoptera phylogenetics. Dashed edges in polytomies represent clades whose position is uncertain. HAS termite clade consists of Hodotermitidae, Archotermopsidae andStolotermitidae. Note that we use "Anaplectidae" in a sense equal to that of Anaplectinae[START_REF] Roth | The cockroach genera Anaplecta, Anaplectella, Anaplectoidea, and Malaccina (Blattaria, Blattellidae: Anaplectinae and Blattellinae)[END_REF] and consider it a family-level lineage as did[START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF]. Conversely, we use "Nocticolinae" whereas other authors consider this a family, despite its position within Corydiidae s.l.[START_REF] Bourguignon | Phylogeny of cockroaches (Insecta, Dictyoptera, Blattodea), with placement of aberrant taxa and exploration of out-group sampling[END_REF][START_REF] Djernaes | Identifying possible sister groups of Cryptocercidae+Isoptera: a combined molecular and morphological phylogeny of Dictyoptera[END_REF][START_REF] Legendre | Phylogeny of Dictyoptera: dating the origin of cockroaches, praying mantises and termites with molecular data and controlled fossil evidence[END_REF][START_REF] Wang | Reconstructing the phylogeny of Blattodea: robust support for interfamilial relationships and major clades[END_REF].

Fig

  Fig. 2. Focal tree topology (1geneMin tree) for major lineages in Blattodea and support. Tree was inferred in RAXML from molecular data for six loci obtained for 575 taxa. Support values are bootstrap frequencies and internode certainty values. We see strong support for the monophyly of Blattodea, Blaberoidea, termites + Cryptocercidae, and termites; moderate support for Lamproblattidae + Cryptocercidae + termites, Tryonicidae + Blattidae, and Blaberidae; low support for Anaplectidae + sister group, Blattoidea, Corydioidea, and most relationships in Blaberoidea. Full tree with support values is available in Appendix D. See Fig. 1 caption for family-name usages.

Fig. 4 .

 4 Fig. 4. Overall locus saturation overall.Pairwise genetic distance with GTR model corrected distance on the horizontal axis and uncorrected genetic distances on the vertical axis. A dashed line is included so that the genetic distances can be compared with parity. Parity would indicate that there is no locus saturation (no homoplasy, and high signal) and high deviation from parity indicates that there is high saturation (high homoplasy and noise). For consistency, all plot ranges are given as between 0 and 1. Plots indicate that COI, COII and 12S are the most saturated and 16S is moderately saturated. 18S and 28S are not saturated but have low genetic distances overall, indicating high conservation. For the protein coding genes, COII first codon positions and COI and COII 3rd codon positions are highly saturated.

Fig. 5 .

 5 Fig. 5. Locus saturation per time-scale. Pairwise branch distance on the horizontal axis plotted against pairwise GTR model corrected distance on the vertical axis. Trend lines calculated from a best fit linear regression. Equations for each trend line are as follows: 12S, 0.048 + 0.077x; 16S, 0.063 + 0.053x; 18S, -0.019 + 0.023x; 28 s, 0.01 + 0.035x; COI, 0.155 + 0.040x; COII, 0.160 + 0.084x. Plots show that 18S and 28S are more conserved at all ages and have phylogenetic utility at deep phylogenetic scales. 16S and COI are moderately fast evolving. 12S and COII are the least conserved and may be limited in their use at deep phylogenetic scales.

Fig. 6 .

 6 Fig.6. Homoplasy among loci. Plots of consistency index (CI) and retention index (RI) for the six loci considered in this study. Protein coding genes are also assessed for each codon position. CI scores show greater differences than RI scores. CI scores for 16S, 18S, and 28S show that these loci have low levels of homoplasy. Contrarily, 12S, COI and COII are highly homoplasious. However, COI and COII second codon positions and COI first codon position appear less homoplastic, while third codon positions for both and first codon positions for COII appear highly homoplastic.

Fig

  Fig. 7. Summary of Tests 1-4 showing that multiple tests are needed to assess phylogenetic support. Results from Tests 3 to 4 taken from "focal tree". Tests 1 and 2 results are the percentage of relevant trees in which the relationship was recovered (trees lacking the taxa in question were omitted from this calculation). "NA" indicates the relationship could not be tested (see Appendix H). Test 3 values are weighted averages. See Section 2.2 for explanation. Test 4 results indicate if a relationship was uniquely supported by an AU test, or if other relationships were also plausible (not unique). Black cells indicate results supporting the relationship in question, white cells are results showing lack of support, and grey cells are ambiguous.

  . 2. Focal tree topology (1geneMin tree) for major lineages in Blattodea and support. Tree was inferred in RAXML from molecular data for six loci obtained for 575 taxa. Support values are bootstrap frequencies and internode certainty values. We see strong support for the

	monophyly	of	Blattodea,	Blaberoidea,
	termites + Cryptocercidae, and termites; mod-
	erate	support	for	Lamproblattidae +
	Cryptocercidae + termites,	Tryonicidae +
	Blattidae, and Blaberidae; low support for
	Anaplectidae + sister	group,	Blattoidea,
	Corydioidea, and most relationships in
	Blaberoidea. Full tree with support values is
	available in Appendix D. See Fig. 1 caption for
	family-name usages.

  Summary of Tests 1-4 showing that multiple tests are needed to assess phylogenetic support. Results from Tests 3 to 4 taken from "focal tree". Tests 1 and 2 results are the percentage of relevant trees in which the relationship was recovered (trees lacking the taxa in question were omitted from this calculation).

	"NA" indicates the relationship could not be
	tested (see Appendix H). Test 3 values are
	weighted averages. See Section 2.2 for ex-
	planation. Test 4 results indicate if a relation-
	ship was uniquely supported by an AU test, or if
	other relationships were also plausible (not
	unique). Black cells indicate results supporting
	the relationship in question, white cells are re-
	sults showing lack of support, and grey cells are
	ambiguous.

Acknowledgements

Funding: This work was supported by the National Science Foundation [award#: 1608559] and Labex BCDiv for France Thouzé.

The authors would like to acknowledge Jessica L. Ware and others for assistance with rRNA secondary structure alignment. Also, thanks to Karen Meusemann and Christoph Mayer for advice on software and analyses. Finally, a deep thank you to the editors handling this paper, and the reviewers as well.

Author contributions: FT and FL conceived the study and laid the groundwork; FT did preliminary analyses, DE and MKK composed the final alignment and did rRNA structural alignment; DE completed the final analyses, and composed the figures; DE and FL wrote the paper; FT, and MKK also contributed to the text of the paper and figures.

Appendices. Supplementary material

Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.ympev.2018.05.007.