In December 2018, we visited intertidal zones in Sanya City, Hainan Province, China. During low tides, we carefully searched coral crevices in the intertidal zone and collected 13 specimens of Desis jiaxiangi (see Additional file 2: Fig. S1) in the location (N18.21999°, E109.51128°, 3 m a.s.l.). Voucher specimens were deposited in the Marine Bank of the China National GeneBank (CNGB voucher nos. Des_001 to Des_013), Shenzhen, China.
Assembly and annotation of the complete mitogenome of Desis jiaxiangi
We extracted total DNA from the whole body of one specimen (Des_002) for traditional whole-genome shotgun sequencing using a Puregene Tissue Core Kit A (Qiagen, Germantown, MD, USA). We constructed the library with a short-insert size of 500 bp using standard protocols (Illumina, San Diego, CA, USA). We sequenced with paired-end reads (150 bp in length) on an Illumina Hiseq X Ten platform at BGI-Wuhan, China. The raw sequencing data was filtered by SOAPnuke  to trim adapter, low quality, high N base ratio and etc. Finally, we acquired about 100 Gb of clean reads after the sequencing and data filtering. We used all of these clean data as an input file to the Python3-based mitogenome assembly toolkit MitoZ  for a primary mitogenome assembly of Desis jiaxiangi.
Subsequently, we performed a BLAST search to compare this mitogenome assembly with corresponding mitogenome data of four relatively close spider species from RTA clade: Agelena silvatica (GenBank Accession no. KX290739.1), Argyroneta aquatica (No. NC_026863.1), Dolomedes angustivirgatus (No. NC_031355.1), and Telamonia vlijmi (No. KJ598073.1). We manually annotated the conserved regions of our assembled Desis mitogenome based on the other four mitogenomes. For the uncertain regions of the assembled Desis mitogenome, we designed two pairs of primers (see Additional file 1: Table S2) using Primer Premier 5 (Premier Biosoft, Palo Alto, CA, USA) to obtain the missing sequences using PCR. To validate the sequences of the uncertain regions, we extracted genomic DNA from the leg tissue of the specimens using a Puregene Tissue Core Kit A (Qiagen).
We performed PCR in 50 µL volume tubes, each containing 25 µL 2 × Taq PCR MasterMix (Tiangen Biotech, Beijing, China), 2 µL of each primer (10 µM), and 2 µL genomic DNA (100 ng/µL), with 19 µL double-distilled water. The PCR protocol was as follows: pre-denaturation at 94 °C for 5 min, 35 cycles of denaturation at 94 °C for 30 s, annealing at 48 °C for 30 s, elongation at 72 °C for 1 min 40 s, and a final elongation at 72 °C for 10 min. A Veriti Thermal Cycler (Applied Biosystems, Carlsbad, CA, USA) was used for the PCR. We checked the PCR products by electrophoresis using 1% agarose gels.
After successfully amplifying the uncertain regions, we obtained the final assembly of the complete mitogenome of Desis jiaxiangi, which was annotated with MITOS2  on MITOS WebServer (http://mitos2.bioinf.uni-leipzig.de/index.py). The final mitogenome sequence was visualized with MitoZ in the visualize module .
Comparative analysis of mitochondrial gene orders is a powerful method of revealing ancient events in the process of species evolution . After completely annotating the mitogenome of Desis jiaxiangi, we compared its gene order with the available 45 complete spider mitogenomes (see Additional file 1: Table S1) from the NCBI GenBank.
Phylogenetic analyses and estimation of divergence time
We prepared a data set comprised the novel mitogenome of Desis jiaxiangi and the complete mitogenomes of 45 other spider species (see detailed species names in Additional file 1: Table S1) downloaded from the NCBI GenBank. We performed phylogenetic analyses of this data set using Bayesian inference (BI) and maximum likelihood (ML) methods. We aligned the nucleotide and amino acid sequences of the 13 PCGs from each mitogenome using the multiple sequence alignment program Clustal Omega in EMBL-EBI . We performed a best-fit model selection of amino acid replacement using ProtTest-3.4.2 . Based on Akaike’s information criterion, MtREV + I + G was chosen as the best model for both inference methods. We reconstructed the phylogenetic tree with the ML method using PhyML 3.0 . The node support values were estimated with 100 replicates and other parameters as the default. We also performed BI using MrBayes v3.2.6  to compare the topologies of the ML phylogenetic trees. MCMC algorithm parameters were set for two independent runs with four chains (one cold chain and three heated chains) for 10,000,000 cycles. The sample frequency parameter was set at 1000 for sampling each chain every 1000 cycles. The first 25% of the runs were discarded as burn-in. We used FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/) to visualize the derived BI and ML trees.
We used the Bayesian MCMCTree program in PAML package v4.9j  to estimate divergence times. Optimized parameters were as follows: clock = independent rates, model = HKY85, nsample = 20,000, burn-in = 2,000. We used calibration points from recommended fossils and a related time in a recent spider fossil review . Given the limited spider species with available mitogenomes, we used three fossils to calibrate the phylogenetic tree: Palaeothele montceauensis (299–304 Ma) for the Mesothelae stem, Eoplectreurys gertschi (164–175.1 Ma) for the Synspermiata stem, and Montsecarachne amicorum (125–129.4 Ma) for the Synspermiata crown. As the fossil of Almolinus ligula (43–47.8 Ma) in the Salticidae crown can be placed in the superfamily Hisponinae , and no sample mitogenome is available for this superfamily, this fossil was not used for calibration in the present work.
To evaluate potential adaptive evolution in the mitochondrial genes of intertidal Desis spiders, we performed positive selection analyses using PAML4.9j . Synonymous substitutions in protein-coding sequences cannot cause changes in amino acids, which are typically found in the third, or sometimes first, position of a codon . We thus used a gene-level approach based on the ratio of nonsynonymous (Ka) to synonymous (Ks) substitution rates to detect potential positive selection signal on PCGs across closely related or divergent species [52, 53]. A Ka/Ks ratio of 1, <1, or> 1 in protein-coding sequences may be interpreted as a neutral mutation, a negative (purifying) selection, or a positive (diversifying) selection, respectively . To investigate the variation in nucleotide substitution rates in spider mitogenomes, we retrieved all 13 mitochondrial PCGs (nd2, cox1, cox2, atp8, atp6, cox3, nd3, nd5, nd4, nd4l, nd6, cytb, and nd1) from each annotated mitogenome of the 46 spider species examined. We aligned each gene separately with the codon-based model in the Muscle module of MEGA7 . Ambiguous regions in each alignment were removed with Gblocks v0.91b . Ka, Ks, and the Ka/Ks ratio across all 13 PCGs were calculated with KaKs_Calculator v2.0 .
To examine positive selection pressure on individual PCGs of the 46 spider species, we generated conserved blocks from codon alignments of each PCG using Gblocks v0.91b. We ran both the branch and branch-site models in the CodeML program of PAML on those codons from the conserved blocks. We used the phylogenetic topology inferred from PhyML as the guide tree for PAML analyses. In branch models, the free Ka/Ks ratio model is allowed to vary among branches to detect positive selection on the foreground branch. In branch-site models, the one-ratio model assumes that all branches have an identical Ka/Ks ratio, whereas the two-ratio model assumes that the foreground branch has a different Ka/Ks ratio from the background branches. The Ka/Ks ratios in branch-site models are thus allowed to vary both among sites and across branches to detect positive selection on a few sites along the foreground branch. The branch-site model A null fixed all Ka/Ks ratios to 1, whereas the branch-site model A (positive selection model) did not fix the Ka/Ks ratio. The branch-site model A was used to detect positive selection sites along the lineages of aquatic spiders (i.e., the foreground branch). The presence of a site with Ka/Ks ratio > 1 is suggested when the positive selection model A fits the data significantly better than the corresponding model A null.
We first used the free-ratios model (branch model) to calculate the average Ka/Ks ratio for all 13 PCGs in each lineage to represent their evolutionary rate of mitogenomes. We paid particular attention to two aquatic species, the intertidal spider (Desis jiaxiangi) and the freshwater spider (Argyroneta aquatica), which have aquatic habitats, whereas the other species are almost all limited to land. For this purpose, we marked the ancestral branch containing both Desis jiaxiangi and Argyroneta aquatica as the foreground branch and the rest of the species as the background branches. To compare the results, we marked the branches of Pirata subpiraticus, Hypochilus thorelli, Argyroneta aquatica, Desis jiaxiangi and Agelena silvatica separately to test these five species independently. We ran the branch models (one-ratio model vs. two-ratio model) and branch-site models (null model A vs. model A) as two pair models using likelihood ratio tests and chi-square tests, respectively, to detect whether the positive selection signals were significant. The null hypothesis assumed that all branches had a common Ka/Ks ratio. An alternative hypothesis assumed that the ancestral branch had an independent Ka/Ks ratio that differed significantly from the background branches with a common Ka/Ks ratio. If this hypothesis has a better fit than the null hypothesis, we considered the occurrence of positive selection in the ancestor of aquatic spiders in certain gene(s) when p < 0.05. Positive selection does not usually generate the function on the whole length of the target genes, and only a few sites can reflect the positive selection signal. The branch-site models further identified the positive site(s), and the Bayes empirical Bayes analysis was used to calculate posterior probabilities to detect those sites under positive selection.