Comparative mitogenomic analysis of subterranean and surface amphipods (Crustacea, Amphipoda) with special reference to the family Crangonyctidae

Abstract

Mitochondrial genomes play important roles in studying genome evolution, phylogenetic analyses, and species identification. Amphipods (Class Malacostraca, Order Amphipoda) are one of the most ecologically diverse crustacean groups occurring in a diverse array of aquatic and terrestrial environments globally, from freshwater streams and lakes to groundwater aquifers and the deep sea, but we have a limited understanding of how habitat influences the molecular evolution of mitochondrial energy metabolism. Subterranean amphipods likely experience different evolutionary pressures on energy management compared to surface-dwelling taxa that generally encounter higher levels of predation and energy resources and live in more variable environments. In this study, we compared the mitogenomes, including the 13 protein-coding genes involved in the oxidative phosphorylation (OXPHOS) pathway, of surface and subterranean amphipods to uncover potentially different molecular signals of energy metabolism between surface and subterranean environments in this diverse crustacean group. We compared base composition, codon usage, gene order rearrangement, conducted comparative mitogenomic and phylogenomic analyses, and examined evolutionary signals of 35 amphipod mitogenomes representing 13 families, with an emphasis on Crangonyctidae. Mitogenome size, AT content, GC-skew, gene order, uncommon start codons, location of putative control region (CR), length of rrnL and intergenic spacers differed between surface and subterranean amphipods. Among crangonyctid amphipods, the spring-dwelling Crangonyx forbesi exhibited a unique gene order, a long nad5 locus, longer rrnL and rrnS loci, and unconventional start codons. Evidence of directional selection was detected in several protein-encoding genes of the OXPHOS pathway in the mitogenomes of surface amphipods, while a signal of purifying selection was more prominent in subterranean species, which is consistent with the hypothesis that the mitogenome of surface-adapted species has evolved in response to a more energy demanding environment compared to subterranean amphipods. Overall, gene order, locations of non-coding regions, and base-substitution rates points to habitat as an important factor influencing the evolution of amphipod mitogenomes.

Introduction

Caves and other subterranean habitats, such as groundwater aquifers and superficial subterranean habitats (SSHs; [1]), represent some of the most challenging environments that exist on Earth. The primary characteristic of all subterranean habitats is the lack of light and associated photosynthesis [1, 2]. Though some subterranean ecosystems are supported by chemoautotrophic production by microbial communities [3, 4], chemoautotrophy rarely provides enough energy to support several trophic levels in most subterranean ecosystems [1, 5]. The primary source of energy input for many cave systems is the organic matter transferred from the surface hydrologically or by animals that frequently enter and exit caves [1, 6], which drive the structure and dynamics of subterranean communities [7,8,9]. Although most subterranean ecosystems are largely thought to be energy-limited [10], food availability can be highly variable both among and within cave systems [11, 12]. Previous studies have shown that many subterranean organisms living in such energy-limited habitats have undergone several physiological and metabolic adaptations to sustain themselves during extended food shortages [13, 14]. Among these troglomorphic traits, low metabolic rate is a key adaptation that occurs in both terrestrial and aquatic fauna of subterranean communities [15, 16].

Mitochondria are the primary sites of energy production in cells, generating ~ 95% of the adenosine triphosphate (ATP) required for everyday activities of life through oxidative phosphorylation [17,18,19]. The mitochondrial genome—mitogenome—encodes 13 essential proteins including two ATP synthases (atp6 and atp8), three cytochrome oxidases (cox1, cox2, and cox3), seven NADPH reductases (nad1, nad2, nad3, nad4, nad4l, nad5, and nad6), and cytochrome b (cytb) subunits. All mitochondrial protein-coding genes (PCGs) play a vital role in the electron transport chain [20,21,22]. Due to the unique characteristics of mitochondria, including maternal inheritance, small genomic size, absence of introns, and their surplus availability in cells, the use of mitochondrial DNA (mtDNA) loci and mitogenomes has a long history in population genetics, phylogenetics, and molecular evolution studies [23,24,25]. Previous studies have demonstrated a close association between mitochondrial loci and energy metabolism [18, 26, 26, 27]. Although considered to largely evolve under purifying selection, there is growing evidence that mitogenomes may undergo episodes of directional selection in response to shifts in physiological or environmental pressures [28, 29] leading to improved metabolic performance under new environmental conditions [26, 30, 31]. For example, previous studies that investigated varying selective pressures acting on mitochondrial PCGs of insects and mammals have revealed significant positive selective constraints at several loci that have comparatively increased energy demands [18, 19, 32]. Similarly, other studies have shown the various adaptive mitochondrial responses of organisms surviving in extreme environments including the deep sea and Tibetan Plateau [29, 32, 33]. However, these adaptations can occur at different metabolic levels, not just mitochondrial metabolism [34, 35]. Thus, variation in mitogenomes of species inhabiting different environments may reflect only a small portion of these adaptive metabolic changes. Despite this limitation, previous studies have detected signals of directional selection in the mitogenomes of organisms dwelling in contrasting habitats with varying energy demands [36,37,38].

Amphipods (Class Malacostraca: Order Amphipoda) are one of the most ecologically diverse crustacean groups including over 10,000 species [39, 40], occurring in a diverse array of aquatic and even terrestrial environments globally, from aphotic groundwater aquifers and hadal depths to freshwater streams and lakes in temperate and tropical forests, among other habitats [41, 42]. Several studies have demonstrated the genetic basis of subterranean adaptation in several taxa, including dytiscid diving beetles [43], cave dwelling-cyprinid fishes [44, 45], anchialine cave shrimps [46], and cave isopods [47]. However, we still have a limited understanding of the mechanisms of subterranean adaptations in amphipods. Although physiological adaptations to challenging environments like cave and groundwater ecosystems have been well-studied in amphipods [13, 16], no studies to date have addressed the selective pressures and the molecular evolution mechanisms of mitochondrial energy metabolism loci in amphipods occupying caves and other subterranean habitats. Subterranean amphipods likely experience different evolutionary pressures on energy management due to lower levels of predation, lower food resources, and more stable environments compared to surface-dwelling taxa that generally experience higher levels of predation and energy resources [48, 49].

