Whole‐genome resequencing reveals the pleistocene temporal dynamics of Branchiostoma belcheri and Branchiostoma floridae populations

Abstract Global climatic fluctuations governed the ancestral demographic histories of species and contributed to place the current population status into a more extensive ecological and evolutionary context. Genetic variations will leave unambiguous signatures in the patterns of intraspecific genetic variation in extant species since the genome of each individual is an imperfect mosaic of the ancestral genomes. Here, we report the genome sequences of 20 Branchiostoma individuals by whole‐genome resequencing strategy. We detected over 140 million genomic variations for each Branchiostoma individual. In particular, we applied the pairwise sequentially Markovian coalescent (PSMC) method to estimate the trajectories of changes in the effective population size (N e) of Branchiostoma population during the Pleistocene. We evaluated the threshold of sequencing depth for proper inference of demographic histories using PSMC was ≥25×. The PSMC results highlight the role of historical global climatic fluctuations in the long‐term population dynamics of Branchiostoma. The inferred ancestral N e of the Branchiostoma belcheri populations from Zhanjiang and Xiamen (China) seawaters was different in amplitude before the first (mutation rate = 3 × 10−9) or third glaciation (mutation rate = 9 × 10−9) of the Pleistocene, indicating that the two populations most probably started to evolve in isolation in their respective seas after the first or third glaciation of the Pleistocene. A pronounced population bottleneck coinciding with the last glacial maximum was observed in all Branchiostoma individuals, followed by a population expansion occurred during the late Pleistocene. Species that have experienced long‐term declines may be especially vulnerable to recent anthropogenic activities. Recently, the industrial pollution and the exploitation of sea sand have destroyed the harmonious living environment of amphioxus species. In the future, we need to protect the habitat of Branchiostoma and make full use of these detected genetic variations to facilitate the functional study of Branchiostoma for adaptation to local environments.

seawaters was different in amplitude before the first (mutation rate = 3 × 10 −9 ) or third glaciation (mutation rate = 9 × 10 −9 ) of the Pleistocene, indicating that the two populations most probably started to evolve in isolation in their respective seas after the first or third glaciation of the Pleistocene. A pronounced population bottleneck coinciding with the last glacial maximum was observed in all Branchiostoma individuals, followed by a population expansion occurred during the late Pleistocene. Species that have experienced long-term declines may be especially vulnerable to recent anthropogenic activities. Recently, the industrial pollution and the exploitation of sea sand have destroyed the harmonious living environment of amphioxus species. In the future, we need to protect the habitat of Branchiostoma and make full use of these detected genetic variations to facilitate the functional study of Branchiostoma for adaptation to local environments.

