Molecular phylogeny of the megadiverse insect infraorder Bibionomorpha sensu lato (Diptera)

The phylogeny of the insect infraorder Bibionomorpha (Diptera) is reconstructed based on the combined analysis of three nuclear (18S, 28S, CAD) and three mitochondrial (12S, 16S, COI) gene markers. All the analyses strongly support the monophyly of Bibionomorpha in both the narrow (sensu stricto) and the broader (sensu lato) concepts. The major lineages of Bibionomorpha sensu lato (Sciaroidea, Bibionoidea, Anisopodoidea, and Scatopsoidea) and most of the included families are supported as monophyletic groups. Axymyiidae was not found to be part of Bibionomorpha nor was it found to be its sister group. Bibionidae was paraphyletic with respect to Hesperinidae and Keroplatidae was paraphyletic with respect to Lygistorrhinidae. The included Sciaroidea incertae sedis (except Ohakunea Edwards) were found to belong to one clade, but the relationships within this group and its position within Sciaroidea require further study.

In terms of biodiversity, Bibionomorpha is a megadiverse group due to the inclusion of the fungus gnats (Sciaroidea, comprising the very large families Mycetophilidae and Sciaridae) and gall midges (family Cecidomyiidae), the latter presumably even being the most diverse and species-rich family of Diptera (cf. Hebert et al., 2016). The number of extant species of Bibionomorpha sensu lato currently described has been estimated at 15,000 (Pape, Bickel & Meier, 2009), although an inestimable number of species in this group still remain uncollected and undescribed. Moreover, fungus gnats and gall midges are notoriously abundant in trap catches (e.g., Malaise traps) from terrestrial habitats, especially mesic forests. Various subgroups of Bibionomorpha are the most speciose among fossil Diptera, being well represented in the fossil record since the Mesozoic and impressively documented from different ambers (Evenhuis, 1994;Blagoderov & Grimaldi, 2004;Grimaldi, Engel & Nascimbene, 2002;Hoffeins & Hoffeins, 2014).
The larval diets of Bibionomorpha are diverse, including detritophagy, saprophagy, predation, mycophagy and phytophagy. Mycophagy has been considered to be ancestral in Sciaroidea, and predation ancestral in Keroplatidae (Matile, 1997). However, these conclusions were based on relatively little empirical evidence and the biology of most Bibionomorpha, even on a generic level, remains understudied. As a notable exception, the biology of many phytophagous Cecidomyiidae has been studied in great detail (e.g., Gagné & Moser, 2013). As for adults, fungus gnats are certainly the most conspicuous bibionomorphs, since they are both abundant (usually aggregating in large numbers at the trunks of fallen, rotten trees, along stream banks, and at similar moist, shady places) and big enough to be noticed with the naked eye. Species of Bibionidae occurring in enormous numbers during spring are widely known, even among general naturalists, as March flies, or lovebugs.
There are a number of unresolved questions regarding Bibionomorpha phylogeny that we have aimed to address in this study: Firstly, the delimitation of the infraorder (sensu stricto versus sensu lato) remains unclear, especially regarding the question whether families, such as Anisopodidae sensu lato, Scatopsidae, Canthyloscelidae, Axymyiidae, and Perissommatidae, belong to the Bibionomorpha or not.
Secondly, the phylogenetic position of non-sciaroid families, such as Bibionidae, Hesperinidae, and Pachyneuridae is still unresolved, in part due to the fact that hesperinid and pachyneurid specimens are rarely collected and poorly represented in collections. Representatives of these three families were included in the morphological analyses by both Fitzgerald (2004) and Amorim & Rindal (2007), which resulted in similar phylogenetic hypotheses, but these hypotheses still need to be tested, especially as the latter study was criticized for fundamental methodological shortcomings (Jaschhof, 2011).
Thirdly, the delimitation of the Sciaroidea is still a controversial issue, especially regarding the inclusion of the Cecidomyiidae and/or Ditomyiidae. Several studies, both molecular (Bertone, Courtney & Wiegmann, 2008) and morphological (Fitzgerald, 2004), have been unable to provide support for the monophyly of Sciaroidea (including Ditomyiidae) and most of the studies have not recognized Cecidomyiidae as a part of Sciaroidea.
Fourthly, and perhaps the most puzzling issue of all, is the phylogenetic position and assignment to family of the Sciaroidea incertae sedis, which is closely related to the question of how the family Rangomaramidae should be appropriately defined (cf. Jaschhof & Didham, 2002;Amorim & Rindal, 2007;Jaschhof, 2011). The authors of the present paper do not follow the concept of the Rangomaramidae as proposed by Amorim & Rindal (2007), a practice in common with the most recent Diptera Manuals (Brown et al., 2009;AH Kirk-Spriggs & BJ Sinclair, 2016, unpublished data) as well as papers specifically related to Sciaroidea incertae sedis (e.g., Hippa & Ševčík, 2014). Nevertheless, Amorim & Rindal's (2007) results and proposals are discussed here when appropriate.
The present paper aims to provide answers, or partial answers, to the four questions outlined above using nuclear and mitochondrial gene markers and taxon sampling designed to address these Bibionomorpha-specific issues.