In this study, we compared the mitogenomes of surface and subterranean amphipods, including the 13 mitochondrial PCGs involved in the OXPHOS pathway to understand the potential molecular mechanisms of energy metabolism in this diverse crustacean group. Our aims were to test whether the mitochondrial PCGs showed evidence of adaptive evolution in subterranean environments in amphipods. We tested the hypothesis that the mitogenome of surface-adapted amphipods will be imprinted by mitogenomic adaptations to the energy demanding environment with greater signal of directional selection when compared to their subterranean counterparts. We compared base composition, codon usage, gene order rearrangement, conducted comparative mitogenomic and phylogenomic analyses, and examined evolutionary signals using publicly available amphipod mitogenomes. In particular, we focused on the amphipod family Crangonyctidae, a diverse family that comprises species inhabiting a variety of surface and subterranean habitats and for which several mitogenomes have been sequenced and annotated recently [50, 51].

Materials and methods

We generated new mitogenomes recently for the following crangonyctid species: Stygobromus pizzinii, S. tenuis potomacus, Bactrurus brachycaudus, Stygobromus allegheniensis, and Crangonyx forbesi [51].

DNA extraction, library preparation, and sequencing

Whole genomic DNA from five crangonyctid species was isolated using the Qiagen DNA Easy Blood and Tissue kit and libraries were prepared using the Illumina TruSeq DNA Library Prep Kit (Illumina Inc., California). Libraries were then paired-end sequenced (2 × 150 bp) on an Illumina HiSeq 4000 platform at the Vincent J. Coates Genomics Sequencing Laboratory at the University of California, Berkeley. We assessed the quality of the raw reads using FastQC v0.11.5 [52], and the reads were trimmed and filtered using Trimmomatic v0.33 [53]. De-novo assembly was carried out using NOVOPlasty v2.6.3 assembler [54]. We then annotated the protein-coding genes, transfer RNAs (tRNAs), and ribosomal RNAs (rRNAs) for each of the five mitogenomes using the mitochondrial genome annotation web server MITOS [55]. The secondary structures of tRNAs were inferred using MITFI [56], a built-in module in MITOS. The locations of start and stop codons of protein coding genes were confirmed using NCBI ORFfinder [57] and by visual comparison to other published amphipod mitogenomes. The location of the control region was confirmed by the presence of a large intergenic spacer region with a string of thymines found immediately after rrnS and before trnl. We then downloaded from GenBank the annotated mitogenomes of 30 additional amphipod taxa that occupy aquatic habitats, including groundwater and springs, and three isopods that served as outgroups for comparative analyses.

Mitogenome analyses

Nucleotide composition, amino acid frequencies, and codon usage were calculated in PhyloSuite v1.1.15 [58, 59]. The web-based program CREx (http://pacosy.informatik.uni-leipzig.de/crex, [60]) was used to perform pair-wise comparison of the gene orders in the mitogenome to determine rearrangement events. CREx comparisons were based on common intervals, and it considers common rearrangement scenarios including transpositions, reversals, reverse transpositions, and tandem-duplication-random-losses (TDRLs). In addition, phylograms and gene orders were visualized in iTOL (https://itol.embl.de/, [61]) using files exported from PhyloSuite. AT and GC skew of entire mitogenomes and individual loci were calculated in PhyloSuite using the formulae: AT-skew = (A – T)/(A + T) and GC-skew = (G – C)/(G + C). Welch two sample t-tests were performed between the surface and subterranean amphipods for different mitogenomic features, including mitogenome length, AT content, AT and GC skew, and rRNA length using R [62]. Visualization of AT-skew, GC-skew, AT-content, and amino acid frequencies were generated in R.

Phylogenetic inference

The amino acid sequences of 13 PCGs of the five new mitogenomes [51], 30 previously published amphipod mitogenomes, and three isopod mitogenomes were aligned using MAFFT version 7 [63]. A total of 38 sequences with 350 positions in the alignment file were trimmed using Gblocks version 0.91b [64] to yield 255 positions (72%) in 6 selected blocks (parameters used: Supplementary Table S1). The alignment was partitioned by gene and then the best-fit partitioning strategy and evolutionary models for each partition were inferred using PartitionFinder v2.1.1 [65]; Supplementary Table S2). Phylogenetic relationships of the 35 amphipod mitogenomes and three isopod mitogenomes using the concatenated 13 PCG alignment were determined using Bayesian inference in MrBayes v3.2 [66]. The analyses contained two parallel runs with four chains each and were conducted for 5,000,000 generations (sampling every 100 generations). After dropping the first 25% “burn in” trees to ensure stationarity after examination of log-likelihood values for each Bayesian run using Tracer v1.7 [67], the remaining 37,500 sampled trees were used to estimate the consensus tree and the associated Bayesian posterior probabilities. All analyses were conducted within PhyloSuite.

Positive selection and selection pressure analyses of mitochondrial PCGs

We performed base-substitution analyses on entire mitogenomes as well as for each of the 13 PCGs individually to compare surface versus subterranean amphipod taxa. The non-synonymous to synonymous rate ratio (dN/dS or ω) for each PCG was estimated using the free-ratio model using the CodeML program implemented in PAML v4.8a [68]. The ω values were estimated for surface and subterranean species separately and visualized in R for comparison. To estimate the probability of positively selected sites in each PCG across all amphipod species, we implemented site models (M1 and M2, M8a and M8), where ω was allowed to vary among sites [69]. To further investigate if some lineages and sites in particular lineages have undergone positive selection, we conducted maximum likelihood analyses on all PCG using the branch model and branch-site model in EasyCodeML v1.21, a visual tool for analysis of selection using CodeML [70].

