Reinvestigating the phylogeny of Myriapoda with more extensive taxon sampling and novel genetic perspective

Taxon sampling

Building upon previous works by Fernández, Edgecombe & Giribet (2018) and Szucsich et al. (2020), where more than 36 species representing the four major groups of myriapods were included in taxon sampling, we sequenced three additional species (one chilopod: Scolopendra sp.; and two diplopods: Epanerchodus sp., Skleroprotopus sp.) in this study. Our sampling was designed to maximise the representation of myriapod groups. Information on sampling localities and accession numbers in the Sequence Read Archive (SRA) database for each transcriptome is shown in Table 1, including four genomes from http://metazoa.ensembl.org. Twenty-one outgroups were also included: eight chelicerates (Liphistius malayanus, Centruroides vittatus, Damon diadema, Archegozetes longisetosus, Araneus diadematus, Egaenus convexus, Euscorpius sicanus, and Nymphon gracile), two onychophorans (Peripatopsis capensis and Peripatoides novaezealandiae), and 11 pancrustaceans (Daphnia pulex, Folsomia candida, Drosophila melanogaster, Eubranchipus grubii, Triops cancriformis, Nebalia bipes, Anaspides tasmaniae, Hemidiaptomus amblyodon, Tisbe furcata, Vargula hilgendorfii, and Xibalbanus tulumensis).

RNA extraction and sequencing

Following the manufacturer’s instructions, total RNA was extracted using a commercial RNA extraction kit (TAKARA). Samples were treated with Ambion turbo DNA-free DNase to remove residual genomic and rRNA contaminants during mRNA purification. The quantity and quality (purity and integrity) of mRNA were assessed using a NanoDrop ND-2000 UV spectrophotometer (Thermo Fisher Scientific).

For mRNA sequencing library preparation, mRNA was first enriched and purified with oligo (dT)-rich magnetic beads and then broken into short fragments, followed by paired-end sequencing on an Illumina Hiseq 4000 platform.

Data processing and de novo assembly