Taxon sampling
A total of 94 terminal taxa are included in the dataset (Table 1 and Table S1). The ingroup contains 60 species from all the families traditionally placed in Bibionomorpha Table 1 List of specimens used for DNA extraction, with GenBank accession numbers. More information about the specimens is listed in the Table S1.
Nannochorista sp. n/a n/a n/a n/a n/a KC177108
As outgroup taxa, we analysed 20 species from 13 families of non-bibionomorph lower Diptera, 13 species from 12 families of Brachycera, and one species of Mecoptera. Samples containing these species were collected throughout the world, usually by means of Malaise traps, in the years 2006-2015 (Table S1). Most of the specimens were collected in localities where no field study permission is required. No taxon included in this study is specially protected by law. All the specimens used for DNA extraction were identified by the authors (JŠ, SF, MJ), according to the most recent taxonomic literature. The voucher specimens are deposited in the reference collection of the Ševčík Lab, University of Ostrava, Czech Republic (JSL-UOC), with several specimens also deposited in the collection of the Silesian Museum, Opava, Czech Republic (SMOC); see Table S1.
Some sequences, mostly for outgroup taxa, were taken from the GenBank database. We were not able to obtain museum specimens of Rangomaramidae or Perissommatidae for our analyses so any efforts in the future to illuminate the phylogenetic position of these enigmatic groups by molecular approaches will depend on freshly collected material.
Sequencing was carried out with BigDye Terminator ver.3.1 (Applied Biosystems, Foster City, USA) on an ABI 3100 genetic analysis sequencer (Perkin Elmer Applied Biosystems, Norwalk, CT, USA), or PCR products were sequenced by Macrogen Europe (Netherlands). All sequences were assembled, manually inspected, and primers trimmed in SEQUENCHER 5.0 (Gene Codes Corporation, Ann Arbor, USA). Sequences used to build an alignment for phylogenetic analysis were deposited in the GenBank database, accession numbers are listed in Table 1.
The identity of all the sequences was confirmed by BLAST similarity searches against NCBI database and double-checked in single gene trees to eliminate possible contaminations or other misleading results. Suspicious-looking sequences, where possible, were sequenced repeatedly from different specimens of the same species as a positive control (e.g., in Asiorrhina parasiatica, Diadocidia globosa, Insulatricha hippai, Lygistorrhina spp., Ohakunea bicolor, Protaxymyia thuja, Sciarosoma nigriclava).