To determine if all 13 PCG are free of functional constraints in subterranean lineages, we compared alternative branch selection models on each PCG tree. First, we tested a model (M0) where a single ω was estimated for all branches. This model was compared to a model (M1) with two ratios, a background ω for surface branches and a separate foreground ω for subterranean branches. In addition, we included a two-ratio model where ω was fixed at 1.0 in subterranean branches (M1fixed) to determine if estimates of ω differed from rates of neutral evolution and a model similar to the M1 model but where each subterranean lineage (B. brachycaudus, B. jaraguensis, P. daviui, and the clade containing Stygobromus and Metacrangonyx) was allowed to have a separate ω (M1a). We also examined a saturated model (M2) where each branch had its own ω. Akaike’s information criterion (AIC) was used to compare models.

For both the branch and branch-site models, a likelihood ratio test (LRT) was conducted for each PCG to test whether ω was homogeneous across all branches. In the branch model, the null hypothesis assumes that the average ω values of branches of interest (ωF) is equal to that of other branches (ωB), whereas the alternative hypothesis assumes the opposite ωF ≠ ωB. If the alternative hypothesis is supported and ω > 1, the foreground lineage is under positive selection. The branch-site model allows ω to differ among codon sites in a foreground lineage when compared to background lineages. We implemented the branch-site model to identify sites on specific lineages regulated by positive selection. Selected sites were considered positively selected only if they passed a Bayes Empirical Bayes (BEB) analysis with a posterior probability of > 95%.

We performed selection pressure analyses on the concatenated 13 PCGs dataset aligned using the codon mode as well as on each PCG with the Bayesian topology (see Fig. 4) as the guidance tree using several approaches available from the Datamonkey Webserver [71]. First, we implemented aBSREL (Adaptive Branch-Site Random Effects Likelihood), an improved version of the commonly used “branch-site” models, to test if positive selection has occurred on a proportion of branches [72]. We implemented BUSTED (Branch-site Unrestricted Statistical Test for Episodic Diversification) to test for gene-wide (not site-specific) positive selection by querying if a gene has experienced positive selection in at least one site on at least one branch [73]. Finally, we implemented RELAX [74] to test whether the strength of selection has been relaxed or intensified along a specified set of test branches.

Results and discussion

We compared the mitogenomes for 35 surface and groundwater amphipods, including mitogenomes of one spring-dwelling and six groundwater species in the family Crangonyctidae by Aunins et al. [50] and Benito et al. [51], to determine whether subterranean species show evidence of adaptive evolution in subterranean habitats. Our study examined whether features of mitogenomes (e.g. base composition, codon usage, gene order) differed in relation to dominant habitat (surface vs. subterranean) and inferred the evolutionary forces potentially shaping mitogenome evolution in amphipods, with an emphasis on crangonyctid species.

Mitogenome length and content

Mitogenome sizes ranged from 14,113 to 18,424 bp for all amphipods and 14,661 to 15,469 bp for crangonyctid amphipods (Table 1). Mean mitogenome size of surface amphipods (15,770 ± 1206 bp; mean ± 1 standard deviation) was higher than that of the subterranean amphipods (14,716 ± 297 bp) (phylogenetic paired t-test: t = 0.586, df = 33, p-value = 0.562; Supplementary Figure S1). All crangonyctid amphipod mitogenomes possessed 13 PCGs, two rRNA genes, 22 tRNA genes, a control region, and intergenic spacers of varying number and lengths (Supplementary Figure S2, annotations of the genomes are presented in Supplementary Table S3) like other arthropods [75]. The length of the crangonyctid mitogenomes was similar to lengths reported for other amphipod families including Allocrangonyctidae, Caprellidae, Eulimnogammaridae, Gammaridae, Hadziidae, Lysianassidae, Metacrangonyctidae, Micruropodidae, Pallaseidae, Pontogeneiidae, Talitridae. Variation in mitogenome length within Crangonyctidae appears to be related to length variation in the nad5, rrnL, and rrnS loci, which were all notably longer in the C. forbesi mitogenome.

Table 1 Summary of mitogenomic characteristics, location, and habitat of subterranean and surface amphipods included for comparative mitogenome analyses

Base composition and AT- and GC-skews

Mitogenome AT% in all amphipods ranged from 62.2 to 76.9% (Table 1). Mean AT% of the subterranean amphipods (71.8 ± 3.6%) was higher than that of the surface amphipods (67.6 ± 3.4%) (phylogenetic paired t-test: t = 0.926, df = 33, p-value = 0.361). Mean AT% of all 13 PCG of the subterranean amphipods was significantly higher than that of the surface amphipods (Supplementary Figure S3a). Variation in AT% across crangonyctid amphipod taxa ranged 63.9–69.3%, with a mean of 67.9 ± 1.93% (Table 1). Across loci, AT% ranged from a minimum of 60.0% at the cox1 locus and a maximum of 75.5% at the nad4l locus (Fig. 1A). Variation in AT% across all PCGs combined ranged from 61.9% (B. brachycaudus) to 69.0% (S. indentatus). Genes encoded on the negative strand had a slightly higher AT-content values than those on the positive strand. The nad6 locus showed the greatest variation in AT-content across species. Bactrurus brachycaudus displayed the outlier lower AT% values for most of the PCG (Table 2; Fig. 1A). Similarly, Bactrurus brachycaudus had the lowest AT content (63.9%) among the crangonyctid mitogenomes, while all other mitogenomes had comparatively typical AT content reported for other arthropods [86, 87]. This could indicate that the evolution of the B. brachycaudus mitogenome has occurred under different evolutionary pressures (adaptive or non-adaptive) than other subterranean crangonyctids.

Fig. 1

Crangonyctidae mitochondrial nucleotide composition. Box plots showing values of nucleotide composition (A + T percentage) (a), AT-skew (b), and GC-skew (c) across mitogenomes, protein coding genes (PCG), and ribosomal (rRNA) and transfer ribosomal (tRNA) RNA. The same features are shown for each protein-coding gene and pooled by codon position and coding strand. Genes coded on the (-) strand are represented by a “-“ sign and genes coded on the (+) strand are represented by “+” sign at the end of the gene label

Table 2 Comparison of mitogenomic characteristics of 35 amphipods discussed in this study

