This study was carried out in accordance with the recommendations of the Guide for the Care and Use of Laboratory Invertebrate Animals. The protocol was approved by the Ethical Committee of Researches of the Nanjing University (NJU).

Adult species of B. japonicum were collected from Qingdao, Shandong province, China, and kept in laboratory for 5 days to facilitate emptying of the digestive tract. Total RNA was isolated from several male and females using the TRIzol reagent (Invitrogen, United States), and genomic DNA was removed by RNase-Free DNase Set (Qiagen, Germany) according to the manufacturer’s instructions. RNA quality was determined using Bioanalyser 2100 (Agilent, United States), and the concentration was measured using the NanoDrop 1000 (Thermo Scientific, United States). Standard cDNA libraries were constructed using standard Illumina libraries prep protocols and TruSeq kit (Illumina, United States). After homogenization treatment of sequencing libraries, RNA sequencing (RNA-seq) was conducted on an Illumina HiSeq^TM 2000 with 100-bp paired-end reads. RNA-seq was performed in BGI (Shenzhen, China).

We firstly used FastQC (Brown et al., 2017) to control the quality of the raw sequencing reads by checking for over-represented and potentially contaminant sequences following stringent criteria: (1) reads with adaptors were discarded; (2) reads with unknown bases > 10% were discarded; (3) reads with a length < 20 bp were discarded. We used the program Sickle to remove or trim low-quality reads [percentage of low-quality bases (bases with sequencing quality score ≤ 5) > 50%] (Cabili et al., 2011). Quality paired-end reads (clean reads) were obtained and then used to assemble the transcriptome using Trinity (Grabherr et al., 2011). In addition, we obtained clean reads of a further adult amphioxus (A. lucayanum) (Accession: SRX437621) from the Sequence Read Archive (SRA) database of the NCBI, and assembled its transcriptome following the above pipelines. Trinity can detect potential isoforms from alternative splicing and label them with the same prefix. For multiple isoforms, the longest unigene was selected from each isoform group as a unique representation for that group. We also downloaded the B. lanceolatum transcriptome from the NCBI Transcriptome Shotgun Assembly (TSA) database (Accession numbers: JT846176-JT905674) (Silvan et al., 2012). Potential coding sequences (CDSs) of three amphioxus species were predicted via Blastx (default parameter) (Mount, 2007) search of the NCBI non-redundant (NR) protein database. The optimal alignment was extracted as template for the determination of the CDSs of unigenes. Then, ESTScan (E-value = 10^-5) (Iseli et al., 1999) software was used to predict the CDS of unigenes that failed to produce any hit.

The protein, mRNA, and CDS sequences of 16 sampled species with sequenced genome were retrieved from online databases (see Supplementary Table S1 for details). Among these genomic sequences, the CDS and protein sequence of 15 species (except for A. mississippiensis) could not be completely matched, thus they were corrected. First, CDS sequences below 120 bp were removed; then, protein sequences were aligned to mRNA using Exonerate software (Guy and Birney, 2005). Based on amino acid-nucleic acid alignments, incompatible CDS and protein sequences were corrected, discarding uncorrected sequences. For the protein sequences translated from corrected CDSs, we used the Blastp tool to search for them in Nr, Swiss-Prot, and gene ontology (GO) databases (Table 1), respectively, performing functional annotation of all genes or unigenes.