| INTRODUC TI ON
The Quaternary Period followed the Neogene Period and spans from approximately 2.58 million years ago (Mya) to the present. It can be divided into two epochs: the Pleistocene (2.58 Mya to 11.7 kya) and the Holocene. During the Quaternary Period, global climatic fluctuations influenced the demography and distribution of biodiversity on Earth (Farrand, 1984;Nadachowska-Brzyska, Li, Smeds, Zhang, & Ellegren, 2015). In glacial periods, many species suffered from extreme contraction or disappearance of their survival habitats (Hewitt, 2000;Hewitt, 2004;Holm & Svenning, 2014). As a consequence, these species went extinct or had to migrate to new areas where they survived in glacial refugia and adapted to the new conditions (Arenas, Ray, Currat, & Excoffier, 2012). In interglacial periods, large areas of suitable habitats became available again giving rise to the re-expansion of many species. The frequent species migration in glacial and interglacial periods likely led to several different genetic consequences (Arenas et al., 2012;Hewitt, 2004). For example, many previously large populations may have undergone a decrease in genetic diversity due to genetic bottlenecks in glacial refugia or the "founder effect" during postglacial expansions (Provine, 2004;Templeton, 1980). Additionally, the geographical separation of populations in glacial refugia likely accelerated the population differentiation and, in some cases, led to allopatric speciation (Hewitt, 2004).
Some newly formed species came into contact with other species during interglacial expansions and sometimes formed stable hybrid regions (Hewitt, 2010).
Understanding ancestral demographic histories and climatic changes allow us to place the current population status into a more extensive ecological and evolutionary context. Many studies have demonstrated that genetic variation in abundance leaves unambiguous signatures in the patterns of intraspecific genetic variation in extant species since the genome of each individual is an imperfect mosaic of the ancestral genomes (Hewitt, 2004;Li & Durbin, 2011;Sheehan, Harris, & Song, 2013). As a consequence, we can infer the historical effective population sizes (N e ) using a whole-genome resequencing strategy. In the first decade of the 21st century, many approaches were developed to analyze the demographic history of species using genetic data (Beaumont & Zhang, 2002;Csillery, Blum, Gaggiotti, & Francois, 2010;Emerson, Paradis, & Thebaud, 2001;Jody & Machado, 2003). However, these methods are restricted to special molecular markers, that is, a limited number of genomic loci instead of whole-genome data (Drummond, Rambaut, Shapiro, & Pybus, 2005;Heled & Drummond, 2008;Jody & Rasmus, 2004;Nikolic & Chevalet, 2014;Pybus, Rambaut, & Harvey, 2000). To accurately estimate the demographic history and timing of speciation, it is important to sample large numbers of genomic loci throughout the whole genome (Arbogast, Edwards, Wakeley, Beerli, & Slowinski, 2002;Ballard & Whitlock, 2004;Edwards & Beerli, 2010).
In 2011, the pairwise sequentially Markovian coalescent (PSMC) model was developed to estimate the changes in effective population size through time using whole-genome data (Li & Durbin, 2011).
The PSMC model uses a hidden Markov model to detect the historical recombination events across a single diploid genome. The hidden state of a given position corresponds to the coalescence time of the two lineages at that position, while the observed state corresponds to the observed genotype (homozygous/heterozygous) at that position. Since species evolve with the genome sequence, the coalescence time may change as a result of recombination. Moreover, the spatial distribution of homozygous and heterozygous loci will provide information on the distribution of coalescence times, which depends on the past population sizes. Although the PSMC model can produce a reasonably accurate estimation of the effective population sizes, estimating population sizes in the most recent past (<20 kya) or more anciently than 3 Mya ago is hampered, because only few recombination events are left in the present sequence.
The lack of sufficient recombination events in the most recent past reduces the power of the PSMC model to satisfactorily infer the species demographic history below 20 kya (Li & Durbin, 2011).
The PSMC model has been widely used in many vertebrate genome sequencing projects. For example, in the Tibetan macaque, giant pandas, and crested ibis, PSMC was used to investigate the demographic histories of these endangered species (Fan et al., 2014;Nadachowska-Brzyska et al., 2015;Zhao et al., 2013); in brown rat (Rattus norvegicus), PSMC was used to detect the ancestral population declines or expansions to assist in the interpretation of population genomic statistics (Deinum et al., 2015); in pied flycatcher, collared flycatcher, and dogs, PSMC was used to complement other demographic history inferences (Freedman et al., 2014;Nadachowska-Brzyska et al., 2013). However, in spite of its wide usage in vertebrates, PSMC approach has not been applied in cephalochordates to infer demographic history.
Cephalochordate, commonly called amphioxus or lancelet, has long been regarded as the phylogenetic model organism for understanding the evolution of vertebrate development (Holland et al., 2008;Huang et al., 2014;Marletaz et al., 2018;Putnam et al., 2008;Yuan, Zhang, Guo, Wang, & Shen, 2015). Lancelets, small vaguely fishshaped creatures (Figure 1), are mainly distributed in shallow subtidal sand flats in temperate, subtropical, and tropical seas around the world (Carvalho, Lahaye, & Schubert, 2017). They live in sandy bottoms whose granulometry depends on the species and the K E Y W O R D S Branchiostoma belcheri, Branchiostoma floridae, climatic fluctuations, demographic history, effective population size, genomic variations site, and they are usually found half-buried in sand (Escriva, 2018).
Unlike the much-reduced metamorphosis of the urochordates, amphioxus species seem to have been maintained in a morphological stasis for several hundred million years. This can be evidenced by the fossil record, which has brought to light a number of Cambrian soft-body fossils that closely resemble the extant amphioxus (Garcia-Fernandez & Benito-Gutierrez, 2009). Therefore, studies on amphioxus species are useful for investigating the conserved patterning mechanisms for chordates (Schubert, Escriva, Xavier-Neto, & Laudet, 2006). Currently, over 30 amphioxus species have been recorded worldwide and can be subdivided into three genera: Branchiostoma, Epigonichthys, and Asymmetron. The majority of amphioxus species are in the genus Branchiostoma. In some areas, the amphioxus populations were previously quite large. Before the 1950s, they were commercially harvested from fisheries in Xiamen, China. However, during the industrial era, pollution and exploitation of sea sand destroyed the harmonious living environment for the largest populations of amphioxus, and thus, amphioxus has been listed in the registry of "State Second-Class Protected Animals" in China and "Endangered Animals of Japanese Marine and Fresh Water Organisms" in Japan (Hantao, Yuanyuan, Chen, Fan, & Yuwu, 2005;Kubokawa, Azuma, & Tomiyama, 1998).