Mitogenome AT-skew in all amphipods ranged from − 0.062 to -0.037. Mean AT-skew of the surface amphipods (0.001 ± 0.02) was positive and slightly higher than that of the subterranean amphipods (-0.004 ± 0.02) (phylogenetic paired t-test: t = 0.045, df = 33, p-value = 0.965). Mean AT-skew of four PCG (cox1, cox2, nad2, nad3) of surface amphipods was significantly higher than that of the subterranean amphipods, whereas the mean AT-skew of nad4 of the subterranean amphipods was significantly higher than that of the surface amphipods (Supplementary Figure S3b). Among crangonyctid amphipods, mean AT-skew was 0.022 ± 0.02 (range 0.004 to 0.061), with all mitogenomes displaying positive skew. Mitogenome GC-skew ranged from − 0.431 − 0.120. Mean GC-skew of the subterranean amphipods (-0.129 ± 0.15) was negative and higher than that of the surface amphipods (-0.236 ± 0.05) (phylogenetic paired t-test: t = 0.349, df = 33, p-value = 0.729). Mean GC-skew of seven PCG (atp6, atp8, cox1, cox2, cox3, nad2, nad3) of subterranean amphipods was significantly higher than that of the surface amphipods, whereas the mean GC-skew of nad4 of the surface amphipods was significantly higher than that of the subterranean amphipods (Supplementary Figure S3c). Among crangonyctid amphipods, mean GC-skew was − 0.264 ± 0.01 (range − 0.275 to -0.248) with all mitogenomes displaying negative skew (Table 1). Strand asymmetry is commonly observed in mitogenomes [88, 89], however, at times it can hinder phylogenetic reconstruction and yield false phylogenetic artefacts especially when unrelated taxa display inverted skews [90, 91]. Bactrurus brachycaudus exhibited the lowest AT skew among the crangonyctid mitogenomes (0.004), while S. tenuis had the lowest GC skew (− 0.275). Crangonyctid amphipod mitogenomes exhibited positive GC-skew values for genes encoded on the (-) strand and negative GC-skew for genes encoded on the (+) strand (Fig. 1C), whereas all PCGs exhibited negative AT-skew values (Fig. 1B). Except the six loci (nad1, nad4, nad4L, nad5, rrnL, and rrnS) which were encoded on the (-) strand, most PCG had negative GC skews. Such strand bias is typical for most mitochondrial genomes in metazoan [81, 83]. This is consistent with the hypothesis that strand asymmetry is caused by spontaneous deamination of bases in the leading strand during replication [88]. All other mitogenomes had comparatively typical AT and GC skew values like other amphipod species [76, 92]. The only outlier to this pattern was the positive GC skew value of tRNAs encoded on the (+) strand of B. brachycaudus (0.012). In general, crangonyctid mitogenomes exhibited relatively consistent skews.

Rearrangements of mitochondrial genome

Comparisons of crangonyctid mitogenomes revealed at least six conserved gene blocks (Fig. 2B). The gene orders in subterranean species (genera Stygobromus and Bactrurus) are identical except for the transposition of tRNA-G,W. However, a few unique gene order arrangements were observed in the spring-dwelling C. forbesi. The gene order of C. forbesi differs from the four subterranean species in the locations of the conserved gene blocks (tRNA-H-nad4-nad4l and nad6-cytb-tRNA-S2 and tRNA-L1-rrnL and rrnS-tRNA-I and tRNA-Y,Q), seven tRNAs (P,T,M,V,G,C, and W), and two protein-coding loci: nad1 and nad2. Compared to the conserved mitogenome gene orders of other crangonyctid mitogenomes, another unique feature in the rearranged C. forbesi mitogenome was the presence of at least two long (~ 50 and 70 bp) non-coding regions (Supplementary Table S3). The locations of rRNA genes in all crangonyctid mitogenomes are mostly similar compared to the pancrustacean ground pattern except for C. forbesi where the rRNA genes had altered positions (Fig. 2A and B). Rearrangements in the mitogenome is common especially when it involves only tRNA-coding genes [93]. In case of ribosomal RNA genes or PCGs, rearrangements occur much less frequently, and they are commonly referred to as major rearrangements, as they might potentially affect the differential regulation of replication and transcription of mitogenomes [94].

Fig. 2

Mitochondrial phylogenomics and gene orders: (a) Bayesian phylogram inferred using amino acid sequences of all mitochondrial PCGs (left) and gene orders (right). Three isopod outgroups are not shown. GenBank accession numbers are included as suffix next to the species names; (b) gene orders of mitochondrial genomes in three genera of crangonyctid amphipods, including Stygobromus, Bactrurus, and Crangonyx. Conserved gene clusters are indicated by different colors and gene rearrangements are highlighted by red border lines

CREx analysis indicated that transpositions and TDRL may have been responsible for the evolution of mitogenomes in crangonyctid amphipods. Two transpositions of tRNA-R,N,S1,E and two steps of TDRL from the ancestral pan-crustacean pattern were needed to generate the gene order observed in Stygobromus species. In addition to the same two transpositions, one TDRL, and a transposition within a second TDRL from the ancestral pattern were required to generate the gene order in Bactrurus. However, four different transpositions (tRNA-N,S1, tRNA-T,P, tRNA-W,C and gene block tRNA-H-nad4-nad4L-tRNA-P,T-nad6-cytb-tRNA-S2) and three steps of TDRL from the ancestral pattern were needed to generate the gene order observed in C. forbesi (Supplementary Figure S4).

Similar to C. forbesi, other surface amphipods including Gmelinoides fasciatus (Micruropodidae) and Onisimus nanseni (Lysianassidae) exhibited a highly rearranged gene order. Other surface amphipods that exhibited a moderate to highly rearranged gene order include Gondogeneia antarctica (Pontogeneiidae), Platorchestia parapacifica and P. japonica (Talitridae), Pallaseopsis kessleri (Pallaseidae), and the two basal amphipods Caprella scaura and C. mutica (Caprellidae) (Fig. 2A). Interestingly, a subterranean amphipod Pseudoniphargus daviui (Allocrangonyctidae) also exhibited a moderate rearranged gene order. The stark contrast between the highly conserved gene order in most subterranean amphipods and the highly volatile gene order in many of the surface amphipods may support the hypothesis that evolution of mitogenomic architecture could be highly discontinuous. A long period of inactivity in gene order and content could have been interspersed by a rearrangement event, this destabilized mitogenome is much more likely to undergo subsequent accelerated rate of mitogenomic rearrangements [95]. Thus, it is appealing to examine mitogenomes of surface amphipod families represented by just a single taxon in our dataset.

