Combining molecular datasets with strongly heterogeneous taxon coverage enlightens the peculiar biogeographic history of stoneflies (Insecta: Plecoptera)

Extant members of the ancient insect order of stoneflies exhibit a disjunct, antitropical distribution, with one major lineage exclusively occurring in the Southern Hemisphere and the other, with few exceptions, on the Northern continents. Here, we address the biogeographic distribution and phylogenetic relationships of stoneflies using a phylogenetic workflow that combines both transcriptomic and Sanger sequence datasets with heterogeneous taxon coverage. We used a dataset comprising 2997 genes derived from the transcriptomes of 30 species and Sanger sequences of seven genes for 498 species. The backbone phylogeny was mainly inferred from the transcriptomic data, whereas the Sanger nucleotide sequence data provided high species density for divergence time estimation and diversification analyses. Our results show that the biogeographic pattern we observe today is primarily more likely shaped by long‐distance over‐land dispersal than by vicariance. We inferred that the ancestors of extant stoneflies originated in the Northern Hemisphere approximately 265 Ma and were presumably restricted to this area due to climatic and geographic boundaries. Our analyses suggest that with the break‐up of Pangaea around 200 Ma and the associated climatic and geographical changes, two groups of stoneflies, the Anarctoperlaria and the Notonemouridae, dispersed to Gondwana and subsequently went extinct on the northern continents. Both groups likely dispersed across Gondwana before its break‐up into the modern continents. At least one member of another group of ‘northern’ stoneflies, the Acroneuriinae, seems to have migrated from North America to South America around 67 Ma. We found four major net diversification rate shifts, indicating rapid radiation patterns that hampered a robust phylogenetic placement of these stonefly groups. Our study provides the first conclusive evolutionary explanation for the unique distribution pattern of stoneflies.


Introduction
Stoneflies exhibit a disjunct biogeographic distribution and their highest diversity is found in the higher latitudes. The evolutionary genesis of this unusual pattern has not been fully resolved and has puzzled researchers for decades (e.g., Illies, 1965;Bȃnȃrescu, 1991;Zwick, 2000;McCulloch et al., 2016;Cui et al., 2019;Ding et al., 2019;South et al., 2021). The order is split into two suborders: (i) Antarctoperlaria, comprising 371 species and found exclusively in southern South America and Australasia, and (ii) Arctoperlaria, comprising 3400 species and primarily distributed in the Holarctic and Oriental biogeographical regions (DeWalt et al., 2020). Three small groups within Arctoperlaria are exceptions: (i) the Notonemouridae that include 120 species and are exclusively found in Australasia, southern Africa, Madagascar and the southern parts of South America; (ii) ten genera of Acroneuriinae that occur in South America; and (iii) some Afrotropical species of the otherwise North American and Southeast Asian genus Neoperla (Perlinae) (Zwick, 2000;DeWalt et al., 2020). This anti-tropical distributional pattern led to multiple studies examining their biogeographic history (e.g., Illies, 1965;Hynes, 1988;Zwick, 2000;McLellan, 2006;Fochetti & Tierno De Figueroa, 2008;Fochetti et al., 2011;McCulloch et al., 2016;Cui et al., 2019;Ding et al., 2019). Considering the poor flight ability of stoneflies and resulting slow dispersal rates (Marden & Kramer, 1994;Thomas et al., 2000;Zwick, 2000;Schultheis et al., 2002;McCulloch et al., 2009), it has generally been assumed that long-distance dispersal played no major role in explaining the biogeographic distribution of stoneflies. Instead, vicariance (i.e., the isolation of previously connected ancestral groups as a result of plate tectonics or other geographic barriers) was considered to be the driving force behind stonefly biogeography. The extreme divergence between the northern and southern groups was explained either by the two suborders having evolved after the break-up of Gondwana ('rifting vicariance hypothesis') or by the extinction of each suborder on one hemisphere while persisting on the other ('extinction vicariance hypothesis') (summarized by Zwick, 2000). However, for some cases, such as the distribution of Notonemouridae and certain subgroups of Antarctoperlaria, molecular dating approaches retrieved dates incompatible with both hypotheses and suggested instead a trans-oceanic long-distance dispersal ('trans-oceanic dispersal hypothesis') (e.g., McCulloch et al., 2016). Based on the analyses of fossil stoneflies, Cui et al. (2019) provided another possible explanation for the biogeographical pattern displayed by Notonemouridae that includes a northern origin, the subsequent dispersal to the southern continents and, finally, the group's extinction in the northern continents ('dispersal and extinction hypothesis'). Unresolved plecopteran relationships and the unclear divergence times of the major splits within the group remain significant hurdles in explaining the order's biogeography. Previous studies recovered conflicting phylogenetic hypotheses depending on the chosen taxa, morphological or molecular data, and outgroup, and/or on the selected tree reconstruction methods (e.g., Nelson, 1984;Thomas et al., 2000;Terry, 2003;McCulloch et al., 2016;Ding et al., 2019;South et al., 2021).
In the present study, we first infer the phylogenetic relationships and divergence times among the families of stoneflies combining both transcriptomic and Sanger sequence datasets. We subsequently evaluate whether divergent levels of species richness among the families can be explained by variations in lineage-specific diversification rates. Finally, we trace the biogeographical history of the families and estimate the relative roles of dispersal and vicariance on the spatial patterns observed in extant stoneflies. In order to address all these points, the analyses require both a large-scale species-level taxon sampling and a well-resolved phylogeny among the principal lineages (Davis et al., 2013;Maddison & FitzJohn, 2015). Therefore, we combine a transcriptomic dataset, including 17 stonefly species representing 12 out of 16 recognized families, with a Sanger sequence dataset, comprising the nucleotide sequences of seven genes from 479 species representing all known extant families. For doing so, we established a workflow using supertree methods to combine 2997 single gene trees inferred from the transcriptomic sequence data with one phylogenetic tree inferred from concatenated Sanger nucleotide sequences. Using this workflow, most deeper nodes (i.e., the backbone) in the tree are inferred from the transcriptomic data, which overrules, due to the much larger number of genes, the phylogenetic information from the Sanger dataset. The Sanger sequence dataset provides phylogenetic information for establishing shallow nodes and deeper nodes for which no phylogenetic information from the transcriptomic data is available. This approach was necessary, as a backbone-constraint method to combine different kinds of datasets produced phylogenetic artefacts (Supplementary information S3).

