Systellognatha is a monophyletic group within stoneflies containing seven families in two superfamilies. However, the family–level phylogeny within Systellognatha is largely unresolved, with some controversial hypotheses generated by molecular data. In this study, we sequenced one complete and three nearly complete mitochondrial genomes (mitogenomes) representing three systellognathan families (including one newly described family), to recover the phylogeny of all seven extant systellognathan families. We tested the effect of heterogeneous sequence divergence within mitogenome sequences on phylogenetic reconstruction. To study incongruent topologies supported by different datasets, we also utilized four–cluster likelihood mapping to elucidate these phylogenetic questions. We found a low degree of compositional heterogeneity in systellognathan mitogenomes. Our results showed that phylogenetic analyses with various datasets under site–homogeneous models failed to recover the Perloidea and Pteronarcyoidea as monophyletic groups. However, Bayesian inference with a site–heterogeneous model (CAT + GTR) recovered Perloidea as monophyletic. Our results confirm a sister–group relationship between Chloroperlidae and Perlodidae. Despite all phylogenetic analyses failed to recover the Pteronarcyoidea as monophyletic, the analyses with site–heterogeneous model based on PCG12, PCG12R and PCGR datasets reconstructed the same relationships among Pteronarcyidae, Peltoperlidae, Styloperlidae and Perloidea. Our study shows the power of the site–heterogeneous models in resolving phylogeny within Systellognatha and provides a new insight into the phylogeny of this group. [ABSTRACT FROM AUTHOR]