HaMSTR (Ebersberger et al., 2009) and reciprocal best hits (BRH) (Moreno-Hagelsieb and Latimer, 2008) methods are generally used for the construction of orthologous genes. The results obtained from HaMSTR and BRH methods were intersected as a final ortholog set of 19 sampled species, to avoid bias caused by using a single method only. For the HaMSTR approach, according to previous description in transcriptomic analyses (Misof et al., 2014; Li et al., 2015), a set of core-orthologs was constructed from nine vertebrate genomes including human (Homo sapiens), mouse (Mus musculus), chicken (Gallus gallus), zebrafish (Danio rerio), Coelacanth (Latimeria chalumnae), fugu (Takifugu rubripes), elephant shark (Callorhinchus milii), lamprey (Petromyzon marinus), and Western clawed frog (Xenopus tropicalis). All 6,793 one-to-one core-orthologous proteins were obtained from the OrthDB v9.1 database (Kriventseva et al., 2015) and were inputted into the multiple alignment tool Muscle (v3.8) (Edgar, 2004) using default sets. Then, a hidden Markov model “primer taxa” was built from the multiple alignment of the core-orthologs using the hmmbuild tool of the HMMER3 package (Finn et al., 2011). The “primer taxa” served as input to produce the core-ortholog database for the program HaMStR v.13.2.6 (Ebersberger et al., 2009) to search for the corresponding orthologs in the remaining 10 species. If one species contained multiple corresponding orthologs in co-orthologous, only the optimal hit was retained. HaMStR v.13.2.3 was run using strict parameters (-sequence file, -est, -hmmset, -refspec, -representative, and -ublast) according to a previous description (Li et al., 2015). Furthermore, Proteinortho (v5.13) (E-value = 10^-10) (Lechner et al., 2011), software compiled based on the BRH method that is widely used for large-scale comparative genomic analysis, was employed to identify co-orthologous among all 19 species. The PhyloTreePruner pipeline was employed to pick the unique representative for each species in co-orthologous based on the gene tree (Kocot et al., 2013). According to the gene labels, the obtained orthologous gene groups were collectively included in the results from HaMSTR and BRH methods. All orthologous genes, including protein and CDS sequences, were aligned using the Prank tool (“-codon” option) (Löytynoja and Goldman, 2005), and then further trimmed via Gblocks (Talavera and Castresana, 2007) with the parameter “-t = c” to remove poorly aligned regions. Trimmed alignments that contained shorter than 60 bp/20 amino acids were removed.

Concatenated alignments used for the construction of the phylogenetic tree were constructed from all orthologs using the FasParser package (Sun, 2017). The best amino acid substitution model was identified for the further phylogenetic analysis using ProtTest (v3.4) (Abascal et al., 2005). The LG+I+G+F model was recommended as the best model based on both the Akaike information criterion (AIC) and Bayesian information criterion (BIC). Phylogenetic analyses were conducted for all aligned protein dataset using site-homogeneous (ML and Bayes) and -heterogeneous models (PhyloBayes). RAxML 7.0.4 (Guindon et al., 2010) was used to constructed ML trees, and 1000 bootstraps (BS) were used to estimate node reliability. The Bayesian phylogenetic tree was reconstructed using MrBayes 3.2.2 (Ronquist et al., 2012). In MrBayes analysis, the LG model was not available; thus, the JTT+I+G+F model was used as the secondary best choice (Yue et al., 2014). Two independent runs with four chains (three heated and one cold) were conducted concurrently for 1,000,000 generations, sampling every 100 generations. When the estimated sample size (ESS) value exceeded 100 and the potential scale reduction factor (PSRF) was close to 1.0, stationarity was considered to be reached, as recommended by the MrBayes software (Ronquist et al., 2012). Next, the first 25% samples were discarded as burn-in, and the branch lengths and posterior probabilities (PP) were calculated in a consensus tree. Bayesian analyses with a site-heterogeneous model were implemented using PhyloBayes 4.1b (Lartillot et al., 2013). After the removal of constant sites from the alignment, two independent chains, starting from a random tree, were run under the CAT-GTR model.

Next, we also used ASTRAL II, a statistically consistent algorithm to estimate the species tree topology under the multi-species coalescent model (Mirarab et al., 2016). Support values were obtained by calculating the local posterior probability, a advantage of ASTRAL II that presented high precision compared to multi-locus bootstrapping on a wide set of simulated and biological datasets (Sayyari and Mirarab, 2016). We conducted the analysis without the “species map” set, which means that it is unnecessary to assign multiple individuals to the same species to one taxon. Species tree obtained from RAxML gene trees of each of the 3070 orthologous loci, summarized with ASTRAL II. As species tree analyses do not require outgroups (Heled and Drummond, 2009), all 19 species thus were included as ingroups in the ASTRAL II analyses.