| Ethics statement
All sample collecting and processing were performed according to the local laws governing the welfare of invertebrate animals and were approved by the Southeast University (Nanjing, China) ethical committee.

| Sequence data preprocessing and mapping
In order to filter the low-quality paired reads, the sequencing raw reads were first filtered with the following criteria: (a) Reads with unidentified nucleotides (N) exceeding 10% were removed, and (b) reads with more than 50% bases having a phred score ≤ 5 were removed. After quality control, all of the clean paired-end reads of B. belcheri and B. floridae were mapped to B. belcheri v.18h27 and B. floridae v.2.0 reference genomes, respectively, using BWA-MEM (v0.7.15) with the following parameters: -M, -k 19 (Li, 2013).
Then, we used sambamba to filter secondary and supplementary F I G U R E 1 The amphioxus of Branchiostoma belcheri from Xiamen, China. The picture was taken in the laboratory alignments from the above-generated SAM files with the following parameters: -F "not (secondary_alignment or supplementary)" -p -l 9 (Tarasov, Vilella, Cuppen, Nijman, & Prins, 2015). The filtered SAM format files were sorted and converted into BAM format files using SAMtools v1.3.1 (Li et al., 2009), and then imported to the GATK MarkDuplicates module (v.4.0.2.1) to mark duplicate reads (McKenna et al., 2010). Finally, we used the Qualimap bamqc tool to estimate the mapped reads and their genomic properties (Okonechnikov, Conesa, & Garcia-Alcalde, 2016).

| Variant detection, filtration, and annotation
In order to improve the quality of variants reported, we performed variant calling following the GATK's best practice pipeline, includ- With the generated dbSNP dataset, we applied GATK v.4.0.2.1 "BaseRecalibrator" and "ApplyBQSR" modules to generate quality scores recalibrated bam files for every individual. Finally, we applied GATK v.4.0.2.1 "HaplotypeCaller" module to separately detect variation from the bam file generated from BQSR process for each individual.
After obtaining the genotype calls from the GATK's best practice pipeline, we applied several data quality filters to control the data quality. We first applied the GATK v.4.0.2.1 "SelectVariants" module to split the variants into SNVs and indels from the generated VCF (SNVs) and indels with much lower or higher read depth than expected were removed because lower depth and large copy-number variation might lead to overcalling and miscalling of variants. We set the minimum and maximum read depths (DP value) for variant calling as 0.1-fold and twofold of the average depth of sequencing, respectively . Additionally, we also filtered the sites that were an "N" in the reference genome.
We downloaded the genome annotation file in gff3 format and protein and transcript files in fasta format of B. belcheri v.18h27 and B. floridae v.2.0 genomes from the NCBI Genome Database. After removing duplicate records, pseudogenes, noncoding RNA genes, and genes without a transcript ID, the final annotation files for B. belcheri v.18h27 and B. floridae v.2.0 contained 25,063 and 28,621 genes, respectively. The final genome annotation files were then used in ANNOVAR v.2016-02-01 software (Wang, Li, & Hakonarson, 2010) and SNPeff v.4.3 software (Cingolani et al., 2012) with default settings to classify the called variants into intergenic regions, untranslated regions (UTRs), upstream and downstream regions (within 1 kb of a transcription start or end site), intronic regions, splicing sites (within 2 bp of a splicing junction), and exonic regions. The SNVs located in exonic regions were further classified into nonsynonymous, synonymous, stop-codon gain, and stop-codon loss variants. The indels located in exonic regions were further classified into nonframeshift, frameshift, stop-codon gain, and stop-codon loss variants.