Codon usage and amino acid frequencies

In addition to the regular start codons (ATA and ATG) and uncommon start codons (ATT, ATC, TTG, and GTG), surface amphipods, particularly Caprella scaura, possessed one rare start codon CTG, whereas subterranean amphipods possessed three rare start codons including CTG, TTT, and AAT to initiate the mitochondrial PCGs (Supplementary Table S4). Codon usage analysis of the five crangonyctid amphipods mitogenomes identified the existence of all codon types typical for any invertebrate mitogenome. In addition to the regular start codons (ATA and ATG), uncommon start codons (ATT, ATC, TTG, and GTG) were also present to initiate the mitochondrial PCG. Such unusual start codons have been reported previously in other arthropods [96, 97]. A few PCG in the crangonyctid mitogenomes possessed truncated or incomplete stop codons (TA- and T–) that have been described in other crustaceans (Supplementary Table S3). These are presumably completed after a post-transcriptional polyadenylation [98,99,100]. Among the crangonyctid mitogenomes, the most frequently used codons are TTA (Leu2; 5.64–8.49%) and TTT (Phe; 5.94–6.78%). Other frequently used codons include ATT (Ile; 4.92–6.85%) and ATA (Met; 4.13–5.34%) (Supplementary Table S5). These four codons are also among the most abundant in non-crangonyctid amphipods included in this study. This bias towards the AT-rich codons is quite typical for arthropods [101]. Among crangonyctid amphipod mitogenomes, relative synonymous codon usage (RSCU) values, which is the measure of the extent that synonymous codons depart from random usage, showed a high prevalence of A or T nucleotides at third codon positions (Fig. 3). This trend was also observed in other subterranean and surface amphipods. This positive correlation between RSCU and AT content at third codon positions has been reported in mitochondrial genomes of the abalone and oyster [102,103,104].

Fig. 3

The relative synonymous codon usage (RSCU) of the mitogenome of all crangonyctid amphipods. The RSCU value are provided on the Y-axis and the codon families are provided on the X-axis

In PCGs, the second copy of leucine (8.86–10.01%) and cysteine (0.95–1.17%) are the most and the least used amino acids, respectively. Amino acid frequency analysis of both surface and subterranean amphipods indicated that five amino acids (leucine, phenylalanine, isoleucine, methonine, and valine) account for more than half of the total amino acid composition and exhibited greater variation among species (Supplementary Figure S5; Supplementary Table S6).

Transfer RNA genes

All 22 tRNA genes were identified in the mitogenomes of crangonyctid amphipods. However, the locations of tRNA genes were highly variable among these mitogenomes, and they also displayed altered positions relative to the pancrustacean ground pattern (Fig. 2; Supplementary Figure S3). The secondary structures of all mitogenome-encoded tRNAs belonging to crangonyctid amphipods were predicted and ranged in length from 50 to 66 bp. Most of the tRNAs displayed the regular clover-leaf structures, however, a few displayed aberrant structures. The tRNA-Ser1 (UCU) lacked the DHU arm in all crangonyctid species. Similarly, the tRNA-Ser2 (UGA) lacked the DHU arm in all crangonyctid species except S. allegheniensis where tRNA-Ser2 (UGA) possessed the DHU arm. The DHU arm was also missing in the tRNA-Cys and tRNA-Arg of B. brachycaudus and tRNA-Arg of C. forbesi. The tRNA-Gln lacked the TψC arm in all crangonyctid species except C. forbesi where tRNA-Gln possessed the TψC arm. In addition to lacking the TψC arm, tRNA-Gln of B. brachycaudus lacked the acceptor stem as well (Supplementary Figure S6). The presence of aberrant structures in tRNAs have been observed in several other crustaceans and invertebrates [79, 105,106,107], which may be the result of replication slippage [108] or selection towards minimization of the mitogenome [109].

Ribosomal RNA genes

The length of rrnL genes in all amphipods ranged 976–1,137 bp and that of rrnS genes ranged 618–1,631 bp. rrnL length of the subterranean amphipods (1,055 ± 26 bp) was higher than that of the surface amphipods (1005 ± 46 bp) (phylogenetic paired t-test: t = 0.921, df = 33, p-value = 0.364). rrnS length of the surface amphipods (738 ± 258 bp) was higher than that of the subterranean amphipods (684 ± 16 bp) (phylogenetic paired t-test: t = -0.558, df = 33, p-value = 0.581). The length of rrnL genes in crangonyctid amphipods ranged 1,034–1,090 bp and that of rrnS genes ranged 671–695 bp. The length of rRNA genes in crangonyctid amphipods was similar to that of other amphipod mitogenomes except C. forbesi, which had long overhangs (~ 50 bp and ~ 25 bp) on the 5’ end of the rrnL and rrnS genes, respectively. AT content ranged 67.8–72.8% in the rrnL genes and 71.5–77.2% in the rrnS genes of crangonyctid species, respectively. GC-skew values for rRNA genes were positive (0.259 to 0.426) and comparable to that of PCGs encoded on the (-) strand (Supplementary Table S7).

Control region and intergenic spacers

In the mitogenome of S. pizzinii the putative control region (CR) was identified as a 1,021 bp sequence between the rrnS gene and the trnl-trnM-trnC-trnY-trnQ-nad2 gene cluster. Similarly, CR was observed in the other crangonyctid mitogenomes, including S. tenuis (556 bp), S. allegheniensis (991 bp), B. brachycaudus (531 bp), S. indentatus (535 bp), and S. tenuis potomacus (773 bp). The CR was similarly located between the rrnS and nad2 genes in some of the other mitogenomes of non-crangonyctid amphipods, including G. duebeni [81], O. nanseni [83], G. antarctica [84], P. daviui [77], and for the pancrustacean ground pattern. However, the adjacent tRNA genes were often different. In G. fasciatus, the CR region was located between the rrnS and nad5 genes [76]. In contrast, a control region 843 bp was observed in C. forbesi which is located between the nad1 and trnM-trnV-nad2 gene cluster and separated by few intergenic spacers was identified as the CR (Supplementary Figure S2; Supplementary Table S3). The only other surface amphipod that had a similar CR location to C. forbesi was P. kessleri with the CR located between nad1 and nad2 genes, although the adjacent tRNA genes were different [76]. Thus, the variable location of the CR in C. forbesi was in concordance with a few surface amphipods, while the subterranean amphipods mostly followed the universal pattern between rrnS and nad2 genes.