Based on the ML tree reconstructed above, Bayesian estimation was performed for the divergence time of the sampled 19 species to construct a time frame of cephalochordate evolution. The nine divergence nodes were signed as calibration points according to previously summarized fossil records (Benton and Donoghue, 2007; Yue et al., 2014). A calibration point (520.00–? Mya) from our early research was used at Olfactores (n) (Chen et al., 1999). In molecular dating analysis, a safe upper bound for the root age is essential in Markov chain Monte Carlo Tree (MCMCTree). Since no reliable fossil record exists that could be used here, the secondary calibration approach was used to obtain the root age (q, Deuterostomia). Currently, the divergence time at the Deuterostomia node was estimated based on the tree with the root node of Bilateria (Deuterostomia+Mollusca) via molecular clock. However, various values have been reported for the maximum constraints of the crown bilaterian divergence time, including 581.5 Mya (Benton and Donoghue, 2007), and 640–730 Mya (Peterson et al., 2008). Yue et al. (2014) explored the divergence time of crown bilaterian via MCMCTree analysis, including well-described fossil calibrations, under fixed root lower bounds (531.50) and a variety of upper bounds (600, 700, 800, and 900 Mya). Then, the authors proposed that 600 and 700 Mya were reliable and calculated their respective divergence time of Deuterostomia based on both calibration constraints (resulting in 532.67–598.27 and 585.50–698.50 Mya) (Yue et al., 2014). Consequently, we employed both results as different secondary calibration dates for the Deuterostomia (q) node in the current investigations. All calibration constraints of divergence nodes are presented in Table 2.

Molecular dating analysis was conducted using the MCMCTree subprogram of the PAML package (v4.8) with the concatenated gene matrix used in the phylogenetic analysis (Dos and Yang, 2011). As fixed topology from the guide tree, we first obtained the ML estimates of the branch lengths, the gradient (G) vector, and the Hessian (H) matrix using the codeml from the PAML package (v4.8) (Dos and Yang, 2011), then, these estimated values were fed into control files of MCMCTree together with the LG substitution model. We set a time unit of 1 million years, and used the same gamma prior G(1, 100) to specify the substitution rate and rate drift parameters. An autocorrelated and lognormal relaxed clock model (Rannala and Yang, 2007) was used to estimate the posterior distribution of the divergence time, given these priors. The burn-in for the MCMC chain was run 10,000,000 generations, with sampling parameters every 1,000 generations. Two independent runs were conducted to ensure convergence of MCMC chains. The program summarized the mean and the 95% confidence intervals (CIs).

The penalized-likelihood (PL) approach was used in the program r8s V1.8 (Sanderson, 2003), in combination with the truncated-Newton algorithm (Sanderson, 2002). The cross-validation method in r8s was used to determine the optimal level of rate-smoothing of the PL analyses with smoothing parameters varying from 1 to 1,000 according to a previously described procedure (Mulcahy et al., 2012). The optimal smoothing parameter (S = 100) was first determined using the same fixed topology (ML) used as input for MCMCTree, since there is evidence that branch lengths are more accurately estimated by ML than by BI methods (Schwartz and Mueller, 2010). CIs for the PL age estimates were obtained by replicating PL analyses of 1,000 trees. The mean divergence time and 95% CIs were summarized using the r8s bootstrap kit.

