{'fraction_non_alphanumeric': 0.03891545118925135, 'fraction_numerical': 0.030472674453791683, 'mean_word_length': 5.150567653080216, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The White-Spotted Bamboo Shark Genome Reveals Chromosome Rearrangements and Fast-Evolving Immune Genes of Cartilaginous Fish

Summary Chondrichthyan (cartilaginous fish) occupies a key phylogenetic position and is important for investigating evolutionary processes of vertebrates. However, limited whole genomes impede our in-depth knowledge of important issues such as chromosome evolution and immunity. Here, we report the chromosome-level genome of white-spotted bamboo shark. Combing it with other shark genomes, we reconstructed 16 ancestral chromosomes of bamboo shark and illustrate a dynamic chromosome rearrangement process. We found that genes on 13 fast-evolving chromosomes can be enriched in immune-related pathways. And two chromosomes contain important genes that can be used to develop single-chain antibodies, which were shown to have high affinity to human disease markers by using enzyme-linked immunosorbent assay. We also found three bone formation-related genes were lost due to chromosome rearrangements. Our study highlights the importance of chromosome rearrangements, providing resources for understanding of cartilaginous fish diversification and potential application of single-chain antibodies.


INTRODUCTION
The white-spotted bamboo shark, Chiloscyllium plagiosum, (hereinafter referred to as bamboo shark) belongs to the class of Chondrichthyes, which is one of the oldest extant jawed vertebrate groups (McKenna, 1988). Cartilaginous fishes including Elasmobranchii and Holocephali shared a common ancestor with other vertebrates about 460-520 Ma, and then evolved independently to distinct lineages (Inoue et al., 2010). The phylogenetic evolution of cartilaginous fishes has been disputed for a long time (Cao et al., 1998;Janvier, 1996;Rasmussen and Arnason, 1999;Zardoya et al., 1998), especially the evolutionary relationships with bony fishes. Furthermore, most cartilaginous fishes have various chromosome karyotypes (2n = 66-104) (Rocco et al., 2003;Schwartz and Maddock, 1986), revealing interesting chromosome evolution processes. As known, immunoglobulins or lymphocyte receptors-based adaptive immunity is restricted to vertebrates (Litman et al., 2010). And as one of the extant early branching jawed vertebrates, cartilaginous fishes developed special immunity (for example, sharks comprise heavy-chain immunoglobulins, which are different from canonical antibodies consisting of both heavy and light chains, Kö nning et al., 2017), which makes them ''immunologist's delight.'' The bamboo shark is a nocturnal reef-dwelling species and widely distributed in the Indo-West Pacific from India to Indonesia, southern China, and Japan (Kyne and Burgess, 2006). Its biological features, including docile nature, small body size (24-37 inches in length), convenient reproductivity, and longevity make it ideal for research. Its special immunoglobulins are superior in biological and medical applications and have been proposed for developing antibody drugs (Wesolowski et al., 2009;Zielonka et al., 2014). Despite its biological and application importance, previous researchers mostly focused on addressing its biology of hematology, reproduction, muscle activity, liver regeneration, and anatomy (Alexander et al., 2016;Maia and Wilga, 2013;Straube et al., 2016). Limited whole-genome sequencing of cartilaginous fishes including elephant shark (2014)  ) prevents our further understanding of genetic mechanisms for these species. Furthermore, these five genomes were all assembled to scaffold level, which also limits our investigation of chromosome rearrangements.
To better understand the evolution and special immunity of cartilaginous fishes, we sequenced and assembled a chromosome-level genome of a female bamboo shark, identifying dynamic chromosome rearrangement events and related evolutionary consequences. We carefully analyzed evolution of immune-related genes, which will provide new resources to understand their biology and applications in immunology. We also found that chromosome rearrangements delete three important bone formation-related genes, which may interpret chondrified endoskeleton of cartilaginous fishes.

Genome Assembly and Annotation
We assembled a 3.85-Gb genome assembly with 51 chromosomes supported by chromatin interaction relationships with Hi-C sequencing data ( Figure 1A and Tables S1 and S2) and karyotype analysis (Ma et al., 2008), and we annotated 19,595 protein coding genes (Table S3) and 63.53% of repeat content (Table S4) in this genome assembly. Comparison of repeat content among cartilaginous fishes and bony fishes shows that cartilaginous fish genomes contain higher proportion of repeated sequences (Table S5). The GC content and repeat and gene density distributed in 51 chromosomes are shown in Figures S1 and S2. And 95.8% of the Benchmarking Universal Single-Copy Orthologs (BUSCOs) (Simao et al., 2015) were identified to be complete in this genome (Table S6). Syntenic relationships revealed unambiguous alignments of 41 bamboo shark chromosomes to 29 chromosomes of chicken (the tetrapod species with most stable karyotypes, Ellegren, 2010, Figure S3), whereas the alignments between bamboo shark and zebrafish chromosomes ( Figure S4) are disordered. Extensive inter-chromosomal rearrangements have been described previously in zebrafish genome (Kasahara et al., 2007). Therefore, these results suggest the bamboo shark genome also kept relatively conserved chromosome karyotypes without many inter-chromosomal rearrangements.
The assembled genome size of bamboo shark is~3.85 Gb, larger than elephant shark (~974 Mb) (Venkatesh et al., 2014), whale shark (~2.93 Gb) (Read et al., 2017), and most bony fishes (340 Mb-2.97 Gb). Wholegenome duplication (WGD) is one of reasons that result in larger genome sizes (Grover and Wendel, 2010). Thus we investigate whether there is a third WGD event in bamboo shark apart from the common two rounds of WGD of vertebrates (Grover and Wendel, 2010). First, we found only one peak on the 4-  iScience Article fold synonymous third-codon transversion rates (4DTv) distribution of bamboo shark genome, which represented the recent common WGD event of all vertebrates ( Figure S5). We also checked the 4DTv distribution of elephant shark and zebrafish, finding that elephant shark also has only one peak, which is similar to bamboo shark, but zebrafish has another peak representing the third WGD of teleost fish (Glasauer and Neuhauss, 2014), indicating the reliability of our results. Second, we checked HOX (homeobox) genes, which are highly conserved in vertebrates and always clustered together (Santini et al., 2003) and thus have become reliable markers of WGD events (Kuraku, 2011). We only identified three HOX clusters in bamboo shark, compared with seven clusters in zebrafish, which experienced a third WGD. Similar HOX clusters were also found in cattle, whale shark and elephant shark genomes ( Figure S6). These two results suggest that bamboo shark genome did not experience a third WGD event. Therefore, the larger genome size of bamboo shark should result from the burst of repeated sequences.

Molecular Phylogeny
In consideration of the important evolutionary position of cartilaginous fishes, we inferred phylogenetic relationship of bamboo shark with other five cartilaginous fishes, five representative bony fishes (four ray-finned fishes and one lobefin fish), and humans. Based on coding sequences of 823 single-copy orthologous genes identified using TreeFam (Li et al., 2006), we constructed phylogenetic trees using both maximum-likelihood and Bayesian methods and generated identical results ( Figure 1B

Reconstruct Ancestral Chromosome
This chromosome-level genome makes it possible to study chromosome evolution of cartilaginous fishes. Thus we reconstructed ancestral chromosome karyotypes of cartilaginous fishes by identifying paralogous and orthologous genes between the bamboo shark and elephant shark genomes (Venkatesh et al., 2014) following a previously described method (Salse et al., 2009) (Table S7). Finally, we constructed 16 putative ancestral chromosomes and illustrated an evolutionary scenario during which eight fission and five fusion events occurred (Figure 2A, colored arrows), possibly for all cartilaginous fishes. As for the bamboo shark, nine fission and four fusion events (black and dotted arrows) occurred, resulting in six candidate daughter chromosomes (Chr8, Chr29, Chr38 and Chr39, Chr45, Chr48) (Figures 2A and S7 and Table S7). All these rearrangements ultimately gave rise to 51 chromosomes of the bamboo shark genome.

Fast-Evolving Chromosomes and Immune Genes
Plenty of chromosome rearrangements play a role in fast-evolving gene families and in fostering largescale changes in gene order (Eichler and Sankoff, 2003). To identify potential causes and consequences of dynamic chromosome rearrangements in cartilaginous fishes, we further analyzed distribution of conserved protein-coding genes of cartilaginous fishes along bamboo shark chromosomes. We identified 2,323 orthologous genes (~12.90% of total genes) shared among bamboo shark, elephant shark, whale shark, brown-banded bamboo shark, cloudy catshark, and white shark ( Figure S8). After exclusion of genes shared among these six cartilaginous fishes and representative bony fishes (medaka, Kasahara et al., 2007, Figure S9, andspotted gar Braasch et al. 2015, Figure S10), we finally identified 1,359 genes conserved only in cartilaginous fishes ( Figure 2B). Interestingly, we found those genes to be unevenly distributed along bamboo shark chromosomes with conserved genes on chromosomes 8,37,39,41,43,44,45,46,47,48,49,50, and 51-notably fewer than (average: 2.6 genes) those of other chromosomes (average: 34.9 genes, (Mann-Whitney U test, p value < 0.001) (Figures 2B and Table S8). We then evaluated the evolutionary rate by calculating K S (synonymous substitutions per synonymous site) values of orthologous genes on these 13 chromosomes (mean K S value: 2.79), which was significantly higher than that of other chromosomes (mean K S value: 1.54, Mann-Whitney U test, p value < 0.001, Figure 2C). In addition, we found heterozygous SNPs in the genome of this individual we sequenced to be notably more frequent on these 13 chromosomes (except Chr43) than other chromosomes (Mann-Whitney U test, p value < 0.001, Figure 2B). All these findings suggest that these 13 chromosomes are fast-evolving. Enrichment analysis (according to the Kyoto Encyclopedia of Genes and Genomes [KEGG]-assigned gene functions and pathways) showed that genes on these 13 fast-evolving chromosomes are significantly enriched in immune-related pathways with 171 immune-related genes (p value < 0.01, Tables S9 and S10). These include allograft rejection, antigen processing and presentation, as well as intestinal immune network for IgA production.

Analysis of MHC-Related Genes
Among these 13 fast-evolving chromosomes, we found that Chr37 and Chr44 likely underwent a special self-fusion event after a possible chromosome or large segment duplication event ( Figure 2A). We also found that major histocompatibility complex (MHC) genes (11 class I and 3 class II genes) are notably enriched on Chr37 (11 genes, Figure 3A), except for those on unanchored scaffolds. MHC genes were not found in the amphioxus genome, whereas one fragment of a possible MHC class II gene was found in sea lamprey (Gene ID: PMZ_0007681-RA; KEGG function: K06752 MHC, class II) (Smith et al., 2018). Upon further investigation of MHC gene numbers in other species, we found both MHC class I and class II genes in cartilaginous fishes and bony fishes except for the elephant shark genome, which lacked MHC class II genes according to our analysis (Tables S11, S12, and S13 and Figure S11). These results suggest that the innate immune system accompanied with adaptive system based on variable lymphocyte receptors (Pancer et al., 2004;Saha et al., 2010) played a major role in defending against infections in amphioxus and sea lamprey, whereas cartilaginous and bony fishes evolved with acquiring a complete MHC-based adaptive immune system. The differences in these immune systems may have arisen from the fast-evolving chromosomes. Moreover, we suggest that MHC class II genes were likely acquired before MHC class I genes based on our identification of an MHC class II-like fragment in the starlet sea anemone genome (NCBI Accession: XP_001628845.1, identified by aligning PMZ_0007681-RA using BLAST, Altschul et al., 1990) and sea lamprey genome (Smith et al., 2018), potentially resolving a long debate about MHC evolution (Flajnik et al., 1991;Kaufman, 1988Kaufman, , 2011Kaufman, , 2018Kaufman et al., 1984;Rock et al., 2016;Zhang et al., 2014). In addition, we found that tripartite motif-containing protein 69 (TRIM69) gene family was expanded significantly in cartilaginous fishes (average 18 copies) compared with ray-finned fishes (less than 3 copies) ( Figure S12). Also, in bamboo shark genome, 13 copies of TRIM69 were also located on

Identification and Amplification of IgNAR
In contrast to MHC genes found on Chr37, we identified that the immunoglobulin new antigen receptor (IgNAR) (Feige et al., 2014) gene loci (four complete IgNAR structure, V-CH1-CH2-CH3-CH4-CH5, and two incomplete IgNAR) was located on Chr44 ( Figure 3B). To obtain and verify the complete sequence assembly of IgNAR region, we sequenced~124-fold new single-tube long fragment reads (stLFR) (Wang et al., 2019) to re-assemble it. The good syntenic relationship ( Figure 3C) and the distribution of pairedend reads ( Figure S13) reveal high-quality assembly of IgNAR region. Because of the application potential of single-domain antibodies (sdAbs) in biotechnical and therapeutic use, we tried to check the diagnostic potential of IgNAR in bamboo shark. We first designed primers (Table S14) based on IgNAR sequences of bamboo shark and specifically amplified the variable domain of New Antigen Receptors (vNARs) from peripheral blood leukocytes and spleen tissue of five bamboo shark individuals. Amplified vNARs were inserted into phagemid vector pMECS and then into E. coli TG1 competent cells to produce the vNAR-phage display library. Sanger sequencing of randomly selected~100 clones shows low repetition and high diversity of vNARs, especially the complementarity-determining region 3 (CDR3) ( Figure S14). We then chose several human disease biomarkers: YKL40 (Kastrup, 2012; Rathcke and Vestergaard, 2009) for cardiovascular disease, CD64 for infectious disease (Hoffmann, 2009), and Ki67 for lung cancer (Li et al., 2015) as targets of four rounds screened monoclonal vNARs. Monoclonal phage enzyme-linked immunosorbent assay    iScience Article (ELISA) shows high affinity to targets of those positive clones. In detail, we obtained 13 unique YKL40-binding clones with signal at least~6-fold (highest:~20-fold) than negative control (NC) (Figures 3D), 9 unique CD64-binding clones with signal at least~3-fold than NC ( Figure 3E), and 3 unique Ki67-binding clones with signal at least~5-fold than NC ( Figure 3F). Although further study should be carried out to verify the affinity of those positive clones and their diagnostic function by using real samples, we believe that this work is of significance for using genome data to develop sdAbs.

P2X Gene Family
Chromosome rearrangements would remove genes that may first become pseudogenes because selective pressure acting on them was relaxed when new phenotypic traits arose or they may have very little effect on its adaptations. Syntenic comparison among chicken, zebrafish, medaka, and bamboo shark showed at least four possible genome rearrangement events occurred in bamboo shark that may result in the deletion of gene, p2rx5 (Purinergic Receptor P2X, Ligand-Gated Ion Channel, 5) after this gene became redundant or non-functioning due to selective pressures acting on it ( Figure 4A). The loss of this gene was also supported by checking RNA sequence data of 14 tissues including blood, eye, gill, heart, liver, muscle, spleen, stomach, dorsal fin, tail fin, pancreas, leptospira, two capsulogenous gland, and two kidney samples. And this gene has been previously reported to be involved in bone development and homeostasis (Nicolaidou et al., 2012;Sitcheran et al., 2003;Solle et al., 2001;Sun et al., 2013;Syberg et al., 2012;Takahashi et al., 1988;Thaler et al., 2014). Furthermore, analysis of the whole gene family of purinergic receptor P2X in sea lamprey, six sharks, and representative bony fishes showed that p2rx3, p2rx5, and p2rx7 genes were lost in six cartilaginous fishes, whereas at least six paralogs (p2rx1, p2rx2, p2rx3, p2rx4, p2rx5, p2rx7) with multiple copies were found in bony fishes ( Figure 4B and Tables S15 and S16  The shark's P2X genes are marked red. Although many p2rx5 genes exist in PMZ, all these genes were less than 246 aa in length (10 of 14 less than 200 aa). This tree was constructed with protein sequences of these genes by using software of MUSCLE (v3.8.31) (Edgar, 2004) and FastTree (v2.1.10) (Price et al., 2010).

DISCUSSION
Because of the therapeutic potential of single-domain antibodies, sharks have drawn scientists' interest for many years. With ideal biological features of bamboo shark, we selected this species and were able to obtain high-quality samples for further research. Combining paired-end, mate-paired, stLFR, and Hi-C sequencing strategies, we successfully assembled a chromosome-level reference genome of bamboo shark. In the present study, we mainly focused chromosome evolution and t fast-evolving chromosomes and immune gene families of cartilaginous fishes. Also, we investigated that burst of repeat that caused larger genome size of bamboo shark and inferred phylogenetic topologies between sequenced cartilaginous fishes and bony fishes, which is important for exploring evolutionary process of vertebrates.
Using this genome, we inferred ancestral chromosomes of cartilaginous fishes and found dynamical rearrangements. Based on chromosome evolutionary processes and comparative genomic analysis, we were able to identify fast-evolving chromosomes and immune-related genes. Moreover, chromosome fusions and fissions would also cause DNA damages, deletion of genes, and formation of new genes that may be functionally important and closely associated to species-specific features. Based on this, we found gene loss events associated with phenotypic diversity, for example, chondrification of the endoskeleton. Thus, our study highlights the importance of chromosome rearrangements in the diversification of cartilaginous fishes. With effective methods described in this study, more chromosome-level genomes can be obtained in the future to further elucidate the early evolution of jawed vertebrates as well as extant jawed vertebrate lineages.
Shark-specific immunocompetence always attracts researchers. Investigating immune genes could help to understand evolutionary processes of immune system of cartilaginous fishes compared with jawless species. MHC-I-like and MHC-II-like genes found in cartilaginous fishes revealed the possible time of acquirement of MHC-based adaptive immunity. Bamboo shark chromosomes (Chr37 and Chr44) enriched with immune-related gene may play vital roles in its powerful immunity, and more chromosome-level genomes of cartilaginous fishes should be accomplished to further confirm this conclusion. Shark singledomain antibodies have shown prospects in therapeutic use and our ELISA experiments also proved their potential use in human diseases. Thus, assembly of IgNAR sequences will accelerate development of antibodies for future medicine. In summary, our results provide valuable resources and will be significant for future research about vertebrate evolution and pharmaceutical development.

Limitations of the Study
In this study, we assembled bamboo shark genome, analyzed chromosome rearrangements, and performed ELISA experiments. However, this genome is the only chromosome-level cartilaginous fish genome and more high-quality genomes should be assembled to further verify our chromosome evolution analysis.
Besides, more functional experiments should be performed to further validate candidate functional genes.

Resource Availability Lead Contact
Further information and requests for materials should be directed to and will be fulfilled by the lead contact, Xin Liu (liuxin@genomics.cn).

Materials Availability
This study did not generate new unique reagents.

Data and Code Availability
The accession numbers for the genome sequencing data, RNA sequencing data, and genome assembly reported in this paper NCBI: PRJNA478295.This Whole Genome Shotgun project has been deposited to National Center for Biotechnology Information (NCBI) under the accession: QPFF00000000 referring project: PRJNA478295. Raw RNA sequencing reads have also been uploaded to the SRA database under accession: SRP154403. The assembled genome can also be obtained from CNSA (CNGB Nucleotide Sequence Archive) by assembly ID: CNA0000025.

METHODS
All methods can be found in the accompanying Transparent Methods supplemental file.

ACKNOWLEDGMENTS
We would like thank Dr. Lynn Fink for revising these manuscripts. This work was supported by National Key

DECLARATION OF INTERESTS
The authors declare no competing interests.

DNA, RNA extraction
Fourteen tissue types including blood, eye, gill, heart, liver, muscle, spleen, stomach, dorsal fin, tail fin, pancreas, leptospira, 2 capsulogenous gland and 2 kidney samples were collected from a female bamboo shark individual. Genomic DNA was extracted from one blood sample using the phenol-chloroform method and its quality and quantity were assessed by pulsed field gel electrophoresis and Qubit Fluorometer. 14 RNAs from these 14 tissue types were extracted by using TRIzol ® Reagent and were assessed by Agilent 2100 bioanalyzer system.

Library construction
Firstly, for WGS libraries with average insert sizes of 180 bp and 350 bp, a Covaris E220 ultrasonicator (Covaris, Brighton, UK) was used to shear the extracted high-quality DNA and AMPure XP beads (Agencourt, Beverly, USA) were used to select target fragments. Then, fragment end-repairing and A-tailing were performed by T4 DNA polymerase, T4 polynucleotide kinase and rTaq DNA polymerase. Next, PCR amplification of eight cycles was carried out and the single-strand circularization process was performed using T4 DNA Ligase, generating a single-stranded circular DNA library for sequencing.
Secondly, for mate-pair libraries, a Covaris E220 was used to acquire ~2 kb DNA fragments and Hydroshear (GeneMachines, CA, USA) was used to acquire ~5 kb, ~10 kb and ~20 kb DNA fragments. After further selection and purification of DNA, fragments were end-repaired and biotin-labeled. The modified fragments were circularized and re-fragmented using a Covaris E220. Biotin-labeled DNA fragments were captured on M280 streptavidin beads (Invitrogen, CA, USA), end-repaired, A-tailed and adaptor-ligated. Biotin-labeled fragments were PCR-amplified, purified on streptavidin-coated magnetic beads, size-selected by agarose gel electrophoresis and column purification, single-stranded and re-circularized. The purified PCR products were heat-denatured together with an adapter that was a reverse-complement to a particular strand of the PCR product, and the single-stranded molecule was then ligated using DNA ligase to get a single-stranded circular DNA library.
Thirdly, a blood sample was used for constructing the Hi-C library. The fresh blood cells ( required for the BGISEQ-500 platform was carried out and DNBs were prepared as previously described. Fourthly, for RNA library construction, mRNA was extracted from different tissues using TRIzol ® Reagent, fragmented, and then reverse-transcribed into cDNA by using Hiscript II Reverse Transcriptase (Vazyme Biotech, Nanjing City, P.R. China). Then, all single-stranded circular DNA libraries were constructed by using the same strategy described as above.
Last, the stLFR library was constructed by its standard protocol that was published recently (Wang et al., 2019).

Sequencing for all libraries
All sequencing data were generated using the BGISEQ-500 platform. Libraries with an average insert size of 180 bp and 350 bp were sequenced yielding paired-end reads with 100 bp in length. Mate-pair libraries with average insert sizes of 2k, 5k, 10k and 20k and Hi-C library were sequenced yielding reads with 50 bp in length.
RNA libraries were sequenced yielding paired-end reads with 100 bp in length. StLFR reads was also sequenced generating paired-end reads with 100 bp in length.
TransDecoder (Haas et al., 2013) was then used to predict the candidate complete ORFs. Further, for de novo gene prediction, we employed AUGUSTUS (v3.1 ) (Stanke et al., 2006) to scan the whole genome with a custom training set generated by using 2,000 high quality genes. Subsequently, we combined the homologybased and de novo-predicted gene sets using GLEAN (Elsik et al., 2007) and integrated the GLEAN and transcriptome results with in-house scripts to generate a representative and non-redundant gene set. The final gene set was assigned with a potential function by aligning proteins to databases including KEGG, Swissprot, TrEMBL and InterPro.

Construction of phylogenetic tree
To infer phylogenetic tree of bamboo shark, we selected other five published shark genome including elephant shark (2014), whale shark (2017), brownbanded bamboo shark (2018), cloudy catshark (2018) and white shark (2019), five representative bony fishes including four ray-finned fishes and one lobefin fish as well as human to identify single copy orthologous genes. we used TreeFam (Li et al., 2006), defining 33,306 gene families including 823 single copy orthologous genes. Based on coding sequences of these 823 genes, we constructed trees using both Maximum-likelihood and Bayesian methods with GTR+gamma model, which generated identical trees. Next we estimated divergence time of important branches using PAML mcmctree programme (Yang, 2007;Yang and Rannala, 2006) with the approximate likelihood calculation method.

Evolution of chromosomes
For ancestral chromosomes construction, we identified paralogous genes and orthologous genes by using criteria defined by Salse et al (Salse et al., 2009) with both Cumulative Identity Percentage (CIP) and Cumulative Alignment Length Percentage (CALP) value of 0.5 and selected genes pairs defined as A match B best and B match A best. Then MCSCAN (Tang et al., 2008) was used to generate synteny blocks with default parameters.
We first noted 54 shared duplications (with 414 paralogous gene pairs) on all the chromosomes. After further integration of these duplications and gene pairs, we found 16 duplicated chromosome pairs and 5 single chromosomes. For identification of conserved genes among bamboo shark, elephant shark, whale shark, brownbanded bamboo shark, cloudy catshark and white shark, spot gar, and medaka, the same criteria with both CIP and CALP value of 0.3 were used. The best hit of multiple matches was selected. Ks values of single copy genes were calculated by KaKs Calculator with default parameter. Heterozygosity of each chromosome was calculated by calling heterozygous SNPs generated by BWA (Li and Durbin, 2009) and SAMtools package (Li, 2011).

MHC genes and P2X gene family analysis
We downloaded coding sequences (CDS) and proteins of Callorhinchus milii (NCBI Accession: GCF_000165045.1), Rhincodon typus (NCBI Accession: GCF_001642345.1), Fugu rubripes (NCBI Accession: Danio rerio, Latimeria chalumnae, Oryzias latipes and Gasterosteus aculeatus from Ensembl databsse (release 84), and sea lamprey from (https://genomes.stowers.org/organism/Petromyzon/marinus) and Branchiostoma floridae from JGI Genome Portal (http://genome.jgi-psf.org/Brafl1/Brafl1.home.html). Next we carried out initial filtering by discarding sequences less than 30 amino acids and kept the longest transcript if one gene contains multiple transcripts. Subsequently, gene function annotation of these gene sets was performed by using the same method to bamboo shark by using BLAST (Altschul et al., 1990). Then we summarized the MHC genes and P2X gene families by integrating the above information. The P2X-like genes were searched against NCBI database.

Identification of IgNAR
We downloaded IgNAR (immunoglobin new antigen receptor from cartilaginous) nucleotide sequences from NCBI database and then aligned these sequences to the bamboo shark genome by using BLAST (v2.2.26) (Altschul et al., 1990) with parameters "-F F -m 8 -e 1e-05" to get high scoring pairs (HSPs). Then we clustered the HSPs by integrating the overlaps between HSPs to verify the IgNAR loci. To further verify the IgNAR loci, we re-assembled this region by using single tube long fragment reads (stLFR). Firstly, we mapped the stLFR reads to our initial IgNAR region using BWA (Li and Durbin, 2009) and selected these reads. Then we extracted all reads that were marked with the same barcodes as the mapped reads. Finally, we used Supernova (v2.1.1) (Mohr et al., 2017) to assemble these reads and got the re-assembled IgNAR region with default parameters.

Bamboo shark vNAR phage display library construction
Five bamboo shark individuals (3 male and 2 female, with body length ranging from ~31.50-39.37 inch) were bled and dissected out for spleen, which were submerged in RNAlater for subsequent storage at -80℃. Total RNA was isolated using the TRIzol reagent (Thermo Fisher Scientific) according to the manufacturer's instruction. We reverse-transcribed RNA into cDNA in a total of 20µl volume with the SuperScript III First-Strand Synthesis System (Thermo Fisher Scientific) according to the manufacturer's instruction. Based on the IgNAR sequences of bamboo shark, we designed and synthesized four pairs of primers for four types of the VNAR domain to amplify the VNAR sequence from the cDNA product by PCR. All reactions were pooled on 1.5% TAE agaose gel and purified by QIAEX II Gel Extraction kit 500 (Qiagen). Purified PCR products and pMECS vector were digested by restriction enzyme PstⅠand NotⅠ(NEB). The digested PCR products and vector were purified by QIAEX II Gel Extraction kit 500 (Qiagen). Then the PCR products were inserted into the pMECS vector at a molar ratio of 3:1 using T4 ligase (NEB). Ligation products were purified and transformed into electrocompetent TG1 cells (Lucigen Corporation). The TG1 cells were stored at -80℃ by addition of glycerol to a final concentration of 15% (v/v).

Bamboo shark vNAR phage display library screening
Nonadjacent wells of a Maxisorp 96-well were coated with 100µl of 5µg/ml human recombinant antigen Ki67, CD64 and YKL40 (Sino Biological). The plate was sealed by microplate sealing tape and was incubated overnight at 4℃. The following morning, the antigen-coated wells were washed three times with 200µl of PBST, then 200µl blocking solution of 0.5% BSA (Genview) was added in the plate which was incubated for 2h on a vibrating platform at RT. After the blocking solution removed from the antigen-coated wells, the plate was washed five times with 200µl of PBST. For each antigen, 100µl rescued phage library were dispensed into both a selection and a control well. The plate was incubated for 2h on a vibrating platform, then washed five times with 200µl PBST. 100µl triethylamine (1.4% v/v) was added into each selection well and control well for 15min at RT to elute phage bound on the wells. The phage solution was neutralized by 100µl of 1M Tris-HCl (pH 7.4).
The rescued phage were titrated to estimate the number by sequentially mixing 10µl of the phage into 90µl of LB to 10 6 fold in a low-binding 96-well plate. 10 2 and 10 3 fold infected TG1 were spread onto the other solid selective medium for single colonies.

Phage ELISA
96 individual clones per target in total were picked at 24 candidates per round and inoculated into wells of a 96deep-well plate with 500 µl 2TY containing 100 µg/ml ampicillin and 2% glucose for 3h at 37℃ and 180rpm.
Sufficient volume of M13KO7 helper phage (1µl) were added into each well at 37℃ with shaking overnight.
The plate was centrifugated at 3200rpm for 10min, and the supernatant was transferred to a clean 96-deep-well plate for ELISA. The antigens were used to coat a 96 well plate at 2µg/ml in PBS buffer, 100µl/well, at 4℃      (Salse et al., 2009)       This tree was constructed with protein sequences of these genes by using software of MUSCLE (v3.8.31) (Edgar, 2004) and FastTree (v2.1.10) (Price et al., 2010).
All of the non-N bases were supported by reads with average depth of ~33 times. The blue histogram represents sequencing depth of this region. The green curves represent relationship of pair-end reads.  Tables   Table S1. Summary of the sequencing data from WGS libraries and Hi-C libraries, Related to Figure 1.

Strategy Insert
Size (