Sequence alignment and phylogenetic analyses
The ribosomal genes 12S, 16S, 18S and 28S were aligned using MAFFT version 7 (Katoh & Standley, 2013) on the MAFFT server (http://mafft.cbrc.jp/alignment/server/). Method L-INS-I, recommended for < 200 sequences with one conserved domain and long gaps, was automatically selected by the program according to the alignment sizes. Resulting alignments were visually inspected and manually refined when necessary. The lengths of individual alignments were: 12S = 767 bp, 16 = 399 bp, 18S = 2,554 bp, 28S = 1,399 bp. We also tested some other alignment algorithms (Q-INS-I MAFFT-which considers secondary structure of the molecules; or T-Coffee), but none of them yielded a tree with better resolution.
All unreliably aligned regions were removed in the program GBLOCKS 0.91b (Castresana, 2000) on the Gblocks server (http://molevol.cmima.csic.es/castresana/ Gblocks_server.html). Conditions have been set as follows: allowed smaller blocks, allowed gap positions within the final blocks and allowed less strict flanking positions. Sequences of the protein-coding genes COI and CAD were checked based on amino-acid translations and provided indel-free nucleotide alignments. The single-gene alignments were concatenated and the final concatenated alignment was deposited in the Figshare database (https://dx.doi.org/10.6084/m9.figshare.c.3458082.v1).
We made a comprehensive alignment of 94 terminal taxa comprising almost all major groups of the infraorder Bibionomorpha, as well as a number of outgroup taxa, including several representatives of Brachycera. We selected representatives of diverse lineages and various trophic strategies in each family of the ingroup to make up a balanced dataset. The final data matrix consisted of 5,018 characters: 12S-306 bp, 16S-326 bp, 18S-1,851 bp, 28S-452 bp, CAD-1425, COI-658 bp. Saturation analysis was performed for all six genes using Xia's test implemented in DAMBE (Xia & Xie, 2001). In the case of the protein coding genes COI and CAD, saturation of third codon positions was tested separately.
Trees were rooted by the representative of the order Mecoptera, which is considered sister group to Diptera (e.g., Misof et al., 2014).
The dataset was analysed using Bayesian inference (BI) and maximum likelihood (ML) methods. The node support values are given in the form: posterior probability (PP) / bootstrap value (BV). We do not present here any results from the maximum parsimony analysis, because it provided almost no significant support at the higher taxonomic levels (interfamilial relationships), in concordance with theoretical arguments indicating the limited performance of the MP method in molecular phylogenetics (e.g., Gadagkar & Kumar, 2005;Spencer, Susko & Roger, 2005).
The partitioned Bayesian inference of 50 million generations on the concatenated dataset was implemented in MrBayes version 3.2.2 (Huelsenbeck & Ronquist, 2001) and carried out on the CIPRES computer cluster (Cyberinfrastructure for Phylogenetic Research; San Diego Supercomputing Center, Miller, Pfeiffer & Schwartz, 2010). Burn-in was set as 30%.
The ML analyses were conducted on the CIPRES computer cluster using RAxML-HPC BlackBox 7.6.3 (Stamatakis, 2006) employing automatic bootstrapping on single gene alignments to check potential conflicting topologies and other artifacts. Subsequently, the same algorithm was applied on the concatenated partitioned dataset.
Phylogenetic trees were visualized using Interactive Tree Of Life (iTOL) (Letunic & Bork, 2011). The trees presented are Bayesian or ML topologies with node support values from both the analyses.

Comparison of the Bayesian and maximum likelihood analyses
The results based on Bayesian inference (BI) and maximum likelihood (ML) analyses of the dataset are summarized in Figs. 1-3, Figs. S1 and S2. For BI, standard deviation of split frequencies was in all cases < 0.01. The log likelihood value for the best tree of the dataset was −128022.18. Both the model-based methods (BI, ML) yielded mostly congruent nodes, with several exceptions mostly at the family level (Figs. 1-3). Well-supported relationships were consistent across all the trees obtained.
Saturation analyses revealed that the phylogenetic markers used in this study show a rather low level of saturation, with the exception of the third positions of the protein coding genes COI and CAD. However, ML phylogenetic analysis based on the dataset with excluded third codon positions provided a tree with very similar topology (cf. Figs. 1 and 2) in comparison with the original tree, with small differences only within several family-level clades and also slight changes in node support values. These minor differences are discussed below when relevant.

Delimitation, monophyly and affiliation of Bibionomorpha
The monophyly of Bibionomorpha sensu stricto (containing Sciaroidea, Bibionidae, Hesperinidae, and Pachyneuridae) was established with high support (PP = 1.00/BV = 100) in both the model-based analyses (see Figs. 1-3, Figs. S1 and S2). The monophyly of Bibionomorpha sensu lato (containing Bibionomorpha s. str. plus Anisopodidae sensu lato, Canthyloscelidae, and Scatopsidae, 1.00/99) was also established. These findings are consistent with results from the previous molecular and morphological studies, which usually found Bibionomorpha (in various concepts) to be monophyletic (e.g., Hennig, 1954;Oosterbroek & Courtney, 1995;Amorim & Yeates, 2006;Bertone, Courtney & Wiegmann, 2008;Wiegmann et al., 2011). Amorim & Yeates (2006) even suggested raising the Bibionomorpha to the rank of suborder. The present study shows Bibionomorpha sensu lato to be sister group to Brachycera (1.00/99), which is congruent with results from the molecular studies by both Bertone, Courtney & Wiegmann (2008) and Wiegmann et al. (2011). Hennig (1973, who came to the same result through morphological study, referred to two synapomorphies supporting this relationship: the conspicuous enlargement of the second latero-tergite and the undivided postphragma of the thorax. Fitzgerald (2004) considered the structure of the parameres (presence of a dorsal sclerite and ventrolateral apodemes) as synapormorphies of this taxon and Michelsen (1996) even introduced a new taxon, Neodiptera, for Bibionomorpha sensu lato + Brachycera based on several cervical and prothoracic skeletomuscular characters as synapomorphies.
On the other hand, the Bayesian tree by Beckenbach (2012), based on the complete mitochondrial genome of 24 species of Diptera, revealed Anisopodidae + Tanyderidae to be the sister group to Brachycera, but with relatively low support. The morphological study by Oosterbroek & Courtney (1995) considered only Anisopodidae (rather isolated in their tree from the rest of Bibionomorpha) as the sister group to Brachycera, with Anisopodidae + Brachycera being the sister group to the clade containing Tipulidae and Trichoceridae. Their Bibionomorpha included the family Axymyiidae and formed a sister clade to the whole ''higher Nematocera and Brachycera'' (Perissommatidae, Scatopsidae, Synneuridae, Psychodidae, Tipulidae, Trichoceridae, Anisopodidae, and Brachycera).
It is beyond the scope of this study to discuss the details of all these morphological analyses, but we think that a general failure of the past is underappreciation of the relevance of Catotrichinae, which according to the prevailing opinion (cf. Gagné & Jaschhof, 2014) are sister group to all the other Cecidomyiidae. The subfamily Catotrichinae contains only a few species, which are seldom treated in publications (most recently by Jaschhof & Jaschhof (2008), so too late to be taken in consideration by the studies referred to above).
Even more important, hardly any dipterist has seen specimens of Catotrichinae so as to appreciate their remarkable resemblance with fungus gnats rather than with other gall midges. As another problem, knowledge of Cecidomyiidae morphology is evidently often too vague to interpret characters and character states correctly (cf. Hippa & Vilkamaa, 2005;Amorim & Rindal, 2007).
The molecular analysis by Bertone, Courtney & Wiegmann (2008) recovered a paraphyletic Sciaroidea, with Ditomyiidae excluded, and with moderate support 0.92/76 for a clade consisting of Ditomyiidae + (Bibionidae + Pachyneuridae). The incongruence between the results of Bertone, Courtney & Wiegmann (2008) and our results might possibly be explained by the more extensive taxon sampling (almost three times as many species of Bibionomorpha sensu lato included) in the present study.

Monophyly and interrelationships of bibionomorph families
Although the sampling was not designed to thoroughly answer questions on the monophyly of particular families, most of the previously recognized families of Bibionomorpha were monophyletic with high support in all the analyses performed (Figs. 1-3, Figs. S1 and S2). The paraphyly of Keroplatidae is discussed below. The paraphyly of Mycetophilidae (with respect to Ohakunea Edwards, a genus usually regarded as unplaced to family) as suggested in Fig. 1 is surprising and is also commented on below.
With respect to the position of Axymyiidae, our analysis produced a surprising result. While some previous studies placed Axymyiidae as the sister group to Bibionomorpha sensu stricto (e.g., Oosterbroek & Courtney, 1995;Wiegmann et al., 2011) or in a clade together with Pachyneuridae and Perissommatidae (e.g., Hennig, 1973;Amorim, 1993), our study grouped Axymyiidae with Culicomorpha, though with only moderate support (1.00/56). To exclude the possible long-branch attraction artifact in this case, we tested the ML topology after removal of the long-branch forming Culicomorpha clade. The resulting topology confirmed the position of Axymyiidae within the outgroup taxa. A similar topology was introduced by Bertone, Courtney & Wiegmann (2008), with Culicomorpha being the sister group to Axymyiidae + Nymphomyiidae. The morphological study by Fitzgerald (2004) also found Axymyiidae neither within nor sister group to Bibionomorpha. While several studies seem to indicate that Axymyiidae may not belong within, or as sister group to, Bibionomorpha, the phylogenetic position of Axymyiidae clearly requires further study.
As regards the families of Sciaroidea, BV supports from our ML analysis are generally rather low to infer interfamilial relationships, but PP supports for our Bayesian analysis provide better resolution. Concerning the relationships between families within Sciaroidea, the only well supported node in both the model-based analyses (1.00/100) is Keroplatidae, excluding Arachnocampa Edwards, 1924, and containing Keroplatinae, Macrocerinae and Lygistorrhinidae; see Figs. 1 and 2 and 'Paraphyly of Keroplatidae'.
Ditomyiidae was shown to branch basally from the rest of Sciaroidea, with high support in both the analyses (1.00/96). This result is supported by most previous studies based on morphology (Zaitzev, 1983;Zaitzev, 1984;Wood & Borkent, 1989;Matile, 1990;Matile, 1997;Hippa & Vilkamaa, 2006), pointing to several primitive characters of ditomyiid larvae, such as the presence of the eighth abdominal spiracle and the structure of the mouthparts.
The family Cecidomyiidae is supported as the sister group to the remainder of Sciaroidea (excluding Ditomyiidae) in both the Bayesian and ML analyses (0.97/91). We consider it remarkable that our analyses showed Cecidomyiidae to be firmly nested within the fungus gnats, as discussed above. Interestingly, in the tree based on 18S only, as well as in the tree based on ribosomal genes only, Cecidomyiidae appeared in the clade with two incertae sedis genera (Chiletricha Chandler + Insulatricha Jaschhof) with high support in the ML analysis (BV = 100), indicating possible affinities of these groups. Also, in the ML tree based on the dataset with the saturated third codon positions removed (Fig. 2), Cecidomyiidae and the incertae sedis merged into one clade (albeit poorly supported). It will be very interesting to see how the topology obtained here will change when other unplaced Sciaroidea, as well as Rangomarama (discussed by Jaschhof & Didham (2002) to be the sister group of Cecidomyiidae), can be added to the dataset. The topologies proposed by Matile (1990) and Matile (1997) differ from our results mainly in the positions of Ditomyiidae, Diadocidiidae and Cecidomyiidae. Some previous authors (e.g., Oosterbroek & Courtney, 1995;Wood & Borkent, 1989) considered Cecidomyiidae as a sister group to Sciaridae. This view is not supported by our analyses.
Bolitophilidae is supported as the sister group to the other Sciaroidea (excluding Ditomyiidae and Cecidomyiidae), but with high support only in our Bayesian analysis (1.00/52). This clade was also supported as monophyletic by a number of larval synapomorphies by Fitzgerald (2004: 42) though larval stages of many taxa remain unknown and were therefore not included. Bolitophilidae is a family poorly represented in previous phylogenetic studies, especially molecular ones. As an exception, Wiegmann et al. (2011) included one species of Bolitophila Meigen in their dataset. Their analysis resulted in the sister grouping of Bolitophilidae and Keroplatidae, a topology that differs from ours, indicating the need of further studies into this issue in the future. One of the possible reasons for this discrepancy might be the fact that Bolitophila sp. was represented in their dataset by only three gene markers (28S,TPI, AATS1) and the whole family Keroplatidae only by one rather atypical member, Arachnocampa flava Harrison, 1966. The possible relationship of Arachnocampa to Bolitophilidae is an interesting hypothesis itself, though beyond the scope of this paper, considering the fact that the type species of the genus, Arachnocampa luminosa (Skuse, 1891) was initially described in Bolitophila and removed from this genus subsequently by Edwards (1924), who also discussed similarities and differences of the two genera.
The present study corresponds with an earlier one by Ševčík et al. (2014) which suggested Sciaridae to be the closest relative of Diadocidiidae, although with relatively high support only in Bayesian analysis (0.96/-). However, this clade is better supported (1/74) in the tree without the incertae sedis genera (Fig. S1) and in the ML tree based on the dataset without saturated third codon positions (Fig. 2).
As revealed in this study, the relationships among the families of Sciaroidea still require further testing with the addition of more incertae sedis taxa and with a more extensive gene sampling (e.g., from next-generation sequencing). However, judging from the short branches in the family-level clades, it can be concluded that these taxa underwent very rapid radiation leaving low phylogenetic signal, which will be difficult to resolve. Another serious limitation will always be obtaining fresh specimens with a sufficient amount of DNA/RNA for several key taxa that are extremely rare.

Paraphyly of Keroplatidae
The family Keroplatidae is paraphyletic with respect to Lygistorrhinidae in this study in both BI and ML analyses (the node support for the entire clade is 0.99/47, see Fig. 1) and also in the ML analysis based on the dataset without the saturated third codon positions (Fig. 2). The close relationship of Keroplatidae and Lygistorrhinidae was also indicated in most of the single-gene trees, being best supported in the CAD-based tree.
The sister clade to Arachnocampa, containing all the remaining taxa of Keroplatidae and Lygistorrhinidae, is even better supported (1.00/100). Within the latter clade, the genera of Lygistorrhinidae form a sister group to Keroplatinae, but with low support (0.65/46), see also below ('Keroplatidae and Lygistorrhinidae'). Tuomikoski (1966) was the first to hypothesize a possible relationship between Lygistorrhinidae and Keroplatidae, indicating that Lygistorrhina Skuse, 1890 most probably represents only a specialized group within Keroplatidae. He argued that the reduced wing venation in this genus may be misleading and that the thoracic structure, male terminalia, and other features of Lygistorrhina are very similar to those of Fenderomyia Shaw, 1948 and other taxa of Macrocerinae. However, his views were criticized by Thompson (1975), arguing that some of the characters shared by both groups (e.g., simple type of the male terminalia) are actually symplesiomorphies. Subsequent studies based on morphology have usually placed Lygistorrhinidae close to Mycetophilidae, and Keroplatidae close to Diadocidiidae (Matile, 1990;Matile, 1997;Chandler, 2002;Hippa & Vilkamaa, 2006;Amorim & Rindal, 2007). On the contrary, the close relationship between Lygistorrhinidae and Keroplatidae was indicated by three independent previous molecular studies (Bertone, Courtney & Wiegmann, 2008;Rindal, Søli & Bachmann, 2009;Ševčík et al., 2014) and it has also recently been suggested by Kallweit (2014) based on morphology.
The results of our analyses thus strongly support the inclusion of Lygistorrhinidae in the family Keroplatidae, probably best treated as a subfamily of the latter. The position of Arachnocampa, as discussed above (in 'Monophyly and interrelationships of bibionomorph families'), requires further study, which is beyond the scope of this paper.
These five genera are the only taxa from this enigmatic grouping with molecular data currently available. Our analyses showed Ohakunea to be only distantly related to the other four genera (Fig. 1), even though its position within the Mycetophilidae came as a surprise (see also 'Mycetophilidae'). The other enigmatic genera are united in a single moderately supported clade (0.99/57) and as sister group to the clade Diadocidiidae + Sciaridae, but with lower support (0.95/-). Within this clade, Sciarosoma is branching basally to the well supported (1.00/82) clade of Nepaletricha and Chiletricha + Insulatricha.
As long as the majority of incertae sedis remain unstudied, we think it premature to attach too much importance to the specific positions of the enigmatic genera in the topology presented here, but the well-supported monophyly of (Chiletricha + Insulatricha) + Nepaletricha could be interpreted as giving support to the Chiletrichinae of Amorim & Rindal (2007).
When we excluded all the incertae sedis genera from the dataset (Fig. S1), the support values for the clade Diadocidiidae + Sciaridae increased to 1.00/78 and those for the Keroplatidae clade (including Arachnocampa and Lygistorrhinidae) increased to 1.00/60. This manipulation shows how important it is to include as many of the incertae sedis genera as possible in future molecular analyses.

Anisopodidae
The monophyly of Anisopodidae sensu lato is well supported (1.00/100) in all the analyses (Figs. 1 and 2), as well as the sister relationship of Mycetobia Meigen, 1818 and Mesochria Enderlein, 1910. The topology Olbiogaster + (Sylvicola + Mycetobia) in our analyses is in concordance with the hypothesis based on morphology given by Amorim & Tozoni (1994) while Michelsen (1999) suggested a different relationship, Sylvicola Harris, 1780 + (Olbiogaster Osten-Sacken, 1886 + Mycetobiidae). The latter topology was produced by our datasets with the third codon positions removed (Fig. 2) but with very low support for the clade (Olbiogaster + Mycetobiinae). A comprehensive molecular phylogeny of this group, based on better taxon sampling, is needed to elucidate these relationships.

Bibionidae and Pachyneuridae
The monophyly of Bibionidae sensu lato (Bibionidae sensu stricto + Hesperinidae) was established with high support (1.00/96). The clade containing the genera Bibio Geoffroy, 1762 and Dilophus Meigen, 1803 (Bibioninae) is sister group to the clade Plecia Wiedemann, 1828 + (Penthetria Meigen, 1803 + Hesperinus Walker, 1848), with all nodes well supported (Fig. 1). This topology suggests that either Hesperinus should not be classified in the separate family Hesperinidae because it would render the rest of Bibionidae paraphyletic or, Bibioninae, Plecia, Penthetria, and Hesperinus each need to be treated as distinct families. The position of Hesperinus within (rather than sister group to) Bibionidae is a novel hypothesis (DNA sequences for Hesperinus are published here for the first time) and it conflicts with several previous hypotheses that place Hesperinus as sister group to the remaining Bibionidae (Blaschke-Berthold, 1994;Fitzgerald, 2004;Pinto & Amorim, 2000). The topology found here requires a number of morphological characters that have been found to support bibionids exclusive of Hesperinus (e.g., male eye holoptic, male antennal flagellomeres compressed, larva with fleshy tubercles, larval intersegmental fissures separating abdominal segments seven and eight unaligned laterally, and larval abdominal segment three with three pseudosegments, see Fitzgerald, 2004) to be independently derived multiple times or derived once and then secondarily lost in Hesperinus. Further studies to elucidate these issues are thus warranted. Congruent with previous morphological studies focused on bibionids (Fitzgerald, 2004;Pinto & Amorim, 2000) the clade Plecia + Penthetria was not supported in this analysis thereby further discouraging the recognition of this clade which is commonly classified as Pleciidae or Pleciinae.
Bibionoidea (Bibionidae sensu lato + Pachyneuridae) was strongly supported as a monophyletic group in our analyses (1.00/100) echoing the findings of several earlier morphological (Blaschke-Berthold, 1994;Fitzgerald, 2004) and molecular studies (Bertone, Courtney & Wiegmann, 2008;Wiegmann et al., 2011). Due to the rarity of pachyneurids for study, the monophyly of the family Pachyneuridae sensu lato (Pachyneura Zetterstedt, 1838 + Cramptonomyiinae) has only rarely been assessed (Amorim, 1993;Blaschke-Berthold, 1994;Fitzgerald, 2004) and monophyly of the group has never been tested using molecular characters (Pachyneura is here sequenced for the first time). Fitzgerald (2004: 41) found seven unambiguous characters supporting monophyly of Pachyneuridae s.l., including one adult and three larval characters unique to the family. The present study also found strong support for the monophyly of Pachyneuridae s.l. thereby not supporting the previous hypothesis based on wing venation (Amorim, 1993) that considers Pachyneuridae s.l. to be paraphyletic and treats the genus Pachyneura in Axymyiomorpha with Axymyiidae.

Cecidomyiidae
The monophyly of Cecidomyiidae is here fully supported (1.00/100) in all the analyses (Figs. 1-3, Figs. S1 and S2). All the major lineages (subfamilies) recognised within the family except Winnertziinae are represented in our dataset. The clade Catocha Haliday, 1833 + (Catotricha Edwards, 1938 + Lestremia Macquart, 1826) was found to be the sister branch to the other Cecidomyiidae, which is remarkable insofar as solely Catotrichinae are usually placed at the base of the family (see above). The sister branch is formed of proepimeron and epimeron, and strict endomycophagy of the larvae). This view is also supported by our data (Figs. 1, 2 and Figs. S1, S2).

Sciaridae
A comprehensive molecular phylogeny of this family was recently presented by Shin et al. (2013). The sampling of Sciaridae in our study is necessarily limited to keep the numbers of terminal taxa in proportion for particular families, but the relationships among the genera included in our dataset principally agree with the topology presented in the latter study (Figs. 1, 2 and Figs. S1, S2).

Comments on the selection of molecular markers and perspectives for the future
This study represents the first molecular phylogeny of Bibionomorpha based on a comprehensive taxon sampling. We tried to use easily-amplifiable, traditional gene markers because we were strongly limited by the scarcity of many taxa included in the dataset (and therefore had a shortage of DNA). For the following rare genus-group taxa we publish the first molecular data at all: all the Sciaroidea incertae sedis included (except Nepaletricha), Asioditomyia Saigusa, 1973, Blagorrhina Hippa, Mattison & Vilkamaa, 2009, Catocha, Catotricha, Chiasmoneura Meijere, 1913, Hesperinus, Hyperoscellis Hardy & Nagatomi, 1960, Mesochria, Pachyneura, Protaxymyia Mamaev & Krivosheina, 1966, Robsonomyia Matile & Vockeroth, 1980, Rutylapa Edwards, 1929-and for some clades we thus present the first phylogenetic hypotheses based on molecular data.
We also succeeded in incorporating two regions (of total length more than 1,400 bp) of the protein-coding nuclear gene CAD, which is recommended for use in reconstruction of Diptera phylogeny by Moulton & Wiegmann (2004) and subsequent authors. The amplification of the various fragments of this gene proved to be much more difficult than in traditional markers, so the inclusion of more protein-coding nuclear genes has not been possible with current material, which includes many rare species, available as only one or two specimens, usually collected several years ago. As soon as more fresh material (especially of Sciaroidea incertae sedis) is available for next-gen sequencing, which may take several years or more, it should be possible to design specific primers for Sciaroidea and amplify the supposedly more informative protein-coding nuclear genes, such as MCS and MAC (see Winkler et al., 2015). Our attempts to amplify MCS and MAC markers from most of the specimens of Bibionomorpha using the primers designed for other groups of Diptera failed.
Despite these shortcomings, we consider it important to publish a tree with lower resolution in a few branches because it contributes to our understanding of the performance of these traditional markers across the various taxa of Bibionomorpha, especially when they provide sufficient resolution for most other relationships within the group. Our dataset is largely based (74% of base pairs) on the very little saturated nuclear 18S, 28S and CAD genes, more suitable for deeper relationships, in combination with three short mitochondrial regions included to improve the resolution of terminal branches, and in the case of COI also to enrich the DNA barcode database.

CONCLUSIONS
1. Monophyly of Bibionomorpha is supported in both the strict sense (Sciaroidea + Bibionoidea) and the broad sense (Bibionomorpha sensu strico + Anisopodoidea + Scatopsoidea). Axymyiidae is not resolved as part of, or as the sister group to, Bibionomorpha. 2. The monophyly of the groups Bibionoidea (Bibionidae sensu lato + Pachyneuridae sensu lato), Bibionidae sensu lato (Bibionidae sensu stricto + Hesperinidae), and Pachyneuridae (Pachyneura + Cramptonomyiinae) were all established with high support in all the analyses. Based on the topology recovered here, treating Hesperinus as the distinct family Hesperinidae would render the rest of Bibionidae paraphyletic. Diverging from previous hypotheses that place Hesperinus as the sister group to the remaining bibionids, the position of Hesperinus in this study is a novel hypothesis that requires further testing. As in previous studies, Pleciidae or Pleciinae (Plecia + Penthetria) was not supported as a monophyletic group and its usage should be discontinued. 3. The major lineages of Bibionomorpha sensu lato (Sciaroidea, Bibionoidea, Anisopodoidea, and Scatopsoidea) are well supported as monophyletic groups and all included families are supported as monophyletic with the exception of Keroplatidae and Bibionidae sensu lato. Both Cecidomyiidae and Ditomyiidae were recovered in model-based analyses as a part of Sciaroidea with high support. All the other families previously placed in Sciaroidea were confirmed to belong to this monophyletic group. The Keroplatidae clade (including Lygistorrhinidae) is well supported, but not so much with Arachnocampa. In most other cases support is too low to confidently resolve family level relationships within Sciaroidea. 4. The five genera of Sciaroidea incertae sedis that were included in this analysis do not constitute a monophyletic group thereby not supporting an expanded concept of Rangomaramidae (Amorim & Rindal, 2007) but further studies with more taxa included are much needed.