To identify the domains contained in the protein sequences for all amphioxus, the amino acid sequence sets of five species of amphioxus were searched in Pfam-A database 31.0 using local HMMER 3.1¹. This is the best program for domain identification and it has often been used together with a profile database, such as Pfam (Rekapalli et al., 2009). The domain ID set of each amphioxus was fed into a Venn Diagram² to obtain specific and shared domain types among the five investigated species of amphioxus. It is worth pointing out that only B. belcheri and B. floridae protein sets were generated from the sequenced whole genome, while those of the other three species were obtained from de novo assembled adult transcriptomes that contained assembly defects and missing protein sequences. However, a domain type could be contained in a variety of protein sequences and the domain size was short (generally in the range of 40–50 amino acids), which effectively avoided bias caused by incomplete species protein sets and de novo assembly errors. Domain sequences were annotated to GO databases using the Blastp tool with default E-value. GO enrichment analysis for specific and shared domains among all amphioxus was performed via Fisher’s exact test, that was implemented in the Blast2GO pipeline (Conesa et al., 2005). The Benjamini-Hochberg (BH) method was used for false discovery rate (FDR) correction for the Fisher test (FDR threshold values = 0.05). The list of significant GO terms was further filtered using GO trimming (v2.0) (Jantzen et al., 2011) to discard redundant terms. Pipeline of primary software used in this study was presented in Figure 1A.

Illumina sequencing for B. japonicum generated ∼52 million raw reads. After quality control for the original data, ∼49 million clean reads were obtained, and the Q20 percentage reached 97.09%, indicating a high quality of transcriptome sequencing (Table 1). Clean data have been submitted to the NCBI SRA database: Accession number SRX448311. By Trinity initial assembly, we obtained 161,542 contigs with an N50 value of 816 bp and an average length of 375 bp. Further overlapping clustering for contig sequences, we finally obtained 92,003 unigenes with an N50 value of 1,753 bp and an average length of 785 bp. The assembly results of A. lucayanum were nearly consistent with previous reports that originally sequenced its transcriptomes (Yue et al., 2014). Furthermore, we aligned pair-end reads back to unigene sequences of B. japonicum, 86.15% clean reads could be completely matched, indicating that the assembly quality meets requirements for subsequent analysis. Based on unigene set obtained, 46,145 and 4,171 CDS sequences were obtained by homologous alignment search and ESTScan detection, respectively. 46,540, 31,941, 30,319 unigenes of B. japonicum were annotated into NR, Swiss-Prot, GO database (Table 1). For genomic sequences corrected, their CDS sequences are fully identical to their corresponding proteins, ensuring reliability of subsequent analysis. The summary statistics of gene correction are shown in Supplementary Table S1.

To further clarify the phylogenic relationships among amphioxus, we performed a phylogenetic analysis with five species of cephalochordates (B. belcheri, B. japonicum, B. lanceolatum, B. floridae, and A. lucayanum), 11 vertebrate species (human, mouse, American alligator, green lizard, chicken, western clawed frog, coelacanth, fugu, zebrafish, elephant shark, and lamprey), one urochordate (sea squirt), and two outgroups (acorn worm and sea urchin). A 19-way concatenated multiple protein alignment, based on 397 orthologous gene groups intersected by HaMSTR and BRH methods (442 for HaMSTR and 415 for BRH), was used to construct the phylogenetic tree. Furthermore, the branch length (the expected amino acid substitution rate) was calculated for each species based on ML (Figure 1B), BI (Supplementary Figure S1), and PhyloBayes (Supplementary Figure S2). These three inference methods obtain congruent tree topologies with nearly identical branch lengths (Figure 1B); nodal support values are generally high (PP = 0.95–1.0, BS = 86–100); moreover, these values are stronger in the BI tree than in the ML tree, as has been previously reported (Yuan et al., 2015; Yuan et al., 2016). The species tree estimated by ASTRAL-II is very similar to the tree estimated in the concatenated analysis, with most nodes showing 100% the local posterior probability (Figure 2). Regardless of which analyzed method (ML, BI, or PhyloBayes) was used, the branch lengths of all amphioxus except for B. lanceolatum were consistently shorter than those of vertebrates, even shorter than that of the elephant shark (Callorhinchus milii) as the slowest evolving vertebrate (Venkatesh et al., 2014). The sea squirt showed the longest branch length in the phylogenetic tree. Interestingly, the branch length of B. lanceolatum is longer than those of other amphioxus and vertebrates investigated in this study. Based on this difference in branch length, the amino acid substitution rate of B. lanceolatum was assessed to be ∼1.3 times higher than that of fugu (Takifugu rubripes) with the fastest evolutionary rate among vertebrates. We conducted Tajima’s relative rate test (Tajima, 1993) based directly on the concatenated protein sequence alignment of pairwise tests between B. lanceolatum and fugu, which demonstrated a significant result (P-value below 0.05). In addition, the phylogenetic trees consistently support a phylogeny of [(B. belcheri + B. japonicum) + (B. lanceolatum + B. floridae) + A. lucayanum), regardless of the utilized analytical method.