Sequencing adaptors and low-quality sequences were trimmed using Trimmomatic v 0.36 (Bolger, Lohse & Usadel, 2014) with default parameters. The clean data were assembled with Trinity (release 2.11.0) with 100 GB memory and a path reinforcement distance of 50 (Grabherr et al., 2011). The redundancy of all assembled transcripts was removed using CD-HIT v. 4.8.1 (under the cd-hit-est mode, with default parameters Li & Godzik, 2006). TransDecoder v5.5.0 (https://github.com/TransDecoder/TransDecoder) was then utilised for nucleotide sequence translation and longest ORF selection.

Orthology assignment and phylogenetic matrix construction

Both classical pipeline (OrthoFinder) and single-copy gene selection method (Benchmarking Universal Single-Copy Orthologs; BUSCO) were used in orthology assignment among the 60 selected taxa. Orthogroups were first identified using OrthoFinder v 2.2.7 (Emms & Kelly, 2015) with default settings (BLASTP E value ≤ 1e−5 and MCL inflation parameter of 1.5). Single-copy genes in arthropods were identified in our datasets with BUSCO v4.1.4 with default settings (E value ≤ 1e−6) based on hidden Markov model profiles, where BUSCO dataset arthropoda_odb10 (https://busco-data.ezlab.org/v4/data/lineages) was used as reference (Seppey, Manni & Zdobnov, 2019). For each BUSCO and each taxon, the longest hit of duplicated BUSCO homologous genes was retained for further analysis.

Putative orthogroup filtering was based on the gene occupancy threshold, which means that an orthogroup (or BUSCO) was selected if it could be found in more than or equal to the threshold number of taxa. For example, a 50% gene occupancy threshold would select orthogroups that were present in ≥ 50% of the included taxa. We selected two thresholds of 100% and 90% gene occupancy to obtain information on most species and to minimise the computational burden. Protein sequences of each orthogroup were aligned with MAFFT v 7.305b (maxiterate was set to 1000 in globalpair mode) prior to concatenation (Katoh & Standley, 2013). Then, Aliscore v2.0 was used to perform each orthogroup’s multiple sequence alignment (MSA) for ambiguous or randomly aligned sections’ identification, followed by Alicut v2.2 for error section’s trimming (Kück, 2009). After filtering with 100% and 90% gene occupancy thresholds, two raw matrices were constructed using custom Python scripts. Finally, the matrices were trimmed using MARE (v 0.1.2) to select optimised data subsets from the supermatrices for phylogenetic inference (Meyer, Meusemann & Misof, 2011).

Best partition schemes finding

The best-fit partitioning schemes and models of evolution for phylogenetic analyses were searched and estimated with PartitionFinderV2.0.0 (Lanfear et al., 2017), where amino acid substitution models were restricted to LG, WAG, JTT, and BLOSUM62. The corrected Akaike Information Criterion (AICc) model was selected and was set under greedy search.

Phylogenetic tree inference

Maximum-likelihood tree

Tree searches were performed for the two supermatrices with a ML approach, using IQ-TREE (v1.6.12) with the above best partitioning schemes. Statistical support was derived from 100 non-parametric slow bootstrap replicates. Replicates to perform SH-like approximate likelihood ratio test were set to 1000 and unsuccessful iterations to stop were set to 300. The initial tree searches were set from a completely random tree (Nguyen et al., 2015). The full command was: ‘iqtree -s MSAmatrix.phy -alrt 1000 -b 100 -t RANDOM -spp partition_ file. nex-nstop 300.’

Bootstrap support inference

Bootstrapping analyses were applied with RAxML-NG v1.0.0 (Kozlov et al., 2019), and the autoMRE bootstrap convergence test was set for a sufficient number of replicates. The bootstrap support was then mapped onto the phylogenetic trees using RAxML v8.2.11 (Stamatakis, 2014).

Topology testing

To evaluate support for the different hypotheses concerning the relationship proposed by previous studies among the four major classes within Myriapoda, topology tests were run for each dataset using the corresponding best partition scheme with IQ-TREE (v1.6.12), where the AU-test, weighted KH-test, and weighted SH-test were included, and all tests performed 100,000 resamplings using the resampling of estimated log-likelihoods (RELL) method (Nguyen et al., 2015). Six proposed hypotheses representing the two most controversial phylogenetic relationships (Edafopoda and Dignatha) of four major classes in Myripoda were compared in our topology tests; detailed information is shown in Fig. 1.

Hypothesis Eda.1(topology: ((CHI,DIP), (SYM, PAU));).

Hypothesis Eda.2(topology: (CHI, (DIP, (SYM, PAU)));).

Hypothesis Eda.3(topology: (DIP, (CHI, (SYM, PAU)));).

Hypothesis Dig.1(topology: ((CHI, SYM), (DIP, PAU));).

Hypothesis Dig.2(topology: (CHI, (SYM, (DIP, PAU)));).

Hypothesis Dig.3(topology: (SYM, (CHI, (DIP, PAU)));).

Identification of LIT genes

A modified version of the phylogenetically informed annotation tool pipeline named PIA2 (https://github.com/xibalbanus/PIA2) was applied for the identification of visual opsins in this study, and the parameters were set default (Pérez-Moreno et al., 2018). In this pipeline, 111 genes from the LIT, a collection of genes that underlie the function or development of light-interacting structures in metazoans, representing 13 different parts in visual pathways (photoreceptor specification, retinal determination network, phototransduction, rhabdomeric, phototransduction, ciliary, retinoid pathway, vertebrate, retinoid pathway, invertebrate, melanin synthesis, pterin synthesis, ommochrome synthesis, heme synthesis, crystallins, diurnal clock, and opsin), was taken as reference (Speiser et al., 2014). We applied the pipeline to protein sequences of each species.

Identification of positively selected genes

Evidence of positive selection was indicated by estimating the ratios of nonsynonymous substitutions (Ka or dN) and synonymous substitutions (Ks or dS), also called substitution rates (Ka/Ks or dN/dS value). The coding sequence of each identified LIT genes was aligned between a pair of taxa separately with MAFFT v 7.305b with default settings. And then the substitution rate was calculated using the KaKs_calculator with the following settings, method of calculation: GMYN, genetic code table: The Echinoderm and Flatworm Mitochondrial Code (Wang et al., 2010).

Transcriptome assembly and phylogenomic dataset construction

NGS technologies have empowered phylogenomic analyses in the last few decades. It has dramatically increased the size of datasets applied to phylogenetic questions. Within the framework of combining NGS technologies and phylogenomic techniques, we decided to re-investigate Myriapoda phylogeny with three newly sequenced species. Combined with published data from two outstanding studies on Myriapoda phylogeny, the data from a total of 60 species (39 from Myriapoda, 8 Chelicerata, 11 Crustacea, and 2 Onychophora) were used in this study (Table 1) (Fernández, Edgecombe & Giribet, 2018; Szucsich et al., 2020). Except for the four species with published genome data, raw reads of the remaining 56 species were trimmed and assembled de novo. Orthology assignments of the 60 species were mainly based on BUSCO results, and more than 70% of the ortholog gene set (BUSCO dataset: arthropoda_odb10, comprising 1013 single-copy protein-coding genes or ortholog group, OG) were identified in 56 species (details in Table S1). Additionally, 786 of the BUSCO orthology assignments were confirmed using OrthoFinder (details in Table S1).

Concatenated supermatrices were compiled using a threshold of percentage gene occupancy of 100% and 90% (Fig. 2) (González et al., 2015). We found that 32 OGs were represented in all 60 species (100% gene occupancy), of which 20 were confirmed by OrthoFinder. There were 505 OGs represented in more than 54 species (90% gene occupancy), and 369 OG assignments were confirmed by OrthoFinder. Thus, two datasets comprising 20 and 369 OGs were obtained (Table S1). After MSA, identification, and removal of ambiguously aligned sections in each dataset, two phylogenomic supermatrices on the amino acid levels were constructed, which are hereafter referred to as OCC100 and OCC90. Matrix OCC100 included 20 gene partitions and spanned 8,401 aligned sites with 0.395 overall information content. Matrix OCC90 included 369 gene partitions and spanned 129,085 aligned sites with 0.318 overall information content.

Phylogenetic tree inference and topology analysis

We constructed Maximum-Likelihood (ML) trees based on the best partition schemes and best-fitting substitution models schemes with matrices OCC100 and OCC90. Three main results were found from the inferred trees. All the analyses recovered Myriapoda as the monophyletic sister group of Pancrustacea with high support (Fig. 3). As for the relationships among the four myriapod classes: Symphyla (SYM), Chilopoda (CHI), Diplopoda (DIP), and Pauropoda (PAU), we found a sister group relationship of CHI+DIP, and another sister group relationship of PAU+SYM (Fig. 3). Both were highly supported by the bootstrap result (PAU+SYM: 100%, DIP+CHI: 100%). The three newly sequenced species (Scolopendra sp., Epanerchodus sp., and Skleroprotopus sp.) were positioned in the expected clades.

A variety of groupings of the Myriapoda classes have been proposed, where two hypotheses, Edafopoda and Dignatha, received the most attention (Fig. 1). Edafopoda is a grouping of PAU+SYM which has been supported by shared genetic sequences (Fig. 1). However, in Dignatha, the PAUs were positioned with the DIPs. In this study, all trees inferred were congruent with the unrooted quartet topology with CHI+DIP and PAU+SYM (Hypothesis Eda.1, Fig. 1). We conducted three types of topology tests—AU test, KH test, and weighted SH test—on the quartet topology of Edafopoda and Dignatha, where four different phylogenomic datasets were applied. The results consistently supported the topology Hypothesis Eda.1, which is the only topology to not be rejected in any test (Table 2). Almost all hypotheses derived from Dignatha were rejected with high significance, especially in the phylogenomic matrix OCC90 (Table 2). When comparing the results from the two phylogenomic matrices, we found that all testing results of matrix OCC100 were consistent with and covered by that of matrix OCC90. Under matrix OCC90, the datasets that were different in outgroup selection (PAN or CHE) exclusively showed divergence when rejecting Hypothesis Eda.2, where all three topology tests on the datasets with CHE were not rejected, which was completely opposite to the results of datasets with PAN (Table 2). In other words, we found that the sister group of Edafopoda (PAU+SYM) received less support from CHI than the clade of DIP+CHI, but it cannot be completely denied.

Outgroup dependence of myriapod phylogeny inference

Despite the quartet topology of CHI+DIP and PAU+SYM being recovered in our analyses, the relationships among the four major classes in Myriapoda varied across phylogenomic datasets, with dependence on outgroup selection proposed in previous studies (Fernández, Edgecombe & Giribet, 2018; Szucsich et al., 2020). Given that the phylogeny inference was sensitive to outgroup choice, we conducted topology tests on datasets with different clades of outgroups: one with only CHE as outgroup, and the other with only PAN (Figs. S2–S5). ML tree inference of the former resulted in a sister group relationship of CHI+DIP and another sister group relationship of PAU+SYM, which was congruent with the results inferred from the datasets with the full taxon sampling (Fig. 3). However, we found that the ML tree inferred from the latter datasets resulted in a sister group relationship of CHI and DIP, with SYM as a sister to this clade, followed by PAU. Although these quartet topology results were also obtained in previous studies, negligible support could be obtained from bootstrapping analyses (Szucsich et al., 2020).

LIT genes’ identification in Myriapoda

Using the PIA2 pipeline, we identified 2,001 transcripts in 39 Myriapoda species as putative components involved in the development of light-interacting structures, including 96 LIT genes, which are important components of 11 visual pathways (Fig. 4). A total of 13 visual pathways were compiled in the pipeline, and two were absent in this study (Table 3, retinal determination network and opsin synthesis) which involved 10 LIT genes. In addition, the other five absent LIT genes were Gq_gamma, RBP3, Dat, TYR, and reflectin_1a. We investigated the completeness of the amino acid level of each ortholog by calculating the ratio of the length of the identified peptide to the target reference peptide, as depicted in Fig. 4; the lighter the cell-filling colour, the more incomplete the transcript. As shown in Fig. 4, the distribution patterns of the identified LIT genes among the four classes are very similar; LIT genes from prc (photoreceptor specification), reti (retinoid pathway, invertebrate), heme (heme synthesis), and crys (crystallins) were rarely identified in Myriapoda (shown as a large blank area in Fig. 4). We also found that LIT genes from reti (retinoid pathway, invertebrate, 0), ommo (ommochrome synthesis, 2), and clock (diurnal clock, 2) were rarely identified in PAU (Table 3). As for the common LIT genes among these four classes, we found that 62 LIT genes could be identified in at least one of the four classes (Fig. 5C), and three (GC, TH, KF) could be identified in all 39 myriapods investigated. The following three genes were separately involved in three different pathways, ctrans: phototransduction in ciliary, mel: melanin synthesis, and ommo: ommochrome synthesis. We compared the LIT genes co-identified between DIP and CHI, and PAU and SYM according to the visual pathways in which the genes participated (Fig. 5B and Fig. 5A). We found that clock (diurnal clock), crys (crystallins), ommo (ommochrome synthesis), and reti (retinoid pathway, invertebrate) were more abundant in the sister group of DIP and CHI.

Selection tests on LIT genes

To test whether genes associated with the evolution of light interactions in Myriapoda have undergone potentially adaptive changes, Ka/Ks calculations were conducted. A total of 23,832 aligned LIT gene pairs were calculated, including 8,717 pairs from CHI&DIP, 2,048 from CHI&SYM, 1,220 from CHI&PAU, 1,158 from DIP&PAU, 1,981 from DIP&SYM, 278 from PAU&SYM, 4,241 from CHI&CHI, 3,981 from DIP&DIP, 47 from PAU&PAU, and 161 from SYM&SYM. Positive selection was detected in 27 pairs, indicated by Ka/Ks ≥1, five pairs from CHI&DIP, one from CHI&PAU, three from CHI&SYM, two from DIP&PAU, two from DIP&SYM, seven from CHI&CHI, and seven from DIP&DIP. As depicted in Fig. 5, no evidence of positive selection for LIT genes was found in PAU&PAU, SYM&SYM, PAU&SYM. Values of Ka/Ks in the range of 0.5 to 1.0, which indicates relaxed selection, were observed in 395 pairs, covering all classes combinations (details in Table S2). The remaining 23,435 pairs had Ka/Ks values ranging from 0.0002 to 0.5, representing 98% of the pairs we calculated, which means that most of the genes in the four major classes were under purifying selection. Positive selection was detected in the following 14 LIT genes: clot, Cnga1, CSAD, DAGK, DDC, Galpha_it, GC, Gprk1, Gq_alpha, Pde6abc, PKC, PLC, RBP1, RDH8, timeless, and trp, which cover the clock, ctrans, mel, pter, retv, and rtrans visual pathways (Fig. 6). In addition, the transient receptor potential protein trp, which encodes a component of the rhabdomeric phototransduction pathway, was identified and positively selected in PAU&CHI, SYM&CHI, and SYM&DIP.