Processing of transcriptomic sequence data
Our transcriptomic sequence dataset comprised 2997 single-copy protein-coding genes sampled in a total of 30 species (including outgroup species). The 17 ingroup species represented 12 out of 16 recognized stonefly families. All ingroup species, except Nemoura sp. (Simon et al., 2012), and 13 outgroup species were sequenced in the 1KITE project (Table S1). RNA extraction, cDNA library preparation, paired-end sequencing on the Illumina HiSeq 2000 platform and de novo assembly were performed by BGI Shenzhen, China. For all libraries 150-bp paired-end reads were generated, with approximately 2.5 Gb of raw data each. Raw reads were assembled using the soapdenovo-trans-31kmer v. 1.01  as described in Peters et al. (2017). Quality assessment and subsequent submission to NCBI Sequence Read Archive and Transcriptome Shotgun Assembly databases were conducted according to the procedure given by Peters et al. (2017). Transcriptome data published with this study are deposited in GenBank NCBI under the 1KITE umbrella project Bioproject ID 183205. For the identification of orthologous transcripts, the ortholog set described by Wipfler et al. (2019) was used. All transcriptome assemblies were searched for orthologous transcripts of protein coding genes with the software orthograph v. 0.6 (Petersen et al., 2017), using the following settings: max-blast-searches = 50; blast-max-hits = 50; extend-orf = 1; and substitute-u-with = X. All identified orthologous amino acid sequences were aligned with the L-INS-i algorithm implemented in the software mafft v. 7.402 (Katoh & Standley, 2013, 2014. Amino acid sequence alignments were checked for outlier sequences. These were then realigned and remaining outliers were subsequently deleted from the multiple sequence alignments (Misof et al., 2014). The sequences of Pediculus humanus Linneaus, Acyrthosiphon pisum (Harris), Acromyrmex echinatior (Forel) and Daphnia pulex Leydig, which were used as reference species in the ortholog set, were also removed. Putative alignment ambiguities and alignment sections indistinguishable from randomly aligned sequences were removed with a modified version of the software aliscore v. 1.2 (Misof & Misof, 2009), with a special scoring function for gappy amino acid data (option −e) and the −r = taxa 2 option in order to compare all sequence pairs in each sliding window. The alignments of all genes were subsequently concatenated using the software fasconcat v. 1.0 (Kück & Meusemann, 2010). This resulted in a concatenated supermatrix comprising the amino acid alignments of 2997 genes (1 431 848 aa positions). We prepared three additional supermatrices to evaluate the potential impact of analytical artefacts and confounding phylogenetical signal in the tree reconstructions. First, we assembled a reduced dataset, which only contained the most phylogenetically informative genes according to weighted geometry quartet mapping as implemented in mare v. 0.1.2-rc (Misof et al., 2013). We used default options with the exception that all stonefly taxa were maintained in the resulting supermatrix. This resulted in a second supermatrix that included the amino-acid alignments from 1764 different genes (873 683 aa positions). To further increase the data density with respect to the targeted stonefly relationships, we optimized the 1764 gene supermatrix by keeping only genes that were found across all stonefly species. This resulted in a third supermatrix comprising amino-acid alignments of 278 genes (86 231 aa positions). Finally, compositionally heterogeneous genes were excluded from the 1764 gene supermatrix. Among-species compositional heterogeneity was inferred using the relative composition frequency variation value for each gene (Zhong et al., 2011). Relative composition frequency variation values were calculated by the bacoca v. 1.1.05 software (Kück & Struck, 2014), and all genes with a relative composition frequency variation value >0.1 were excluded (Fernández et al., 2016;Vasilikopoulos et al., 2019). This resulted in a fourth supermatrix comprising amino acid alignments of 389 genes (154 094 aa positions). We applied pairwise sequence comparisons using Bowker's matched-pairs test of symmetry (Bowker, 1948) to estimate overall deviation of global stationary, reversibility and homogeneity (SRH) of the taxa in each supermatrix. Calculations and visualization were done with the software symtest v. 2.047 (Jermiin et al., 2008). The results of the Bowker's test showed that the larger supermatrices (1764 and 2997 genes, respectively) both strongly violate the assumptions of global stationary, reversibility and homogeneity. In contrast, the 278 and 389 genes' supermatrices exhibit significantly fewer model violations ( Fig. S4a-d).

Processing of Sanger sequence data
The Sanger nucleotide sequence dataset consisted of the nucleotide sequences of 498 species and of seven genes (Table S2). All 16 families of Plecoptera were represented in this dataset and all species, whose transcriptomes were analysed, were also included. All nucleotide sequences in the Sanger sequence dataset were retrieved from NCBI GenBank and BOLD, respectively. We included the nucleotide sequences of the nuclear and mitochondrial small and large ribosomal RNA (12S, 16S, 18S and 28S rRNA), as well as of the protein-coding genes cytochrome-c-oxidase I (COI), cytochrome-c-oxidase II (COII) and histone 3 (H3). The nucleotide sequences of the genes were aligned separately with the L-INS-i algorithm implemented in mafft v. 7.402. Ribosomal RNA sequences contain highly variable regions, which cannot be unambiguously aligned. Therefore, we excluded ambiguously aligned positions with the software aliscore v. 1.2, specifying the maximum number of pairwise comparisons as described earlier. After rRNA alignment masking and concatenation of all rRNA and protein-coding gene fragments with fasconcat, the final concatenated Sanger dataset comprised 6250 nucleotide positions.

Tree reconstruction
Transcriptomic gene trees. We performed maximum likelihood (ML) tree reconstructions for each transcriptomic amino acid gene alignment based on the sets of genes used to construct the four transcriptome supermatrices (see earlier). We used these 'transcriptomic gene trees' subsequently in the supertree analyses (Fig. 1a). Therefore, first, all four amino acid supermatrices were split into single gene partitions using amas v. 0.94 (Borowiec, 2016). Individual taxa that had no sequence information in these gene alignments were subsequently excluded using the script 'sequence_cleaner.py' in biopython v. 1.7 (Cook, 2009). All single gene alignments were analysed using the ML tree search in iq-tree v. 1.6.12 (Nguyen et al., 2015;Chernomor et al., 2016). Prior to the tree reconstruction, we selected the best scoring amino acid substitution matrix for each gene with the ModelFinder algorithm implemented in iq-tree (Kalyaanamoorthy et al., 2017), including exclusively substitution matrices appropriate for nuclear markers and the protein mixture models LG4M and LG4X (Le et al., 2012). The best model for each gene was selected by means of the corrected Akaike information criterion score (Burnham & Anderson, 2002). Additionally, we conducted tree reconstruction analyses in iq-tree using the concatenated 278/389/1764/2997 gene supermatrices (see Supplementary information S3).

Sanger sequence tree.
For the ML tree reconstruction analyses of the concatenated Sanger nucleotide sequence dataset (Fig. 1b), we also used iq-tree v. 1.6.12. The three codon  Black boxes refer to the method/dataset that recovered the presenting tree. Dotted lines represent topological differences between the supertree analyses. (b) Dataset compilation and maximum likelihood tree reconstruction results of the concatenated seven Sanger sequenced genes. (c) Supertree analyses of the combined single gene trees from the transcriptomic datasets and the single Sanger sequence tree, using two different supertree approaches. Red dots indicate nodes recovered in all supertree analyses. They are used as constraints in the Sanger dataset analyses. (d) Estimation of tree topology and divergence times in a Bayesian inference framework, with topological constraints derived from the supertree analyses from the combined transcriptome and Sanger gene tree datasets. Subsequent divergence rate and biogeographical history analyses used the resulting dated tree. positions of each protein-coding gene and each rRNA gene were defined, resulting in total in 13 distinct data blocks. To find the best-fitting partitioning scheme and models, we used the ModelFinder algorithm in iq-tree. Due to small partitions, we refrained from using the free-rate model approach in iq-tree (B.Q. Minh, personal communication) and restricted the model search to those models supported by the software package beast (Suchard et al., 2018), which was subsequently used to infer divergence times. To reduce the risk of inferring a ML tree from a local optimum in tree space, we performed 50 independent ML tree searches with random starting trees. Branch support was assessed from 1000 ultrafast bootstrap (UFBoot) replicates (Minh et al., 2013), with the 'bnni' option (Hoang et al., 2018), and an increased maximum number of iterations to stop (−nm 10 000).
Supertree analyses. Supertree analyses were first performed for the four sets of 'transcriptomic gene trees' alone ( Fig. 1a). Second, supertree analyses were performed on each of the four sets of 'transcriptomic gene trees' together with the 'Sanger sequence tree' (Fig. 1c). To compare the performance of different supertree approaches for combining trees with very heterogeneous taxon coverage, two supertree methods were applied as follows: (i) matrix representation with parsimony (Ragan, 1992), as implemented in the phytools v. 0.6-27 package (Revell, 2012) for the r v. 3.5.3 statistical software (R Development Core Team, 2019). In matrix representation with parsimony, phylogenetic information is provided as binary characters specifying the branching pattern of the input source trees. This data matrix is evaluated by parsimony and integrated into a composite (super-) tree. (ii) fastrfs (Vachaspati & Warnow, 2017), which uses dynamic programming to find a binary (super-) tree that has the minimum total Robinson-Foulds distances (Robinson & Foulds, 1981) to the original input trees. Both methods were run with default parameters. We also tested other supertree methods, such as plumist (Kupczok, 2011), clann (Creevey & McInerney, 2005) and l.u.st (Akanni et al., 2014). However, these methods were either not able to calculate supertrees on the combined datasets or produced phylogenetically unreasonable trees. The popular multispecies coalescence method as implemented in astral-iii v. 5.61 (Zhang et al., 2018) could in principle handle the combined datasets, but was not applicable to our approach. astral-iii maximizes the number of quartet trees in gene trees that are shared by the species tree. Thus, the species number raised to the power of four is the weight of each input tree in the astral tree. Consequently, the contribution of the transcriptomic gene trees with only 17 species is too low to affect the supertree analyses of the combined datasets including the Sanger tree with 498 species (Siravash Mirarab, astral developer, personal communication).

Divergence time estimations
Divergence times were estimated using a Bayesian Markov chain Monte Carlo approach implemented in the software beast v. 1.10.4 (Suchard et al., 2018) and using the concatenated Sanger dataset. Tree topology and divergence times were estimated simultaneously. However, the tree search was constrained, and we assumed all groups to be monophyletic that were consistent among all supertree reconstruction results of the combined datasets (cf. Fig. 1c,d). We further specified a relaxed-clock model with an uncorrelated log-normal clock model and a birth-death (bd) tree model prior for divergence times. Substitution model parameters were set according to the ModelFinder results from the Sanger sequence dataset, and all parameter values were unlinked across partitions. Topology and divergence time estimations in beast were run for 180 million generations and sampling was done every 20 000 generations; the first 10% of the samples of the run were discarded (burn-in phase). To graphically access the convergence and mixing of parameters and the effective sample sizes, we used tracer v. 1.7.1 (Rambaut et al., 2014). Post-burn-in samples were used to construct a maximum clade credibility tree with median node heights in treeannotator v. 1.10.4 (Suchard et al., 2018). For the calibration of the plecopteran divergence time estimations, we relied on a set of seven fossils as suggested by Cui et al. (2019). The maximum age of the root was calibrated by the age of the Rhynie chart (Evangelista et al., 2019). For all calibrated nodes within Plecoptera, we specified a uniform distribution for the divergence date, with the lower bound indicated by the minimum age of the fossil layer interval (Table S3) and the upper bound represented by the age of the stem-stonefly †Gulou carpenteri (Béthoux et al., 2011) from the Middle Carboniferous of China (319.9 Ma).

Biogeographical analyses
Biogeographical analyses were conducted in r with the biogeobears v. 1.1.2 package (Matzke, 2013). biogeobears uses three different models to estimate ancestral ranges: (i) the Dispersal Extinction Cladogenesis model (Ree & Smith, 2008); (ii) the dispersal-vicariance analysis model (Ronquist, 1997); (iii) the bayarea model (Landis et al., 2013). biogeobears is further able to model the colonization of new areas with a parameter describing a founder-event speciation (+J). However, we deliberately refrained from using the +J parameter, as its use can potentially distort the ancestral range reconstruction of ancient groups with a proposed widespread distribution (Ree & Sanmartín, 2018), as is the case in Plecoptera (Fochetti & Tierno De Figueroa, 2008;DeWalt et al., 2020). To compare the fit of all three models to the given data, the Akaike information criterion corrected for small sample sizes was used. Ancestral range reconstructions were estimated using the maximum clade credibility tree from the beast analysis. Given the limited vagility of most plecopteran species, we constrained the number of maximum areas per ancestral range to four. We selected the following regions for the biogeographical analyses: (i) Palaearctic, (ii) Afrotropical, (iii) Oriental (iv), Australasia, (v), Nearctic and (vi) Neotropical (Table S4). We also used a time-stratified model with five time slices that reflect periods bounded by major changes in Earth's tectonic and climatic conditions (Table S5) following Rolland & Condamine (2019).

Diversification rate estimations
To infer whether diversification rates in stoneflies vary across time and families, we compared the specific diversification rates across the different clades in the tree. We inferred clade-specific changes in diversification in bamm v. 2.5 (Rabosky et al., 2013;Rabosky, 2014). The prior distributions on speciation (lambda) and extinction (mu) rates were estimated with the r package bammtools v. 2.1  using the 'setBAMMprior' command (lambdaInitPrior = 9.645, lambdaShiftPrior = 0.004, muInitPrior = 9.645). Markov chain Monte Carlo runs were performed with 50 million generations, sampling every 10 000 generations. The first 10% of the samples were discarded as burn-in phase. Subsequently, the mean diversification rates and the best overall shift configuration estimated as the maximum shift credibility configuration were visualized with bammtools. To account for the assumed non-random sampling of our dataset, missing taxa were considered by specifying the percentage of species that have been sampled from each plecopteran family (Table S6). Species richness values for each family were taken from the study by .

Supertrees inferred by transcriptome gene trees
All supertree reconstruction analyses of the transcriptomic gene trees support the monophyly of Antarctoperlaria, Systellognatha and Euholognatha (Figs 1a, S1 and S2). However, the phylogenetic relationships among these groups were inconsistent, with Euholognatha either inferred as sister group of Systellognatha (= Arctoperlaria) or of Antarctoperlaria. All analyses recovered consistent phylogenetic relationships within Antarctoperlaria and Systellognatha. However, no consistent phylogenetic relationships were recovered within Euholognatha. In the latter, the relationships among families differed strongly depending on the applied tree reconstruction method and on the dataset analysed. These results were similar to the ML tree reconstruction analyses of the transcriptome supermatrices (Supplementary information S3): a monophyly of Antarctoperlaria, of Systellognatha, and of Euholognatha was highly supported by all analyses, but the relationships among these groups and relationships within Euholognatha varied depending on which of the four datasets we analysed (Fig. S3).

Trees inferred by Sanger sequence dataset
The results of the best ML tree reconstruction of the concatenated Sanger sequence dataset (Figs 1B and S5) recovered a monophyletic but weakly supported Arctoperlaria (UFBoot = 51; SH-aLRT = 21.8). Within Antarctoperlaria, Austroperlidae was paraphyletic, with Crypturoperla Illies as sister to Acruroperla Illies and the remaining Austroperlidae and Gripopterygidae. Among the families in Systellognatha, the mostly well-supported relationships corresponded to the transcriptome analyses (Fig. 1a). Styloperlidae, which is missing in the transcriptome datasets, was sister group to all other families, and Kathroperla Banks (Chloroperlidae) was sister group to Perloidea. In Euholognatha, Leuctridae (except Megaleuctra Neave) was sister group to all other Euholognatha, whereas Scopuridae was sister group to Taeniopterygidae and Capniidae, including the genus Megaleuctra (Megaleuctrinae), which was missing in the transcriptome datasets. Notonemouridae, which was also missing in the transcriptome datasets, was recovered as sister group to Nemouridae. All relationships among families in Euholognatha were only weakly supported.

Supertrees inferred by the combined datasets
In the supertree analyses of the combined transcriptome gene trees and Sanger sequence tree datasets (Fig. 1c), the matrix representation with parsimony and fastrfs analyses recovered either monophyletic Arctoperlaria as sister group of Antarctoperlaria or a paraphyletic Arctoperlaria with a sister group relationship of Antarctoperlaria and Euholognatha. Taking into account the fact that the trees based on the transcriptomic datasets (with 278/389/1764/2997 genes) were used as guiding information for the supertree analyses of the combined datasets (279/390/1765/2998), all analyses were consistent: the principal relationships among the plecopteran families recovered with the transcriptomic datasets were also recovered in the subsequent supertree analyses of the combined datasets (Fig. 1a, c, Supplementary information S1). In the matrix representation with parsimony analyses, the position of Systellognatha and Euholognatha and the position of Taeniopterygidae + Scopuridae varied (Figs 1a and S1). The latter was either recovered as sister group of Capniidae + Megaleuctrinae (279 and 1765 gene tree datasets), as sister group of Notonemouridae + Leuctridae (390 gene tree dataset), or as sister group of all remaining families, except Nemouridae (2998 gene tree datasets). In the fastrfs analyses (Figs 1a and S2), Nemouridae was either sister group of Capniidae + Megaleuctrinae (279 and 390 gene tree dataset) or sister group of all other Euholognatha (1764 and 2998 gene tree datasets), whereas Notonemouridae + Leuctridae were either recovered as sister to Taeniopterygidae + Scopuridae (279 and 390 gene tree datasets) or as sister to Capniidae + Megaleuctrinae (1764 and 2998 gene tree datasets). The relationships within the families were in agreement with those of the Sanger sequence tree (Figs 1b and S5). However, due to the inconsistencies among the results of the different datasets, we were unable to use a single tree as input for the dating analyses and therefore decided to constrain the tree search within the divergence time estimation analyses only with those nodes recovered in all previous tree reconstruction analyses.

Bayesian tree reconstruction and divergence time estimation results
The phylogenetic tree reconstruction and divergence time estimation based on the constrained Bayesian analysis in beast (Figs 2 and S6) of the Sanger dataset places the origin of Coloured branches indicate relative speciation rates along each branch as inferred by the evolutionary rate analyses in bamm (increasing from blue to red). Red stars indicate the positions of a diversification rate shift in the maximum shift credibility configuration of the bamm analyses. Green dots indicate the topological constraints as inferred from the supertree analyses. Black/white pie charts represent node support in the beast analyses. Eu., Eusthenioidea; Neo., Neogene; Pal., Paleogene; Perm., Permian; Ptero., Pteronarcyoidea.  Table 1, and the beast treeannotator output is presented in Supplementary information S2.

Divergence rate estimations
The bamm phylorate plot of speciation rate and the maximum shift credibility score (Fig. 2) display considerable but different patterns of rate heterogeneity within Euholognatha and Systellognatha. In Systellognatha, a single rate shift leading to an increased speciation was found at the base of the superfamily Perloidea (Chloroperlidae, Perlidae, Perlodidae). In contrast, in Euholognatha, higher speciation rates were inferred within the families Capniidae, Leuctridae and Nemouridae. The rate-through-time plots (Fig. S7), which display rate changes over time in the different groups, also indicate increased net diversification in Perloidea, Capniidae, Leuctridae, Nemouridae and Perloidea.

Biogeographical analyses
The results of the best biogeographical model are shown in Fig. 3. The analyses conducted with the bayarea model are significantly better supported compared to the analyses with the two other models (Table 2). An origin in a combined ancestral area including the Nearctic and Palaearctic regions (Holarctic) was recovered at the root. Within Arctoperlaria, both Systellognatha and Euholognatha radiated in the Holarctic region throughout the Jurassic to the Middle Cretaceous. Most families radiated into the Nearctic later, while the break-up of Laurasia with the opening of the North Atlantic Ocean in the Late Cretaceous led to a radiation in the Palaearctic. In Euholognatha, the origin of the exclusive southern Notonemouridae remains ambiguous. A first splitting group colonized South Africa, but all other groups dispersed in the Neotropical and/or Australasian regions. Antarctoperlaria simultaneously dispersed to the Southern Hemisphere in the Middle Jurassic, but their origin was estimated in the Neotropical and Australasian region, with a later major radiation in the latter.

Stonefly phylogeny and diversification rates
In recent large-scale phylogenetic studies, genomic data have been frequently combined with Sanger-sequenced standard molecular marker genes (e.g. Leaché et al., 2014;Fernández et al., 2018). For this purpose, backbone-constraint methods have been developed or suggested (Stamatakis, 2006;Nguyen et al., 2015;Kjer et al., 2016), which add taxa to a given tree by respecting the restrictions of the backbone-constraint tree. This  method of combining different datasets proved to be problematic for the reconstruction of stonefly relationships (Supplementary information S3), as we retrieved the polyphyly of many traditional families, such as Capniidae, Leuctridae, Perlidae, Perlodidae, Peltoperlidae and Taeniopterygidae, which contradict previous molecular (Terry, 2003;McCulloch et al., 2016;Ding et al., 2019;South et al., 2021) and morphological studies (Nelson, 1984;Zwick, 2000). We thus consider them as methodological artefacts, most likely caused by the extreme heterogeneity of data size and included species between the two employed data types (30 transcriptomes and 278/389/1764/2997 genes) versus 498 (Sanger species and 7 genes). In the supertree approach that we applied to this study, the strong weight of the transcriptomic dataset stems from many thousand gene trees that 'guide' the final supertree analyses, and which stand against the single tree inferred in a concatenated analysis of the Sanger sequence data. We thus recommend our workflow ( Fig. 1) for future application on datasets where similar problems are encountered with a constraint-based tree inference approach. According to our results, the stoneflies separated into the three clades Antarctoperlaria, Systellognatha and Euholognatha (Figs 1a-c, S1-S3). We recovered stable results within Systellognatha that mainly corroborate the results by Terry (2003) and South et al. (2021) by retrieving monophyletic Perloidea nested within paraphyletic Pteronarcyoidea (Fig. 2). For Perloidea, we recovered the pattern suggested by Illies (1965), with Chloroperlidae (except Kathroperla) forming the sister group of Perlodidae, contradicting all previous studies (Zwick, 1973(Zwick, , 2000Nelson, 1984;Thomas et al., 2000;McCulloch et al., 2009;Ding et al., 2019) except that by South et al. (2021). Kathroperla as sister group to all remaining Perloidea is congruent with the findings in the study by Terry (2003). South et al. (2021) retrieved an identical pattern in some of their tree reconstructions, but most of their datasets and analyses recovered Kathroperla as sister to Chloroperlidae + Perlodidae. Although the relationships within Systellognatha were consistent among all our analyses, within Euholognatha only the sister groups Megaleuctra + Capniidae gained consistent support. Illies (1967) also proposed a closer relationship to Capniidae based on characters of the male genitalia. However, this position is challenged by more recent morphological studies on wing venation and larval characters, which frequently recovered the subfamily Megaleuctrinae (= Megaleuctra) as sister to all remaining Leuctridae (Zwick, 1973(Zwick, , 2000(Zwick, , 2006Béthoux, 2005;Béthoux et al., 2015) and support its assignment to Leuctridae by Frison (1935Frison ( , 1942. With comprehensive transcriptome datasets, South et al. (2021) also consistently recovered Megaleuctra as sister group to the remaining Leuctridae.
The inconsistencies among the different supertree analyses and datasets indicate that the phylogenetic information of the combined dataset is insufficient to reliably resolve the relationships between Antarctoperlaria, Euholognatha and Systellognatha and to reliably resolve relationships within Euholognatha. For the latter, our transcriptomic dataset was also missing representatives of Notonemouridae. Thus, the position of this family is exclusively based on Sanger sequences, limiting the amount of data used to resolve relationships among Notonemouridae and the other families in Euholognatha. However, despite the more comprehensive transcriptome dataset, South et al. (2021) also recovered ambiguous relationships within Euholognatha. In our divergence time estimation analysis, the times between the early branching events are generally longer in Systellognatha than in Euholognatha, and we also detected three net diversification rate escalations in the euholognathan families: Capniidae, Leuctridae and Nemouridae (stars in Figs 2 and S7). The short internodes in Euholognatha and the accelerated diversification rates in different lineages imply that rapid radiation events, which take place over relatively short time periods, restrict the accumulation of phylogenetic signal among the targeted clades (Whitfield & Lockhart, 2007;Whitfield & Kjer, 2008). Patterns like these are assumed to hamper phylogenetic analyses and might be the cause for the consistent phylogenetic instability of euholognathan and arctoperlarian relationship estimates. These problems remain in larger phylogenomic datasets: the implementation of ∼100 arctoperlarian transcriptomes by South et al. (2021) recovered a conspicuously higher support for the phylogenetic pattern in Systellognatha, whereas uncertainties among euholognathan family relationships persisted. These problems will therefore not vanish by simply adding more data (e.g., Kubatko & Degnan, 2007). Instead, it will be necessary to carefully check for potential sources of tree reconstruction errors (Scornavacca & Galtier, 2017), such as undetected paralogous genes, introgression and incomplete lineage sorting.

Stonefly biogeography
The colonization of the Northern Hemisphere. Extant stoneflies evolved from stem group representatives that inhabited the northern part of Pangea (Holarctic), considerably before its break-up (Fig. 4a). Gondwana plecopteran fossils from the Permian (e.g., Van Dijk & Geertsema, 2004) thus represent stem lineage stoneflies as suggested previously (Zwick, 2000;Grimaldi & Engel, 2005). The global distribution of the stem group is plausible when considering that stoneflies separated around 350-400 Ma from their sister groups that are also all distributed globally (Tong et al., 2015;Evangelista et al., 2019). The emergence times of Antarctoperlaria, Euholognatha and Systellognatha between 265 Ma (95% HPD: 294-236 Ma) and 256 Ma (95% HPD: 283-228 Ma) are much older than those of previous divergence time estimates that proposed ∼181 Ma (Ding et al., 2019) and∼121 Ma (McCulloch et al., 2016), respectively, but agree with estimates based on the fossil record (Cui et al., 2019). The three major lineages of stoneflies originated in the Holarctic and we thus clearly refute the proposed vicariance theories as explanation for the current stonefly distribution that either implies a complete Pangean distribution and the subsequent extinction in either hemisphere ('extinction vicariance hypotheses') or the evolution of specific subgroups on either Gondwana or Laurasia from a Pangean stem lineage after the break-up of Gondwana ('rifting vicariance hypotheses'). The northern isolation of the early extant stoneflies is concordant with the proposed Earth's climate and geography at the end of the Permian. Warm to temperate humid climates, which stoneflies require, were limited to mid and higher latitudes (Cui & Kump, 2014;Miller & Baranyi, 2019). Low to mid latitudes were characterized by semi-arid to arid areas (Chaboureau et al., 2014), forming a belt-like barrier for the southern distribution of the fresh water-dependent stoneflies (a in Fig. 4a). Additionally, the collision of Gondwana and Laurasia created a large range of fold mountains, the Variscan mountain range (Kroner & Romer, 2013) that formed another natural boundary (v in Fig. 4a). These factors created an effective barrier that impeded the first representatives of modern plecopteran from spreading South, despite the fact that the climatic conditions on the southern continents allowed the existence of plecopterans stem group representatives (see earlier). Within the families of Systellognatha and Euholognatha (except Notonemouridae), our results show a typical vicariance pattern over their ancestral range of the Holarctic. When the Palaearctic started to separate from the Nearctic in the Palaeogene, most extant genera already existed and are thus found nowadays in North America and Eurasia (Fig. 3).

The first conquest of the Southern Hemisphere.
At the beginning of the Jurassic, the disintegration of the Pangaean supercontinent was correlated with deep changes in the Earth's climate, leading to a more humid and warmer climate with increased sea level and precipitation (Chaboureau et al., 2014). The decreasing arid zones resulted in an increased connectivity between humid tropical and warm temperate zones. Additionally, the drifting apart of Gondwana and Laurasia led to the flattening of the Variscan mountain range. Two groups of stoneflies, the Notonemouridae and Antarctoperlaria (Fig. 3), dispersed South and started to diversify on Gondwana around 180 Ma (Antarctoperlaria: 185 Ma, 95% HPD 228-146 Ma; Notonemouridae: 181.6 Ma, 95% HPD 208-161 Ma). Both groups colonized the South independently, as the Notonemouridae are nested deep within Euholognatha (Fig. 2). The timing of these Laurasia-Gondwana dispersal events must have occurred between the above-described climatic improvements and mountain flattening and the beginning of their Gondwana radiations. Palaeontological data for the Notonemouridae (Cui et al., 2019) suggested a date around 200 Ma, which is in agreement with our results (Fig. 4c). We thus conclude that the distinct geographic pattern in the extant stoneflies with their peculiar North-South separation is explained by long-distance dispersal over land, while the continents were still connected, and we refute a trans-oceanic migration proposed by McCulloch et al. (2016). Instead, we can show that the pattern proposed by Cui et al. (2019) for the Notonemouridae also apply to the Antarctoperlaria. Both groups colonized Gondwana rapidly and appeared in Australasia -the most distant region from the northern continents by that time -around 165 Ma (Fig. 4d) and thus before Gondwana started to break apart. This challenges the theory that stoneflies are poor dispersers due to their restricted flight ability (Brundin, 1967(Brundin, , 1972Zwick, 2000;Schultheis et al., 2002;Fochetti & Tierno De Figueroa, 2008;McCulloch et al., 2009). The rapid and long-distance dispersal ability of stoneflies has been shown for some species that, for example, migrated along drainage courses after the last glaciation period (Pessino et al., 2014) or that might have used connected streams and rivers to disperse over North America (Kauwe et al., 2004). We cannot trace the exact migration route of the two groups over Gondwana due to their subsequent extinction in Antarctica, India and for Antarctoperlaria also in southern Africa, which was caused by cooling in the former and increased temperature and aridity in the Cretaceous in the latter (Fochetti & Tierno De Figueroa, 2008;Cui et al., 2019). But within Notonemouridae, the Afrotropical species form the sister group of the remaining taxa, which might indicate that this was the origin of their southern dispersal. This is also supported by the palaeoclimate, as the eastern part of the Laurasia-Gondwana border, which included Africa, rapidly became humid and thus less inhabitable by stoneflies than the western part between North America and South America (Fig. 4c,d). It appears likely that these stoneflies relied on a system of connected flowing waters to achieve this dispersal as was reported for Pteronarcys californica Newport in western North America (Kauwe et al., 2004). The hypothesis of long-distance dispersal as driving force of stonefly biogeography is further supported by the distribution of two pairs of closely related antarctoperlarian species [Stenoperla helsoni McLellan -Neuroperla schedingi (Navás) and Austroheptura Illies -Klapopteryx armillata Navás]. Both pairs include one species from Australasia and one species from South America, both of which diverged around 65-75 Ma from their respective sister group (Fig. 2). At that time period, the only possible exchange route would be via Antarctica (Fig. 4g), which connected these two continents until the opening of the Drake Passage around 41 Ma (Scher & Martin, 2006). These cases provide additional evidence for two independent cases of long-term dispersal between South America and Australia.
The second conquest of the Southern Hemisphere. The driving force behind the developments that allowed the Notonemouridae and Antarctoperlidae to disperse to the South -the separation of Gondwana and Laurasia -eventually resulted in the opening of the Tethys Ocean around 160 Ma (s in Fig. 4d). This oceanic boundary then prohibited any exchange between the Laurasia and Gondwana groups until the end of the Mesozoic (Fig. 4d-f). This long-lasting North-South separation of stoneflies was broken approximately 67 Ma (95% HPD: 88-48 Ma), when one group of Perlidae, the Acroneuriinae colonized South America (Fig. 3). The timing and migration pattern of this event was strongly debated, and possible explanations include vicariance with a slow acclimatization of originally cool-adapted animals to tropical conditions, a Miocene dispersal from North America, or independent Cretaceous dispersals of the subgroups Acroneuriini and Anacroneuriini from Asia or Europe via Africa (see Zwick, 2000 and the literature therein). Neotropical Acroneuriinae form a monophyletic lineage that colonized South America from North America via long-distance dispersal (Fig. 4g). As there was no land bridge between the continents at that time, they had to disperse over the ocean. However, a chain of meso-American islands existed that might have been used as stepping stones for the dispersal (Fig. 4g). This potential island-hopping could have decreased the distance to either actively fly or passively disperse via air current to reach land and fresh water. DeWalt & South (2015) showed that stoneflies can bridge distances over water of at least 20 km to reach islands. Small species with long wings are thereby better suited for long-distance dispersal as they can be lifted by air currents (Malmqvist, 2000;DeWalt & South, 2015). For a group of poor fliers, such as stoneflies (Zwick, 2000;Fochetti & Tierno De Figueroa, 2008;McCulloch et al., 2016), this might be an important aspect that is further underlined by the fact that this is the only case where we found stoneflies to cross an ocean. Some species of the genus Neoperla (Perlidae) also left their ancestral northern distributional range and colonized a southern biogeographic region, the Afrotropics (Fig. 4h). Our taxon sampling included three species of Neoperla from North America, but unfortunately no species from Africa. We thus cannot address this southern dispersal in detail, which is assumed to have taken place in the Miocene or Pliocene (Zwick, 2000). However, the genus Neoperla evolved around 92 Ma (95% HPD: 113-71 Ma) and the youngest included species is 40 Ma (95% HPD: 61-21 Ma), which might indicate a much earlier colonization of Africa (Fig. 3).

Conclusion
Our analyses provide the first conclusive picture of the evolution and biogeography of the major stonefly lineages, indicating that vicariance did not play a major role concerning the anti-tropical biogeographic pattern we observe today except in their ancestral northern range. Instead, groups inhabited other biogeographic regions via long-distance dispersal over land as in Antarctoperlaria and Notonemouridae. Our results thus refute the generally accepted idea that vicariance is the driving force behind ancient stonefly biogeography as proposed by original studies (e.g., Illies, 1965;Zwick, 2000) or textbooks (e.g., Resh & Cardé, 2009). These migration patterns are in accordance with the limited flight abilities of stoneflies that impedes any long-distance dispersal by single individuals. Instead, the group apparently dispersed over many generations requiring constant access to fresh water for the aquatic immature stage. This dependency on fresh water makes oceans and arid regions strong natural barriers to stonefly dispersal. However, our results also show that stoneflies can bridge long distances considerably quickly, which add new complexity to the generally accepted hypothesis of stoneflies as poor dispersers. Another common pattern in stonefly zoogeography is extinction. Various groups disappeared from specific geographic areas such as the Notonemouridae and Antarctoperlidae from the northern continents, Antarctica and India, or various groups of Systellognatha that went extinct in parts of the Holarctic. From a methodological standpoint, we show that supertree analyses are a reliable and convenient tool for combining Sanger with transcriptomic data. We expect this approach to become more important as data heterogeneity in terms of included species (which will grow faster for Sanger data) and data size (which will grow faster for genomic datasets) will increase. Using the example of stonefly phylogenies, we presented a reliable tree reconstruction approach that combines the two major trends in molecular phylogenetics, the increased sequencing of transcriptomes and genomes and the Sanger sequencing of large datasets including entire museum collections.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.            Data S1: Contains all trees generated by supertree analyses of the combined datasets.
Data S2: Contains the results of the BEAST Bayesian Inference analyses, showing the relationships and divergence times of stoneflies.
Data S3: Contains the methods and results of the constrained IQ-TREE analyses of the concatenated Sanger sequence genes, with the best trees from the ML analyses of the four transcriptomic datasets.