The non-coding regions or intergenic spacers identified in the crangonyctid mitogenomes varied in number and length. The number of intergenic spacers ranged from 7 to 17 and their lengths ranged from 1 to 99 bp (mean 13.0 bp ± 18.6). Two crangonyctid mitogenomes (S. allegheniensis and C. forbesi) possessed the largest intergenic spacers (a total of 220 and 249 bp, respectively; Supplementary Table S3). Among the non-crangonyctid amphipods, G. fasciatus and G. antarctica possessed relatively large non-coding intergenic spacers (a total of 3,863 bp and 4,354 bp, respectively; [76, 84].

Phylogenetic inference

The phylogenetic analyses of the 13 concatenated PCG from 35 amphipod species using Bayesian Inference (BI) resulted in a well-supported phylogeny, with the crangonyctid species forming a well-supported monophyletic group (Fig. 4). Within Crangonyctidae, Stygobromus species formed a monophyletic group sister to Bactrurus + Crangonyx; however, few crangonyctid taxa were included in our analysis. A previous study based on the cox1 gene found that Stygobromus was not monophyletic, but several relationships had low support [110]. Likewise, Stygobromus was recovered as polyphyletic in a multilocus concatenated phylogenetic analysis of the Crangonyctidae by Copilaş-Ciocianu et al. [111]. In addition, several well-supported clades were recovered within Crangonyctidae but relationships among these clades had low support. Other past studies have not supported monophyly of the widespread genera (i.e., Crangonyx, Stygobromus, and Synurella) in the family based on either morphological [112] or molecular data [113, 114]. A comprehensive phylogenomic study with robust taxonomic sampling is greatly needed to better elucidate evolutionary relationships and test biogeographic and ecological hypotheses regarding the origin and diversification of this diverse family of subterranean and surface-dwelling amphipods.

Fig. 4

Bayesian phylogeny of aligned protein-coding loci (3,607 amino acids) for five new amphipod mitogenomes (Stygobromus allegheniensis, S. pizzinii, S. tenuis potomacus, Bactrurus brachycaudus, and Crangonyx forbesi) in addition to 30 additional amphipod mitogenomes available on Genbank. The three isopods Ligia oceanica, Limnoria quadripunctata, and Eophreatoicus sp.14 FK-2009 are included as an outgroup to root the phylogeny. New mitogenomes generated in this study are highlighted. GenBank accession numbers are included as suffix next to the species names. Values at nodes represent posterior probabilities

Selection in PCGs

Most of the energy required for active movement to escape predation and meet energy demands is supplied by the mitochondrial electron transport chain [99, 100. Mitochondrial genes encode for all of the protein complexes related to oxidative phosphorylation except for succinate dehydrogenase (complex II) [115,116,117]. Because of their importance in oxidative phosphorylation during cellular respiration, it is unsurprising that many studies have shown evidence of purifying (negative) selection in mitochondrial PCG [29, 118, 119]. While we found strong evidence for purifying selection in amphipod mitochondrial PCGs in our selection analyses, we also found signatures of positive selection in a few of the mitochondrial PCGs in the surface amphipods.

Using a free-ratio model (M2; [27]), we calculated the ω values for the 13 PCGs for the terminal branches to estimate the strength of selection between different primary habitats (i.e., subterranean vs. surface). The cox2 locus significantly differed in ω values between the amphipods of the two habitat types (p = 0.020), with higher ω values for the surface amphipods. Similarly, cox1 and cox3 genes also exhibited a similar trend (p = 0.095 and p = 0.057, respectively) (Fig. 5). This could be because the rate at which slightly deleterious mutations (ω) responsible for the mitochondrial gene evolution accumulates comparatively faster in cox gene family of the surface lineages. However, this result is quite contradictory to previous studies showing higher functional constraint and conserved pattern in the genes coding for cox proteins than in other mitochondrial genes [119, 120].

Fig. 5

Ratio of non-synonymous to synonymous substitutions (ω) in the 13 PCGs of subterranean (coral color) and surface (cyan color) amphipods based on the free-ratio model. Boxes include 50% of values; ω is not significantly different between subterranean and surface amphipods for any gene except cox2^* (P value = 0.02)

To test if the 13 PCGs in subterranean lineages evolve at different relative rates compared to surface lineages, we compared a series of ML branch-based selection models (Table 3). For all PCG loci except atp8, nad3, and nad4l the saturated model (M2) where each branch had its own ω was favored. For atp8, nad3, and nad4l the best models were the M0 (single ω for all branches) and M1 (two ω model with one for surface and one for subterranean linages). In addition, the M1a model (six ω model with one for surface and one for each subterranean linage) was included in the set of best models for the atp8 and nad4l loci. To further test if specific branches have undergone variable selective pressures, especially those amphipod branches adapted to surface habitats, we employed the two-ratio branch model. When the ω values for each PCG were compared between each amphipod terminal branch and the other 34 amphipod taxa, several loci in surface amphipod mitogenomes were found to be undergoing positive selection (ω1 > ω0; Fig. 6; Supplementary Table S8). This suggests that many surface amphipods have experienced directional selection in their mitochondrial loci perhaps due to high energy demands and was in accordance to previous studies in other arthropods [19, 32, 121, 122]. In contrast, several loci in subterranean amphipod mitogenomes have undergone purifying selection (ω1 < ω0). Surprisingly, a few loci in subterranean taxa displayed positive selection (ω1 > ω0; Fig. 6; Supplementary Table S8). To test if individual gene codons were subject to positive selection, we implemented two pairs of site models (M1a vs. M2a and M8a vs. M8). The M8 model identified one positively selected site on the atp8 locus (37 N; p = 0) and one positively selected site on the nad5 locus (482 Q; p = 0). Similarly, The M2a model identified two positively selected sites (37 N & 31 S; p = 0.0194) on the atp8 locus (Table 4).

Table 3 AIC scores and ω estimates for various branch-based selection models for the 13 PCGs (one ω for all branches (M0), two-ratio model with background (surface) ω and single ω for subterranean branches (M1), two-ratio model with background (surface) ω and single ω for subterranean branches fixed at neutral evolution (ω = 1) (M1fixed), six-ratio model with background (surface) ω and a single ω for each subterranean lineage, B. brachycaudus (C1), Stygobromus clade (C2), B. jaraguensis (C3), Metacrangonyx clade (C4), P. daviui (C5) (M1a), and ω for each branch (M2)). The best-fit model(s) for each PCG is highlighted in red color

Table 4 Evidence of positive selection on the mitochondrial PCGs of subterranean and surface-dwelling amphipods based on site model

Fig. 6

Results of selective pressure analysis of mitochondrial PCGs with LRT P-value < 0.05 in subterranean and surface-dwelling lineages of amphipods based on branch 2 vs. 0 model. Different colored shapes represent different mitochondrial genes. Squares represent purifying selection and circles represent positive selection. Surface amphipod branches are colored blue and subterranean amphipod branches are colored red

Similar to flying grasshoppers that have evolved to adapt to increased energy demands to maintain the flight capacity [32], the mitochondrial loci of surface amphipods may have evolved mechanisms to meet increased energy demands due to predation, dispersal, and other factors. Although surface amphipods appear to be evolving under selective pressures different from those of the subterranean taxa and their mitochondrial loci have accumulated more nonsynonymous than synonymous mutations compared to subterranean taxa, the branch model tests did not clearly support positive selection on these branches, and we cannot rule out the influence of relaxed selection. Previous studies have demonstrated that positive selection will act on only a few sites for a short period of evolutionary time, and a signal of positive selection often is overwhelmed by continuous negative selection that sweeps across most sites in a gene sequence [123].

In contrast to branch models where ω varies only among branches, branch-site models allow selection to vary both among amino acid sites and lineages. Thus, branch-site models are considered quite useful in distinguishing positive selection from relaxed or purifying selection [123]. Using the more stringent branch-site model, we detected positive selection in 14 branches and 12 loci with a total of 308 amino acid sites under positive selection. Among them, 80 amino acid sites in seven loci (atp6, atp8, cox3, nad2, nad3, nad4, and nad5) were identified on the subterranean terminal branches, whereas 228 amino acid sites in 10 loci (atp6, atp8, cox1, cox2, cytb, nad1, nad2, nad3, nad5, and nad6) were identified on the surface terminal branches. Nearly three times as many positively selected amino acid sites were detected on surface branches compared to subterranean branches. Most of the positively selected loci on surface branches were found in C. forbesi with 114 sites (Fig. 7; Supplementary Table S9). In total, eight positive selected loci (atp6, atp8, cox1, cox2, cytb, nad1, nad4, and nad5) were identified by the branch-site model and by at least one other model on the surface branches, whereas only four positive selected genes (atp6, atp8, cox3, and nad5) were identified on the subterranean branches.

Fig. 7

Evidence of positive selection on the mitochondrial PCGs (LRT P<0.05) and positively selected site (BEB: P≥95%) in subterranean and surface-dwelling lineages of amphipods based on branch-site models. Different colored circles represent different mitochondrial loci. The number within each circle represents the number of positive selection sites detected for the locus. Surface amphipod branches are colored blue and subterranean amphipod branches are colored red

The identification of many positively selected amino acid sites suggests that episodic positive selection has acted on mitochondrial PCGs of surface amphipods. In addition, we also identified a few positively selected sites on subterranean branches primarily in B. brachycaudus with 39 sites and P. daviui with 25 sites (Supplementary Table S9). Bactrurus brachycaudus is usually associated with springs and caves [124], while P. daviui is associated with groundwater wells [77].

Direction and magnitude of selection pressures

Given the crucial role played by the mitochondrial genome in metabolic energy production [125], we hypothesized that the mitogenome of surface amphipods may show evidence of adaptation (directional selection) to life in surface habitats where energy demand is higher relative to subterranean habitats. We found support for directional selection in surface lineages based on three different selection analyses (RELAX, aBSREL, and BUSTED). In summary, all tests confirmed the existence of a moderate signal of positive or diversifying selection, as well as signal for significant relaxed purifying selection in the mitogenome of surface amphipods. This supports a previous study by Carlini and Fong [126] who reported evidence for relaxation of functional evolutionary constraints (positive or diversifying selection) in the transcriptome of a subterranean amphipod Gammarus minus.

We implemented aBSREL on the concatenated 13 PCG dataset comprising all 35 species as test branches and detected episodic diversifying selection in seven species: P. daviui (p = 0), O. nanseni (p = 0.0008), G. fasciatus (p = 0.0298), G. fossarum (p = 0.045), B. jaraguensis (p = 0.0016), C. forbesi (p = 0), and B. brachycaudus (p = 0.0001). We then used aBSREL to conduct independent tests for the crangonyctid species as the test branch and the remaining species as reference branches. We detected evidence of episodic diversifying selection in C. forbesi (p = 0) and B. brachycaudus (p = 0.0001) (Table 5). Using BUSTED, which provides a gene-wide test for positive selection, we detected evidence of episodic diversifying selection in three of the surface species: C. forbesi (p = 0.011), G. fasciatus (p = 0.033), G. antarctica (p = 0.009), whereas evidence of gene-wide episodic diversifying selection was found in just one of the subterranean species, P. daviui (p = 0.020) (Table 5). Using RELAX, which tests whether the strength of selection has been relaxed or intensified along a specified set of test branches, we detected selection evidence of relaxed selection in C. forbesi (p = 0) and other surface species, including O. nanseni, G. fasciatus, G. fossarum, G. antarctica, and P. kessleri. Contrastingly, evidence of intensification of selection was detected in subterranean species, including S. tenuis (p = 0), S. allegheniensis (p = 0.0025), S. indentatus (p = 0), and S. pizzinii (p = 0). Surprisingly, a few of the surface species including C. mutica (p = 0.015), E. cyaneus (p = 0), and P. japonica (p = 0) exhibited intensification of selection and subterranean species including P. daviui (p = 0) and M. dominicanus (p = 0.015) exhibited relaxation of selection (Table 5).

Table 5 Selection signals in the mitogenomes of amphipods inferred using aBSREL, BUSTED, and RELAX algorithms. The dataset comprising all 13 concatenated protein-coding genes with 3,607 amino acid sites in the alignment. K column: a statistically significant K > 1 indicates that selection strength has been intensified, and K < 1 indicates that selection strength has been relaxed. LR is likelihood ratio and D indicates the direction of selection pressure change: intensified (I) or relaxed (R), where * highlights a statistically significant (p < 0.05) result. Mitogenomes with significant LRT P -value < 0.05 are highlighted in red color

In addition to the concatenated 13 PCG dataset, we also conducted selection analyses for each PCG to determine which genes might be evolving under unique selection pressures. We found evidence of directional selection in atp8 of C. forbesi (p = 0.026) and nad3 of S. pizzinii (p = 0.041) using aBSREL and cox3 of B. brachycaudus (p = 0.029) using BUSTED. Atp8 of the surface amphipod C. forbesi exhibited strong evidence of directional selection, which was quite surprising as atp8 is a small locus sometimes missing from metazoan mitogenomes and normally evolves under highly relaxed selection pressures [127]. RELAX analyses uncovered five loci (cox1, cox3, cytb, nad1, and nad3) that exhibited relaxed selection and one gene (atp6) that exhibited intensification of selection in C. forbesi. Similarly, three loci (cox3, nad5, and nad6) in B. brachycaudus showed evidence of relaxed selection. Several loci in other subterranean species, including S. tenuis, S. allegheniensis, and S. pizzinii, exhibited varying levels of intensification of selection, whereas none exhibited relaxed selection (Table 6). Some of these outliers were expected, as nad5 and nad6 are known to evolve faster among the mitochondrial loci [128]. Also, evidence for relaxation of functional evolutionary constraints (positive or diversifying selection) has been reported in the nad family of subterranean Gammarus species adapted to the subterranean environment [126]. Although this may explain outliers in the subterranean B. brachycaudus mitogenome, it remains unclear why cox3 exhibited signatures of relaxed selection. This gene is generally one of the most conserved mitochondrial loci in animals [92, 129, 130], and high levels of purifying selection has been observed in the cox family in other amphipod species [29]. In C. forbesi, atp6 showed signatures of positive selection, which contrasted most other PCGs in its mitogenome that exhibited relaxed selection. Overall, in accordance with the results obtained using the concatenated dataset, individual mitochondrial loci of subterranean amphipods mostly exhibited varying levels of purifying selection, whereas surface amphipods predominantly exhibited more relaxed selection.

Table 6 Selection signals in the mitochondrial PCGs of crangonyctid amphipods sequenced in this study inferred using aBSREL, BUSTED, and RELAX algorithms. K column: a statistically significant K > 1 indicates that selection strength has been intensified, and K < 1 indicates that selection strength has been relaxed. LR is likelihood ratio and D indicates the direction of selection pressure change: intensified (I) or relaxed (R), where * highlights a statistically significant (p < 0.05) result. PCGs with significant LRT P -value < 0.05 are highlighted in red color

To provide further evidence of positive selection, we implemented the RELAX, aBSREL, and BUSTED algorithms on the branch, branch-site, and site models. Eight loci (atp8, cox1, cox2, cytb, nad1, nad4, nad5, and nad6) all involved in the OXPHOS pathway were under positive selection in surface branches by at least two methods. The loci nad1, nad4, nad5, and nad6 encode the subunits of NADH dehydrogenase, also called Complex I, that initiates the oxidative phosphorylation process. Complex I is the largest and most complicated proton pump of the respiratory chain and is involved in electron transfer from NADH to ubiquinone to supply the proton motive force used for ATP synthesis [131], Complex I plays a key role in cellular energy metabolism by pumping gradient of protons across the mitochondrial membrane producing more than one-third of mitochondrial energy [132]. Genes cox1 & cox2 encode the catalytic core of Cytochrome c oxidase also called Complex IV. Complex IV is directly involved in electron transfer and proton translocation [133]. Gene atp8 encodes a part of ATP synthase, also called Complex V, and plays a major role in the final assembly of ATPase [133]. In summary, our selection analyses revealed signals of positive selection in several mitochondrial genes of surface amphipods, which may be associated with increased energy demands in surface environments. In contrast, subterranean amphipods showed signatures of purifying selection, which may be related to maintaining efficient energy metabolism in subterranean habitats.

Conclusion

In this study, we compared mitogenome features including AT/GC-skew, codon usage, gene order, phylogenetic relationships, and selection pressures acting upon amphipods inhabiting surface and subterranean habitats. We described a novel mitochondrial gene order for C. forbesi. We identified a signal of directional selection in the protein-encoding genes of the OXPHOS pathway in the mitogenomes of surface amphipods and a signal of purifying selection in subterranean species, which is consistent with the hypothesis that the mitogenome of surface-adapted amphipods has evolved in response to a more energy demanding environment compared to subterranean species. Our comparative analyses of gene order, locations of non-coding regions, and base-substitution rates points to habitat as an important factor influencing the evolution of amphipod mitogenomes. However, the generation and study of mitogenomes from additional amphipod taxa, including other crangonyctid species, are needed to better elucidate phylogenetic relationships and the evolution of mitogenomes of amphipods, as mitogenomes are available for just a tiny fraction of the more than 10,000 described amphipods. In addition, more evidence is needed to further validate our inferences, such as studying the effects of amino acid changes on three-dimensional protein structure and function. Nevertheless, our study provides a necessary foundation for the study of mitogenome evolution in amphipods and other crustaceans.