Based on the ML tree, we performed MCMCTree analysis for the divergence times among all 19 tested species (Figure 3). Eleven time constraints were used across the phylogenetic tree (see section “Materials and Methods”). The Deuterostomia divergence times (root age) were employed as 532.67–598.27 and 585.50–698.50 Mya, respectively, based on previous investigations (Yue et al., 2014). The results for all age estimates with 95% CI (confidence interval) and mean values are shown in Tables 2, 3. Regardless of which root age we used, the most recent common ancestor of the living amphioxus species (the divergence time between Asymmetron and Branchiostoma) was estimated to be ∼104 Mya; the early splits within each main clade occurred at ∼87 Mya; B. lanceolatum diverged from B. floridae at ∼72 Mya, and the divergence time between B. belcheri and B. japonicum was ∼61 Mya. In addition, to ensure the robustness of estimation results in the divergence time frame, R8S analysis was employed using the same set of time constraints (with a root age of 532.67–598.27 Mya) that were used for the MCMCTree estimation. All divergence results with 95% CIs intervals and mean values are shown in Table 2. Despite the obviously different divergence times in several nodes, e.g., f, k, q, and r nodes, the r8s analyses overall supported the MCMCTree results, particularly within the cephalochordate. Consequently, the divergence times calculated with the MCMCTree were reliable and were directly used for the further discussion.

Domains are independent functional evolutionary units that can be independently folded (Vogel et al., 2004); most of them have ancient original history (Itoh et al., 2007). To investigate the function of ancient and species-specific domains among amphioxus, we firstly retrieved 10,106, 10,708, 8,912, 8,313, and 9,821 domains from A. lucayanum, B. belcheri, B. floridae, B. lanceolatum, and B. japonicum, respectively, by searching for homologous domains in the Pfam-A database. Among these domains, 4,814 (group 1, G1) were shared between the genera Asymmetron and Branchiostoma (Supplementary Figure S3) and were retained and handed down from their common ancestor.

Furthermore, gene ontology (GO) enrichment analysis for G1 found that 162 enriched GO terms belonging to the “biological process” subcategory are primarily related to development, cellular process and function, and metabolic processes (Supplementary Table S2). Furthermore, we detected 492 (group 2, G2), 361 (group 3, G3), 427 (group 4, G4), 255 (group 5, G5), and 366 (group 6, G6) species-specific domains in B. belcheri, B. floridae, B. lanceolatum, B. japonicum, and A. lucayanum, respectively (Supplementary Tables S3–S7). To examine whether these species-specific domains had biological functionality, we also investigated the enrichment of GO terms for those in each of G2-6 (Supplementary Table S8). Among 21 enriched GO terms in G2 (B. belcheri), those associated with the immune, inflammatory, stimulus response, apoptosis, and phagocytosis stood out in particular. For G3 (B. floridae), 42 GO terms were enriched, and those associated with development, immunity, and apoptosis were overrepresented. Next, 33 GO terms enriched from G4 (B. lanceolatum) were primarily involved in lipid storage and regulation, apoptosis, differentiation and movement of immune cells, and response to stimuli. In G5 (B. japonicum), 18 enriched GO terms were primarily related to lipoprotein oxidation and metabolism, antioxidant motions, apoptosis, and immune response. The 37 GO terms enriched in G6 (A. lucayanum) were involved in organics and energy metabolism, muscle development, and immunity system process. In general, GO terms involving innate immunity and apoptosis were most commonly enriched by these species-specific domains, follow by GO terms related to lipid metabolism and regulation as well as tissue development, suggesting a central role of these related domains in amphioxus speciation.