| Inference of demographic history using PSMC
We applied the Pairwise Sequentially Markovian Coalescent (PSMC) method developed by Li and Durbin (2011) to infer the trajectory of the demographic history across time for Chinese and American lancelets (Li & Durbin, 2011). PSMC analysis requires a whole-genome diploid consensus sequence in fastq format. Considering the significant influence of different sequencing depths, generation times, and per-generation mutation rates of species, we used different settings to optimize the time estimation of the PSMC results. To evaluate the influence of sequencing depth for PSMC results, we used the SAMtools view command to subsample the high-depth samples with the following parameters: samtools view -bs FLOAT (integer part as seed). We then presented PSMC plots of 15×, 20×, 25×, 30×, 35×, 40× reads for these high-depth samples (BB_ZJ1, BB_XM1, BF_CN1, and BF_TP1), respectively. For each of these samples, we obtained the consensus sequences using the SAMtools mpileup command and the BCFtools call command with the following parameters: "samtools mpileup -C50 -q25 -uf *.fa (reference genome) *.bam | bcftools call -c -| vcfutils.pl vcf2fq -d minDP -D maxDP> * .psmc.fq". To minimize the impact of collapsed regions in the whole-genome, we excluded all sites for which read depth was less than one-third of the average read depth and more than twice the average read depth across the genome using the parameters -d and -D. After that, over 75% of the genome are reserved for the downstream PSMC analysis of all B. belcheri and B. floridae individuals. We then used the program fq2psmcfa integrated into the PSMC software to transform the consensus sequence (*.psmc.fq) into the fasta-like format files for downstream analysis. The program "psmc" was used to infer the demographic histories from the individual genomes, and the settings were chosen manually according to suggestions given by Li and Durbin (Li & Durbin, 2011). The maximum number of iterations was set to 25 by the -N option, the initial theta/rho ratio value to 5 by the -r option, and the upper limit of the most recent common ancestor (TMRCA) was first set to 13 by the -t option.
The effective population size was inferred across 61 free atomic time intervals by the -p option: "1*2 + 58*1 + 2*2", which means that the first population size parameter spans the first two atomic time intervals, each of the next 58 parameters spans one interval, and the last two parameters span two intervals, respectively. We used the program psmc_plot.pl to visualize the trajectory of demographic history across time for every individual.
The time estimation of demographic history can be significantly influenced by the mutation rate and generation time. We estimated the generation time to be 1.5 years for B. belcheri and 1 year for B. floridae based on previous studies (Igawa et al., 2017;Theodosiou et al., 2011;Yu & Holland, 2009). In this study, all the results assume that the generation time in modern Branchiostoma is representative of the ancestral population. The per-generation mutation rate is another important factor for the time estimation of demographic history. Previous studies have documented that the per-generation mutation rate scales inversely with genome size in DNA-based microbes, but scales positively in cellular organisms (Drake, 1991;Lynch, 2010). The best fitting model for per-generation mutation rate was 2.585324 × 10 −10 × g 0.584 for genome size g in Mb (Crellen et al., 2016;Lynch, 2010), giving the estimates of 9 × 10 −9 for B. belcheri and 1 × 10 −8 for B. floridae. Previous studies of amphioxus genome projects have estimated the mutation rates of B. belcheri and B. floridae to be on the order of 10 −9 to 10 −8 substitutions/site/ generation (Huang et al., 2014;Putnam et al., 2008). To account for the potential bias of mutation rates, we also performed demographic analysis with other two mutation rates (per site per generation) for B. belcheri and B. floridae: 3 × 10 −9 and 6 × 10 −9 . We compared the PSMC results of different mutation rates and discussed their corresponding outcomes in the results section. Finally, we performed 100 bootstrap replicates to assess uncertainty in effective population size estimates for each individual. Bootstrapping was conducted by randomly sampling with replacement 5-Mb sequence segments obtained from the consensus genome sequence using the splitfa script provided with the PSMC software. genes, and about 8% (41,496,865 sites) of the genome were marked as "N" (Putnam et al., 2008

| Detection and annotation of genetic variations in B. belcheri and B. floridae
The alignment files of each individual were further processed with our customized genotyping pipeline to produce the genotype calls.
Then, we applied stringent variant filtrations to further control the genotype calls (Liu et al., 2018;Moran et al., 2019). As shown in Table S1, an average of 12.82% false-positive variants was filtered in ten B. belcheri individuals and an average of 11.44% was filtered in ten B. floridae individuals. After filtering false-positive variants and large indels (>50 bp), we generated the final variants for each Branchiostoma individual. As shown in Table 2

| Estimation of generation time for Branchiostoma
As previously discussed, the growth rate of Branchiostoma varies in different species, thus contributing to different generation times of species (Chen, Shin, & Cheung, 2008;Stokes & Holland, 1995). The body-length ranges of one-, two-, and three-year-old B. belcheri were estimated to be 5-28, 28-38, and 38-45 mm, respectively (Chen et al., 2008). In contrast, the Florida lancelet B. floridae has a rapid growth rate and its body length was observed to be over 30 mm within only one year (Stokes & Holland, 1995, 1996. As described in previous studies, the minimum body length at sexual maturity is about 25-28 mm for B. belcheri and the mature B. belcheri individuals mainly belong to 1-and 2-year-old groups (Chen et al., 2008). In contrast, the Florida lancelet B. floridae has a rapid growth rate and its body length was observed to be over 30 mm within only 1 year.
As described in previous studies, the minimum body length at sexual maturity is about 25-28 mm for B. belcheri and the mature B. belcheri individuals mainly belong to 1-and 2-year-old groups. Holland (1995, 1996) reported that the minimum mature individual of B. floridae had a body length of about 23 mm body and the mature individuals mainly belonged to the one-year-old group (Stokes & Holland, 1995, 1996. From the above, we estimated a generation time of 1.5 years for B. belcheri and 1 year for B. floridae to infer the ancestral N e dynamics in the following PSMC analyses.

| Influence of sequencing depth on PSMC analysis
As our data included individuals for which the mean sequencing varied from 11.48× to 122.87× (Table 1)
As shown in Figure 5a, when applying the mutation rate of Then, the Zhanjiang population maintained a stable N e during the interglacial period between the third and last glaciations. In contrast, F I G U R E 5 Temporal dynamics of effective population size for seven high-depth B. belcheri individuals. (a) The generation time (g) and the mutation rate (μ) are assumed to be 1.5 years and 3 × 10 −9 , respectively. (b) The generation time (g) and the mutation rate (μ) are assumed to be 1.5 years and 9 × 10 −9 , respectively. The four light green areas

| D ISCUSS I ON
The ratio of the number of transitions to the number of transversions (Ti/Tv) for a pair of sequences becomes 0.5 when there is no bias toward either transitional or transversional substitution. The Ti/ Tv ratio of Branchiostoma was ~1.3, which was much less than the criterion of many mammal genomes, but a little higher than that of some fishes (DePristo et al., 2011;Lachance et al., 2012). Although transversion mutations (A ↔ C, A ↔ T, C ↔ G, and G ↔ T) should be twice as much as transition mutations (A ↔ G and C ↔ T), transition mutations are more likely than transversions because substituting a single ring structure for another single ring structure is more likely than substituting a double ring for a single ring. As shown in Figure 3, the Ti/Tv ratio ranged from the lowest of 0.562 in mosquito to the highest of 2.389 in sheep. The extremely low Ti/Tv ratio of Branchiostoma population most probably resulted from the very low levels of CpG methylation in its genome, which has been proved by a recent publication (Marlétaz et al., 2018). In human and some other mammal genomes, CpG methylation is a very common characterized epigenetic modifications that cells use to control gene expression, and Cytosine (C) on the genome is prone to C -> T transitions under methylation modification (Cooper, Mort, Stenson, Ball, & Chuzhanova, 2010;Li & Zhang, 2014). Additionally, transitions in genomes are less likely to result in amino acid substitutions (due to wobble base pair) and are therefore more likely to persist as synonymous SNVs in populations.
The PSMC approach is an excellent method to infer detailed ancestral population size dynamics over time using information from whole-genome sequences and has recently been used in many genome projects. However, as described by Mazet et al. (Chikhi et al., 2018;Mazet, Rodriguez, Grusea, Boitard, & Chikhi, 2016), when considering structured populations, the population structure can mimic the signature of any demographic scenario, because PSMC assumes a panmictic model and the gene flow patterns could result in the changes of estimated N e (Chikhi et al., 2018;Gautier et al., 2016). Therefore, it is not possible to distinguish between a scenario involving habitat fragmentation and a scenario of actual population size changes from the PSMC analyses. It also should be noted that the inferred N e over time largely depends on the sequencing depth and the parameters used in the PSMC analysis (Li & Durbin, 2011;Nadachowska-Brzyska, Burri, Smeds, & Ellegren, 2016). As described by Nadachowska-Brzyska et al. (2016), an appropriate sequencing depth for proper inference of demographic histories using PSMC is ≥18×. However, we found the threshold (18×) of sequencing depth was not enough to get a proper PSMC curve for Branchiostoma. The best threshold of sequencing depth for Branchiostoma was ≥25× ( Figure 4), which most probably due to its high heterozygosity to fail call heterozygous sites with low-depth data ( Figures S3 and S4).
The recombination rate, per-generation mutation rate, and generation time estimates are three important parameters for the PSMC results, and the inference of ancestral N e can be biased if they are under-or overestimated (Nadachowska-Brzyska et al., 2015. The recombination rate will not influence the accuracy of PSMC results, when it was set appropriately. As shown in Figure Figure S6); however, when given a fixed generation time, a doubling of the per-generation mutation rate will halve the estimate of N e and also move the PSMC curve backward in time ( Figure S7). In this study, we assumed that genetic variation in modern Branchiostoma species is representative of the ancestral population. Based on recent studies about amphioxus growth and gonad development, we estimated the generation time to be 1.5 years for B. belcheri and 1 year for B. floridae to infer the ancestral N e dynamics (Stokes & Holland, 1995, 1996. Additionally, previous studies have found that the growth and gonadal development of Branchiostoma are closely related to seawater temperature, and the higher the seawater temperature, the faster the growth and the earlier the gonadal development (Igawa et al., 2017;Theodosiou et al., 2011;Yu & Holland, 2009). Therefore, the effective population size of Branchiostoma should decline during the most glacial periods, especially during the last maximum glaciation. As shown in Figures 5 and 6, only the effective population size generated by a per-generation mutation rate of 3 × 10 −9 suffered from a large reduction during glacial periods and remained stable during the interglacial periods of the Pleistocene. Therefore, the per-generation mutation rate of 3 × 10 −9 is more likely to represent the changes in the historical effective population sizes over time than the other two higher mutation rates.
The PSMC approach has recently been applied to infer demographic history in many genome projects but has never been applied   (Miller et al., 2012). In this study, we also detected population reductions in the Branchiostoma populations at the beginning of the LGP and found a population bottleneck during the LGM (Figures 5 and 6). At the beginning of the LGP, the PSMC curves of the seven B. belcheri individuals from China all showed a severe population reduction, while the three B. floridae individuals only experienced a minor population reduction during this period (mutation rate = 3 × 10 −9 ). This may be due to the fact that the B. floridae population from America was subject to a very long-term population reduction before the LGP, and thus, the population size was very small during the LGP. Moreover, the time points of population bottlenecks between B. belcheri and B. floridae were slightly different, which may be due to the different levels of glaciation in Eastern China and North America. However, it should be noted that the resolution of the PSMC analysis was poor in recent times (including the LGP), especially since 20 kya (Li & Durbin, 2011;Nadachowska-Brzyska et al., 2015). We performed 100 bootstrap replicates to assess the robustness of our PSMC results during the Pleistocene. As shown in Figures 5 and 6, the bootstrapping PSMC curves confirmed the accuracy of our methods in most of the Pleistocene. However, in the more recent past (<20 kya), there is considerable variation in effective population size estimates among bootstrap replicates for the most recent time intervals, confirming the limitations of PSMC analysis to infer demography in the more recent past (Li & Durbin, 2011).
Amphioxus has always been considered as a promising model animal to study the evolutionary and developmental mechanisms of vertebrates because of its unique phylogenetic position and simple body plan. Amphioxus species are distributed throughout the world along tropical and temperate coasts in habitats with a suitable temperature and salinity and a sandy seabed. Before the 1960s, along parts of the coast of China, amphioxus species were so numerous that they constituted the basis of a fishing industry (Chin, 1941;Light, 1923). Although edible, amphioxus species have never been sufficiently abundant to constitute a significant source of food to humans or to be an important part of the food chain in nature. With the severe industrial pollution and the construction of seawalls, the habitats of amphioxus species have been destroyed. At present, amphioxus species have been listed in the registry of "State Second-Class Protected Animals" in China and "Endangered Animals of Japanese Marine and Fresh Water Organisms" in Japan. We are currently in a warm interglacial period of the Holocene, which is suitable for the breeding of amphioxus species. Therefore, we should provide sanctuaries for the amphioxus species rather than destroying their habitats for human development.

CO N FLI C T O F I NTE R E S T
The authors report no conflict of interest.