Single-molecule real-time (SMRT) sequencing facilitates Tachypleus tridentatus genome annotation

Journal Pre-proof Single-molecule real-time (SMRT) sequencing Tachypleus tridentatus genome annotation

facilitates

Fangrui Lou, Na Song, Zhiqiang Han, Tianxiang Gao PII:

S0141-8130(19)40237-7

DOI:

https://doi.org/10.1016/j.ijbiomac.2020.01.029

Reference:

BIOMAC 14333

To appear in:

International Journal of Biological Macromolecules

Received date:

12 December 2019

Revised date:

4 January 2020

Accepted date:

4 January 2020

Please cite this article as: F. Lou, N. Song, Z. Han, et al., Single-molecule real-time (SMRT) sequencing facilitates Tachypleus tridentatus genome annotation, International Journal of Biological Macromolecules(2020), https://doi.org/10.1016/ j.ijbiomac.2020.01.029

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier.

Journal Pre-proof

Single-molecule real-time (SMRT) sequencing facilitates Tachypleus tridentatus genome annotation Fangrui Lou1,2, Na Song1, Zhiqiang Han2*, Tianxiang Gao2* 1

Fishery College, Ocean University of China, Qingdao, Shandong, 266003,

China. Fishery College, Zhejiang Ocean University, Zhoushan, Zhejiang, 316022,

of

2

ro

China.

*corresponding author(s): Correspondence and requests for materials should be to

Z.Q.H

(email:

T.X.G.

(email:

lP na

Abstract

and

re

[email protected]).

[email protected])

-p

addressed

Tachypleus tridentatus is a keystone species in marine ecosystems. Its hemolymph

Jo ur

also provides the limulus amebocyte lysate (LAL) for detection of bacterial endotoxin in human medical service. Here we combined SMRT sequencing and Illumina RNA-seq to characterize the novel isoforms, novel genetic loci, fusion isoforms formation and transcriptome structure and further to unveil the transcriptome complexity of T. tridentatus. We identified 26,705 non-redundancy isoforms form 10,919 genetic loci, including 25,713 novel isoforms, 2,403 novel genes and 170 fusion isoforms. In addition, 1,578 novel genes and 23,172 novel isoforms were annotated in the NR, Pfam, KOG, COG, eggNOG, Swiss-Prot, KEGG and GO databases. Meanwhile, we have obtained 4,671 gene family clustering based on

Journal Pre-proof

genetic loci. Furthermore, there are 17,296, 4,887, 1,054, and 1,435 APAs, AS events, lncRNAs, and TFs were identified in the T. tridentatus long-read transcriptome and the target genes of 1,054 lncRNA sequences were also predicted. Overall, our work firstly provided the long-read transcriptome and these data are very necessary to improve the annotation information of T. tridentatus genome and optimize the

of

boundaries of 12,342 original reference annotated genes. Furthermore, these information are a potential resource to study LAL secretion mechanisms in T.

ro

tridentatus.

Jo ur

na

lP

re

-p

Key words: Tachypleus tridentatus; long-read transcriptome; SMRT sequencing

Journal Pre-proof

1. Introduction The marine living fossil, horseshoe crab is a keystone species in marine ecosystems, it can serve as a bioturbator and its eggs is a fundamental protein source for many marine animals [1,2]. Additionally, its hemolymph provides the limulus amebocyte lysate (LAL) for detection of bacterial endotoxin in human medical service

of

[3]. A large number of studies referred horseshoe crab also is an indicator animal for reflecting the environmental state of coastal habitats, this is the case because

ro

hemolymph can help horseshoe crab have particular sensibility and higher tolerance

-p

to metal contaminants compared to other marine animals [4,5]. Therefore, the

re

horseshoe crab is valued because of both its ecological and economic importance. As

lP

an extant species of the horseshoe crab, Chinese horseshoe crab (Tachypleus

na

tridentatus) has maintained their specific ecological value for several million years in the Indo-Pacific region [6]. However, declining populations of T. tridentatus has

Jo ur

widely been reported in many Asian waters owing to habitat destruction, environmental pollution and overfish for LAL production [4,7,8]. In fact, T. tridentatus has existed in the IUCN Red List of Threatened Species [9]. Additionally, it is worth noting that LAL harvest processes might lead to different mortality rates of male and female T. tridentatus due to LAL primarily comes from female [3,10]. This trend can change the behavioral and physiological effects associated with the breeding activity, especially in the spawning season, and further leading to significantly decrease of T. tridentatus resource and ultimately leading to the significantly decreased of egg abundances [11,12]. Although a previous study has revealed the

Journal Pre-proof

genome information of T. tridentatus, the low accuracy and completeness of genome annotation information have affected the exploration of important biological processes (including LAL secretion) of T. tridentatus. Recently, characterization of the transcriptome provides insight into the analysis of gene function and its regulation mechanism. Highthtought short-read sequencing

of

has revolutionized the transcriptome study and provides a novel method to perfect the genome annotation information [13,14]. However, short-read can result in incorrect

ro

annotation information due to short transcripts often fails to identify alternative

-p

splicing, alternative transcription initiation, or alternative transcription termination

re

sites [15,16]. Additionally, whole-genome duplication event also undermined the gene

lP

annotation ability of short transcripts, this is the case because highly similar

na

transcripts were appeared in that event [17]. Above all, the genome of several organisms is very difficult to obtain complete and accurate annotation information

Jo ur

without full-length (FL) sequences. Recently, the third generation sequencing technology Pacific Bioscience (PacBio) provides the opportunity to directly obtain FL transcripts due to its technology accomplishes single molecule real-time (SMRT) sequencing with a read length up to 20 kb [18]. Additionally, SMRT sequencing can comprehensively identify the novel transcribed regions and disambiguate between distinct but highly similar isoforms of functional genes, ultimately it provides abilities to investigate the polymorphism of protein structure and function [19,20]. While higher sequencing error rate existed in the early SMRT sequencing, self-correction of FL transcripts and correction of short reads have reduced the effective error rate

Journal Pre-proof

[21,22]. Therefore, SMRT sequencing can serve as a valuable tool used in the improvement process of the genome annotation [23,24]. To date, SMRT sequencing has been successfully used to analyze the FL transcripts in many marine animals, such as Haliotis discus hannai [25], Danio rerio [26]. However, no studies have been done in T. tridentatus.

of

In this study, SMRT sequencing was used to produce FL transcripts of T. tridentatus and the sequencing process was based on PacBio RSII platform.We can

ro

optimize the T. tridentatus genome annotation information by the FL transcripts and

-p

ultimately provide fundamental resource for some important biological regulatory

re

mechanism of T. tridentatus, such as LAL secretion.

na

2.1 Animal specimens

lP

2. Materials and methods

T. tridentatus specimens were collected from the coastal water of Zhoushan,

Jo ur

China. Female and male T. tridentatuss were immediately euthanized, and muscle tissues were then extracted using sterilized scissors and forceps. Then, muscle samples were separately snap-frozen in liquid nitrogen and stored at -80°C prior to the subsequent experiments. The total RNA of individual muscle was extracted using a standard Trizol Reagent Kit according to the manufacturer’s protocol. The RNA concentration was measured using NanoDrop 2000 (Thermo). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). 2.2 Illumina RNA-seq library construction and sequencing

Journal Pre-proof

We purified mRNA by depleting rRNA from total RNA using the RNA Purification Beads. The remaining RNA was cleaned three times using Beads Binding Buffer, and eluted RNA was incubated at 94°C for 8 min. First-strand cDNA synthesis was performed by combining the fragmented mRNA with 8 μl of First Strand Synthesis Act D Mix and SuperScript II Reverse Transcriptase and incubated at 25°C for 50 min, dropped to 4°C. For second strand synthesis, we added 5 μl of End Repair

of

Control and 20 μl of Second Strand Marking Master Mix and incubated at 16°C for 60

ro

min. The reaction was then purified using 90 μl of AMPure XP beads. Subsequently,

-p

A-Tailing Control and Ligation Control were applied to A-tailing and adapter ligation

re

of the double stranded cDNA, respectively. In the cDNA fragments enrich reaction

lP

process, we added 5 μl of PCR Primer Cocktail and 25 μl of PCR Master Mix and ran

na

the initial denaturing at 98°C for 30 sec, followed by 15 cycles of denaturing at 98°C for 10 sec, annealing at 60°C for 30 sec, and extending at 72°C for 30 sec, and

Jo ur

finished with a final extending at 72°C for 5 min. RNA-seq libraries were diluted to 10 pM and quantified using an Agilent 2100 Bioanalyzer. Libraries were sequenced on an Illumina HiSeq 2000 across one lane with paired-end 150 bp. 2.3 SMRT sequencing library construction and sequencing The first-strand cDNA was synthesized and amplified using the SMARTerTM PCR cDNA Synthesis Kit (Takara Bio USA, Inc., Mountain View, CA, USA). High-quality large-scale PCR products were purified and eluted to the optimal size (approximately 1-6 Kb) using the BluePippinTM Size-Selection System (Sage Science Inc., Beverly, MA, USA). Subsequently, each SMRTbell template library was constructed after

Journal Pre-proof

repairing the damage, blunt end ligation, and adding the SMRTbell adapter for selected cDNAs. The SMRTbell templates combined with polymerases and then the polymeride was bound to the zero-mode waveguide. In the sequencing process, two SMRT cells were prepared on the PacBio RSII platform. 2.4 Illumina RNA-seq data processing

of

Raw RNA-seq reads in FASTQ format were firstly filtered and then clean reads were obtained by removing reads containing sequencing adapters, ploy-N, and low

ro

quality. The remaining high-quality reads were then mapped to the T. tridentatus

-p

genome sequence using HISAT2 software [27]. Meanwhile, we analyzed the location

re

distribution and cover the depth distribution of RNA-seq data on the different

lP

chromosomes.

na

2.5 SMRT sequencing data processing

All polymerase reads produced by the PacBio RSII platform were filtered by

Jo ur

removing reads with sequences length less than 50bp and accuracy less than 0.80 using SMRT analysis software v2.3 Suite (http:// www.pacb.com/devnet/) with the following parameters: readScore = 0.80 and minLength = 50. The remaining high-quality polymerase reads were broken at the adaptor position. Subsreads were generated by trimming all adaptor sequences of polymerase reads. Effective subreads were further obtained by removing reads with sequences length less than 50bp using SMRT analysis software v2.3 Suite with parameters: minSubReadLength = 50. Circular consensus (CCS) reads were processed using Iso-seq pipeline of SMRTLink software with the following parameters: minFullPass = 0 and minPredictedAccuracy

Journal Pre-proof

= 0.90. Then, these CCS reads containing both polyA tail signal, 5' and 3' primers were considered as FL reads, otherwise, these CCS reads were considered as Non-full-length (NFL) reads. It is worth noting that two cDNA template sequences were always direct connected to generate chimeric reads due to the concentration of adaptors or SMRTbell in the library construction process. These FL reads without

of

chimeric reads were considered as FL non-chimetric (FLNC) reads (Figure 1.A). Next, iterative clustering for error correction (ICE) was used to obtain consensus isoforms

ro

and polished FL consensus sequences from ICE were polished using Quiver [18]. It is

-p

worth noting that polished high-quality consensus sequences were classified with the

re

criteria post-correction accuracy above 99%. In order to improve sequences accuracy,

lP

polished low-quality consensus sequences were corrected using clean female

na

RNA-seq data based on Proovread software [28]. In fact, the 5' end of some FLNC reads will be degraded in the sequencing process, in spite of the 3' end can remain

Jo ur

intact due to the presence of the polyA structure. Therefore, redundant transcripts might be generated in the clustering process of FLNC reads. In order to remove redundant sequences, the calibrated consensus sequences were mapped to the reference genome [29] using Genomic Mapping and Alignment Program (GMAP; [30]) with the following parameters: --cross-species and --allow-close-indels 0. The high percent of identity (PID) and high percent of coverage (PCO) aligned FLNC reads were regarded as the non-redundant sequences and these sequences with low PID (less than 90 ％ ) and low PCO (less than 85 ％ ) were removed by the cDNA_Cupcake

package

(https://github.com/Magdoll/cDNA_Cupcake/wiki).

Journal Pre-proof

Additionally, exon 5' end difference was not considered when collapsing redundant sequences and therefore we incorporated sequences that were only different at the exon 5' end (Figure 1.B). BUSCO software was applied to the quantitative assessment of the completeness of non-redundant sequences [31]. 2.6 Identification of fusion isoforms

of

All consensus sequences were used to find fusion isoforms and the identification criteria were as follows (Figure 1.C; [18]): (i) a single FL consensus sequence must be

ro

mapped to two or more annotation loci in the T. tridentatus genome; (ii) the minimum

-p

coverage for each annotation loci must be at least 5% of FL consensus sequences and

re

the minimum coverage length is 1 bp; (iii) the combined alignment coverage of all

lP

annotation loci must be at least 95% of FL consensus sequences; (iv) distance

na

between two mapped locus must be at least 10 kb.

optimization

Jo ur

2.7 Identification of potential novel genes, novel isoforms and genetic structure

We identified some potential novel genes, novel isoforms based on the mapping information of FL consensus sequences. Specifically, the FL consensus sequence that having a low overlapping ratio (< 20%) with the reference gene site or having a high overlapping ratio (> 20%) but the gene direction was not consistent with reference gene was considered as a novel gene. The filtering criteria of novel isoform were as follows: (i) the final splice site of 3' end changed; (ii) new intron or new exon emerged (Figure 1.D; [18]). All novel genes and novel isoforms were used to analyze the potential functions based on the following databases using the BLAST alignments:

Journal Pre-proof

NR, KOG/COG, eggNOG, Swiss-Prot, KEGG, and GO. Meanwhile, we obtained the gene family clustering based on the Pfam database annotation information. Finally, we also re-optimized the genetic boundaries of the reference annotations based on the

lP

re

-p

ro

of

novel genetic information.

na

Figure 1. PacBio data processing. (A). FLNC. (B). Remove redundancy of consensus

Jo ur

sequences. (C). Fusion isoforms identification. (D). Novel isoforms identification. 2.8 T. tridentatus transcriptome structure analysis FLNCs were applied to identify alternative polyadenylation using TAPIS pipeline [32]. The comprehensive distribution and classification analysis of alternative splicing (AS) structure in non-redundant sequences were identified using Astalavista software [33]. Four computational approaches include coding potential calculator (CPC; [34]), coding-non-coding index (CNCI; [35]), coding potential assessment tool (CPAT; [36]) and Pfam database were combined to sort non-protein coding RNA candidates from putative protein-coding RNAs in the novel isoforms. Putative protein-coding RNAs were filtered out using a minimum length and exon number threshold. Isoforms with

Journal Pre-proof

lengths more than 200 nt and have more than two exons were selected as long non-coding RNA (LncRNA) candidates and further screened using CPC, CNCI, CPAT, and Pfam database that have the power to distinguish the protein-coding genes from the non-coding genes. Additionally, LncRNA target genes were predicted based on the position and effect of LncRNA and mRNA using LncTar software [37]. Finally, iTAK

of

software [38] was used to identify the T. tridentatus transcription factors (TFs) from novel isoforms based on the animalTFDB 2.0 [39] transcription factor database.

ro

2.9 Gender-specific expression analyses

-p

Gene expression levels were obtained by mapping clean RNA-seq reads to the

re

long-read transcripts. Then, differentially expressed genes (DEGs) of female and male

used

as

the

filtering

thresholds

of

DEGs

na

were

lP

T. tridentatus were identified using EBSeq, FDR ＜ 0.01 and Fold Change ≥ 2

(http://www.bioconductor.org/packages/devel/bioc/html/EBSeq.html). Subsequently,

Jo ur

the functional categories and biochemical metabolic pathways of DEGs were predicted based on Blast2GO software, and the threshold parameter of significant difference was Q ≤ 0.05. 3. Results and Discussion 3.1 Sequencing of the T. tridentatus transcriptome using Illumina platform The raw RNA-seq reads of female and male T. tridentatus were released in the NCBI Sequence Read Archive under BioProject number PRJNA576971, with accession number of SRR10270295 and SRR10270295. After filtering, a total of 42,435,306 and 48,654,696 clean RNA-seq reads were obtained in female and male T.

Journal Pre-proof

tridentatus, respectively. Then, we used HISAT2 to sequence align clean RNA-seq reads with the reference genome and obtain the location information of clean reads on the reference genome. In female RNA-seq data, 95.18% of clean reads could be mapped to the T. tridentatus genome, of which 81.71% were mapped to the only loci in the genome. In male RNA-seq data, 94.98% of clean reads could be mapped to the

of

T. tridentatus genome, of which 84.17% were mapped to the only loci in the genome. Meanwhile, we analyzed the coverage depth distribution of reads that have been

ro

mapped to the reference genome (Figure 2). In addition, we counted the number of

-p

successful reads in different regions of the reference genome (exons, introns and

re

intergenic regions). It is undeniable that although a large number of reads are mapped

lP

to the exon region (72.65% in female and 68.91% in male), a small number of reads

na

are mapped to the intron region (11.76% in female and 13.99% in male) and intergenic region (15.59% in female and 17.10% in male). In fact, reads from mature

Jo ur

mRNA should be mapped to the exon region. Additionally, reads mapped to intron region are related to mRNA precursor and intron retention with AS. It is worth noting that the intergenic regions is belonging to the non-coding region, therefore it is not transcribed. In the present study, some RNA-seq reads were mapped to intergenic regions, which may be related to incomplete annotation information of T. tridentatus genome.

Journal Pre-proof

of

Figure 2. The location and coverage depth distribution of mapped reads on the

ro

reference genome. (A). Female T. tridentatus. (B). Male T. tridentatus.

-p

3.2 Sequencing of the T. tridentatus transcriptome using PacBio plotform

re

In the present study, high-quality total RNA was extracted from T. tridentatus

lP

and then sequenced on the Pacific Bioscience RS II platform. After filtering, we obtained 32.73 G effective subread data from 2 SMRT cells and the subreads were

na

released in the NCBI Sequence Read Archive under BioProject number

Jo ur

PRJNA577185, with accession number of SRR10271121 and SRR10271122. Furthermore, we detected 793,144 CCS reads with a mean sequencing depth of 16 passes. Among them, 635,260 (80.09％) were considered as FLNC. ICE was applied to improve the consensus accuracy of FLNCs and a total of 38,762 consensus isoforms were obtained, with the mean length of 3,124bp. Furthermore, 37,756 high-quality consensus isoforms were mapped to T. tridentatus genome using GMAP and 1,001 low-quality consensus isoforms were corrected by short read RNA-seq data. Finally, 26,705 non-redundancy isoforms were yielded and used to the following analyses. The length distribution of CCS, FLNC, and consensus isoform is shown in Supplementary file 1. We also provided a preliminary assessment of the T. tridentatus

Journal Pre-proof

annotation representation in the non-redundancy isoforms based on the aligned results in the reference genome (Figure 3). In reference genome, we found 60,679 isoforms covering 37,712 genetic loci and 96.54% (58,578/60,679) isoforms were > 1kb in length. Out of 37,712 reference genetic loci, 28,494 (75.56%) genetic loci produced more than one isoforms, corresponding to 51,461 isoforms. Our long-reads transcript

of

mapping results showed that 26,705 non-redundancy isoforms were allocated to 10,919 genetic loci and about 84.31% (9,206/10,919) genes were > 1kb in length.

ro

Among 26,705 isoforms, 2,030 (7.60%) were single-exon isoforms and 24,675

-p

(92.40%) were multiple-exon isoforms. Additionally, 6,029 (55.22%) genetic loci

re

produced more than one isoforms, corresponding to 21,745 isoforms. We determined

lP

whether the number of isoform was affected by the gene length and results showed

file

2).

na

that the number of isoform was not correlated with the gene length (Supplementary Interestingly,

we

found

that

the

gene

SLC7A8

(gene

ID:

Jo ur

evm.TU.Hic_chr_13.3049) had the largest number of isoforms in both reference genome and Pac-Bio data, with 72 and 71 isoforms, respectively, which played a vital role in the amino acid transport process of T. tridentatus [40].

lP

re

-p

ro

of

Journal Pre-proof

na

Figure 3. Distributed visualization of genome locus, transcriptome locus, alternative

Jo ur

splicing, fusion transcripts and LncRNA at chromosome level. The curved black line in each circle corresponding to each chromosome. The lines above and below the each curved black line represent the locus information of the sense strand (+) and the antisense strand (-), respectively.

Finally, 86.1% complete BUSCOs in the arthropoda_odb9 database were covered by isoforms and therefore the completeness of T. tridentatus long-read transcriptome was be confirmed (Figure 4).

ro

of

Journal Pre-proof

-p

Figure 4. BUSCO assessment results.

re

3.3 Fusion isoforms

lP

A fusion isoform is a chimeric RNA encoded by one fusion gene or two different

na

genes that are subsequently joined by trans-splicing or somatic chromosomal rearrangements. Gene fusion is a common feature in many eukaryotes, but no

Jo ur

examples of gene fusion have been described in T. tridentatus. In the present study, we identified 170 fusion isoforms covered 1,943 isoforms and were corresponded to 340 genetic loci. Distribution of fusion isoforms revealed that large amounts of fusion isoforms were predominately located in chr8 (31 fusion isoforms), followed by chr13 (27 fusion isoforms) and chr3 (22 fusion isoforms). Perhaps for two reasons: (1) The quantity variance of fusion isoforms may be influenced by the chromosome length. In fact, we found that chr13 has the longest length (275,858,341 bp) based on reference genome [29]; (2) Chromosomal rearrangements (including chromosomal translocation, interstitial deletion and inversion) caused by transposons eventually lead to gene

Journal Pre-proof

fusion, therefore, we speculate that there may be more transposons on chr3, chr8 and chr13 [41]. However, it cannot be denied that the gene fusion identified by the transcriptome sequencing may come from the RNA fusion after transcription of different genes, rather than the actual genomic variation. Future studies still need to address chromosome-specific gene fusion. The pervious study considered that fusion

of

isoforms usually formed by the concatenation of two or more separate isoforms,

ro

whereas all fusion isoforms in this study from two isoforms. 3.4 Novel genes and novel isoforms

-p

3.4.1 Identification and annotation of novel genes and novel isoforms

re

Mapping information showed that all non-redundancy isoforms were classified

lP

into three groups: (i) 992 known isoforms from 567 genetic loci, (ii) 21,242 novel

na

isoforms of 8,516 genetic loci and (iii) 4,471 isoforms were probably from 2,403 novel genetic loci that do not overlap with any annotated loci. It is undeniable that our

Jo ur

non-redundant isoforms only covered about 24.09% (9,083/37,712) of reference annotation loci. A recent zebrafish long-read transcriptome study reported a similarly low proportion (25.4％) of transcripts accurately matched to the reference genome [26]. This may be related to the completeness of genomic annotation information. The low coverage of the non-redundant isoforms found in the genome data also supports the presence of significant previously undiscovered transcriptional diversity [26]. It is worth noting that a large number of new isoforms (25,713) may include transcripts with skipped exons, retained introns, and other types of canonical alternative splicing events [42]. Additionally, two other categories of transcripts were classified as new

Journal Pre-proof

isoforms: (1) these long-read transcripts with 5'- or 3'-end exons missing from the reference transcript; (2) these long-read transcripts mapped to regions that are generally considered less likely to give rise to true transcripts. It is worth noting that the annotation information of T. tridentatus genome may be limited by the previous data used, therefore the potential novel genes and novel

of

isoforms were further to annotation to improve the completeness of T. tridentatus genome annotation information. We searched 2,403 novel genes and 25,713 novel

ro

isoforms in NR, Pfam, KOG, COG, eggNOG, Swiss-Prot, KEGG and GO databases.

-p

In total, 1,578 novel genes and 23,172 novel isoforms were annotated in the

re

aforementioned databases. Specifically, 1,550, 1,098, 973, 352, 1,419, 814, 665 and

lP

607 novel genes could be found in NR, Pfam, KOG, COG, eggNOG, Swiss-Prot,

na

KEGG and GO database, respectively. Additionally, 23,056, 10,051, 18,389, 8,655, 22,158, 16,028, 13,765 and 12,497 novel isoforms could be found in NR, Pfam, KOG,

Jo ur

COG, eggNOG, Swiss-Prot, KEGG and GO database, respectively. We further integrated all annotated novel genes and novel isoforms to the T. tridentatus reference genome annotation and ultimately generated an updated reference genome annotation information that containing 54,450 annotated isoforms (Supplementary file 3) and 32,755 annotated genetic loci (Supplementary file 4). Meanwhile, we have obtained 4,671 gene family clustering based on genetic loci (Supplementary file 5). 3.4.2 Genetic structure optimization It is worth noting that the reference annotation is often not accurate due to the limitations of the software or data, so we optimize the gene structure of the T.

Journal Pre-proof

tridentatus reference annotation based on the novel genes and novel isoforms. In the present study, the map read sequence appeared outside the boundaries of 12,342 original reference annotated genes, and the untranslated region (UTR) upstream or downstream of these reference annotated genes was further extended and the modification of the gene boundaries was finally completed (Supplementary file 6).

of

3.5 T. tridentatus transcriptome structure 3.5.1 Alternative polyadenylation (APA)

ro

Polyadenylation of most mRNAs 3' end is necessary for transport, localization,

-p

stability and translation of RNA in the eukaryotes protein biosynthesis. Previous

re

studies have shown that polyadenylation can generate different transcripts in the

lP

coding region or the 3' end non-coding region, thus complicating the eukaryotic

na

transcriptome and generating different regulatory mechanisms. Additionally, differential polyadenylation of mRNAs may play different roles in T. tridentatus

Jo ur

development. Here we reported the first genome-wide map of global polyadenylation sites of T. tridentatus. A total of 17,296 APAs from 8,470 genes were obtained from T. tridentatus transcriptome (Supplementary file 7). It is worth noting that the establishment of accurate polyadenylate cleavage site also can help us to determine the 3' end location of each gene in T. tridentatus, and finally improve gene annotation. 3.5.2 Alternative splicing (AS) AS is an important transcriptional mechanism that regulates the expression of eukaryotic traits. Or more specifically, the same gene can selectively recognize exons and splicing sites to produce mature mRNA and protein isoforms with different

Journal Pre-proof

functions during splicing. In the present study, we quantified the number of AS events and compared the distribution of AS types in the T. tridentatus transcriptome based on Astalavista software (Figure 5). In total, we have discovered 4,887 AS events in the T. tridentatus transcriptome and the main AS event types were intron retention, exon skipping, alternative 3' splice site, alternative 5' splice site and mutually exclusive

of

exon (Supplementary file 8). Among all AS events, exon skipping predominated, accounting for 41.80% of all AS events, while mutually exclusive exon was the least

ro

frequent (2.82%). The study of AS may be very important to understanding why there

Jo ur

na

lP

re

-p

are differences in the LAL secretion of different genders T. tridentatus.

Figure 5. The types and number of AS events in the T. tridentatus transcriptome. 3.5.3 Long non-coding RNA (LncRNA) and LncRNA target genes LncRNA is a non-coding RNAs with a length of more than 200 nucleotides. LncRNA plays an important role in many life activities such as dosage compensation effect, epigenetic regulation, cell cycle regulation and cell differentiation regulation, and thus LncRNA has become a research hotspot. Therefore, LncRNAs might constitute an important class of regulators of LAL secretion in T. tridentatus. In the present study, we used four computational approaches (CPC, CNCI, CPAT and Pfam

Journal Pre-proof

database) to identify putative lncRNAs from all novel isoforms. By filtering and excluding isoforms with lengths more than 200nt, 1,486, 7,737, 5,015 and 5,010 lncRNAs were found based on CPC, CNCI, CPAT and Pfam database, respectively. Only 1,054 lncRNAs were found in all four computational approaches (Figure 6). According to the position of lncRNA on the annotation information of reference

of

genome, 1,054 lncRNAs were divided into four categories: lincRNA (575), antisense-lncRNA (50), intronic-lncRNA (160) and sense-lncRNA (244). Additionally,

ro

we used two methods to predict the target genes of 1,054 lncRNA sequences and the

Jo ur

na

lP

re

-p

predicted gene results are shown in the Supplementary file 9.

Figure 6. The number and classification of lncRNAs.

3.5.4 Transcription factor (TF) TF is a kind of protein that can bind to specific nucleotide sequence on the upstream of a gene. It can regulate the binding of RNA polymerase to DNA template and finally influence the transcription of genes. Therefore, TFs may be very important in regulating the LAL secretion in T. tridentatus. In the present study, we predicted 1,435 TFs from the novel isoforms of T. tridentatus based on the animal transcription

Journal Pre-proof

factor database (animalTFDB 2.0). Furthermore, all TFs were classified into 47 families and a large number TFs were dominant in zf-C2H2, bHLH, HMG families Supplementary file 10. 3.6 Gender-specific expression analyses of T. tridentatus To elucidate the gender-specific expression pattern of T. tridentatus, we first

of

counted the number of DEGs. Results showed that a total of 1,796 DEGs were obtained between female and male T. tridentatus (Figure 7.A). Among all DEGs, 964

ro

up- and 832 down-regulated genes were founded in female T. tridentatus. To explore

-p

the potential functions and regulatory pathways of these differential genes, we

re

conducted GO terms and KEGG pathways enrichment analysis for 1,796 DEGs.

lP

Results showed that a total of 698 DEGs were successfully mapped to 47 GO terms

na

(Figure 7.B). Among these GO terms, “metabolic process”, “catalytic activity” and “cell” were dominant in “biological process”, “molecular function” and “cellular

Jo ur

component”, respectively. Based on KEGG database, the KEGG classifications of DEGs were shown in Figure 7.C. In brief, some metabolism, immune, genetic information

processing,

environmental

information

processing

and

cellular

processes-related pathways that may participate in the sexual differentiation process were detected in this study. Interestingly, our results indicated that 14 genes are enriched in the Wnt signaling pathway. In fact, the Wnt signaling pathway is a conservative signaling network, which involve in a variety of biological processes, such as embryonic development, growth regulation and others [43]. It is worth noting that Wnt4 gene in the Wnt signaling pathway is regarded as a sex-determining gene

Journal Pre-proof

and it is crucial for female sexual development [44]. This is the case because Wnt4 gene can regulate the formation of the mullerian duct and then inhibit substance and testosterone secretes, which ultimately lead to sexual differentiation [45]. Further research is required to elucidate the physiological function of these DEGs, and ultimately provide fundamental resource for understanding some critical biological

Jo ur

na

lP

re

-p

ro

of

processes of T. tridentatus, such as LAL secretion.

Figure 7. The number (A), GO term (B) and KEGG pathway (C) analyses of gender-specific DEGs. 4. Conclusions In the present study, SMRT sequencing and RNA-Seq were applied to sequencing the T. tridentatus transcriptome. We obtained 26,705 non-redundancy

Journal Pre-proof

isoforms, including fusion isoforms, novel isoforms of previously annotated genes as well as novel genetic loci. We further integrated all annotated novel genes and novel isoforms to the T. tridentatus reference genome annotation and ultimately generated an updated reference sequence database. Additionally, we identified APA, AS, lncRNA and TF in T. tridentatus transcriptome data and predicted the target genes of

tridentatus

of

lncRNAs. These data are very necessary to improve the annotation information of T. genome and optimize the gene structure.

Furthermore, these

-p

LAL secretion mechanisms in T. tridentatus.

ro

non-redundancy isoforms and transcriptome structure are a potential resource to study

re

Acknowledgments

Reference

na

China (2017YFA0604902).

lP

This work was supported by the National Key Research and Development Program of

Jo ur

[1] D.M. Rudkin, G.A. Young, Horsehsoe crabs-An ancient ancestry revealed. In: Tanacredi, J.T., Botton, M.L., Smith, D.R. (eds) Biology and Conservation of Horseshoe Crabs. Springer, US 2009 25-44. [2] M.L. Botton, C.N. Shuster, Horseshoe crabs in a food web: who eats whom? In: Shuster, C.N. Jr, Barlow, R.B., Brockmann, H.J. (ed) The American Horseshoe Crab. Harvard Press, Cambridge 2003 133-153. [3] R.L. Anderson, W.H. Watson, C.C. Chabrt, Sublethal behavioral and physiological effects of the biomedical bleeding process on the American horseshoe crab, Limulus polyphemus. Biol. Bull. 225 (2013) 137-151. DOI: 10.1086/BBLv225n3p137.

Journal Pre-proof

[4] C.P. Chen, H.Y. Yeh, P.F. Lin, Conservation of the horseshoe crab at Kinmen, Taiwan: strategies and practices. Biodivers. Conserv. 13 (2004) 1889-1904. DOI: 10.1023/b:bioc.0000035868.11083.84 [5] M.L. Botton, M. Hodge, T.I. Gonzalez, High tolerance to tributyltin in embryos and larvae of the horseshoe crab, Limulus polyphemus. Estuaries 21 (1998) 340-346.

of

DOI: 10.2307/1352480. [6] K. Sekiguchi, H. Sugita, Systematics and hybridization in the four living species

ro

of horseshoe crabs. Evolution 34 (1980) 712. DOI: 10.2307/2408025

-p

[7] Y. Liao, S. Hong, X. Li, A survey on the horseshoe crabs in the north of South

re

China Sea. Acta. Zool. Sin. 47 (2011) 108-111 (Chinese with English abstract).

lP

[8] P.K.S. Shin, H. Li, S.G. Cheung, Horseshoe crabs in Hong Kong: current

na

population status and human exploitation. In: Tanacredi, J.T., Botton, M.L., Smith,

347-360.

Jo ur

D.R. (eds) Biology and Conservation of Horseshoe Crabs. Springer, US, 2009

[9] Centre, World Conservation Monitoring. Tachypleus tridentatus. The IUCN Red List of Threatened Species 1996. [10] D. Rutecki, R.H. Carmichael, I. Valiela, Magnitude of harvest of Atlantic horseshoe crabs, Limulus polyphemus, in Pleasant Bay, Massachusetts. Estuaries 27 (2004) 179-187. DOI: 10.2307/1353450. [11] L. Hurton, J. Berkson, Potential causers of mortality for horseshoe crabs (Limulus polyphemus) during the biomedical bleeding process. Fish. Bull. 104 (2006) 293-298. DOI: 10.1111/j.1444-2906.2006.01172.x.

Journal Pre-proof

[12] A.S. Leschen, S.J. Correia. Mortality in female horseshoe crabs (Limulus polyphemus) from biomedical bleeding and handling: implications for fisheries management.

Mar.

Freshw.

Behav.

Physiol.

43

(2010)

135-147.

DOI:

10.1080/10236241003786873. [13] F. Ozsolak, P.M. Milos, RNA sequencing: advances, challenges and opportunities. Nat. Rev. Genet. 12: 87-98. DOI: 10.1038/nrg2934.

of

[14] J. Duitama, P.K. Srivastava, I.I. Măndoiu, Towards accurate detection and

ro

genotyping of expressed variants from whole transcriptome sequencing data. BMC

-p

Genomics 13 (2012) S6. DOI: 10.1186/1471-2164-13-s2-s6.

re

[15] T. Steijger, J.F. Abril, P.G. Engstrom, et al., Assessment of transcript

lP

reconstruction methods for RNA-seq. Nat. Methods 10 (2013) 1177-1184. DOI:

na

10.1038/nmeth.2714.

[16] H. Tilgner, D. Raha, L. Habegger, L., et al., Accurate identification and analysis

Jo ur

of human mRNA isoforms using deep long read sequencing. G3 3 (2013) 387-397. DOI: 10.1534/g3.112.004812.

[17] J.H. Postlethwait, Y.L. Yan, M.A. Gates, et al., Vertebrate genome evolution and the zebrafish gene map. Nat. Genet. 18 (1998) 345-349. DOI: 10.1038/ng0498-345. [18] Y. Li, C.C. Fang, Y.H. Fu, et al., A survey of transcriptome complexity in Sus scrofa using single-molecule long-read sequencing. DNA Res. 25 (2018) 421-437. DOI: 10.1093/dnares/dsy014 [19] M.J. Chaisson, J. Huddleston, M.Y. Dennis, et al., Resolving the complexity of the human genome using single-molecule sequencing. Nature 517 (2014) 608-611.

Journal Pre-proof

DOI: 10.1038/nature13907. [20] F. Min, S. Wang, L. Zhang, Survey of programs used to detect alternative splicing isoforms from deep sequencing data In Silico. Biomed. Res. Int. 2015 (2015) 831352. DOI: 10.1155/2015/831352. [21] K.F. Au, V. Sebastiano, P.T. Afshar, et al., Characterization of the human ESC

of

transcriptome by hybrid sequencing. Proc. Natl Acad. Sci. USA 110 (2013) E4821-E4830. DOI: 10.1073/pnas.1320101110.

ro

[22] S.S. Vembar, M. Seetin, C. Lambert, et al., Complete telomere-to-telomere de

-p

novo assembly of the Plasmodium falciparum genome through long-read (>11 kb),

re

single molecule, real-time sequencing, DNA Res. 23 (2016) 339-351. DOI:

lP

10.1093/dnares/dsw022.

na

[23] K.E. Kim, P. Peluso, P. Babayan, P., et al., Long-read, whole-genome shotgun sequence data for five model organisms. Sci. Data 1 (2014) 140045. DOI:

Jo ur

10.1038/sdata.2014.45.

[24] S.E. Abdel-Ghany, M. Hamilton, J.L. Jacobi, et al., A survey of the sorghum transcriptome using single-molecule long reads. Nat. Commun. 7 (2016) 11706. DOI: 10.1038/ncomms11706. [25] M.A. Kim, J.S. Rhee, T.H. Kim, et al., Alternative splicing profile and sex-preferential gene expression in the female and male Pacific Abalone Haliotis discus hannai. Genes 8 (2017) 99. DOI: 10.3390/genes8030099. [26] G. Nudelman, A. Frasca, B. Kent, High resolution annotation of zebrafish transcriptome using long-read sequencing. Genome res. 28 (2018) 1415-1425. DOI:

Journal Pre-proof

10.1101/gr.223586.117. [27] D. Kim, B. Langmead, S.L. Salzberg, HISAT: a fast spliced aligner with low memory requirements. Nat Methods 12(2015) 357-360. DOI: 10.1038/nmeth.3317. [28] T. Hackl, R. Hedrich, J. Schultz, et al., Proovread: large-scale high-accuracy PacBio correction through iterative short read consensus. Bioinformatics 30 (2014)

of

3004-3011. DOI: 10.1093/bioinformatics/btu392. [29] L. Gong, G.Y. Fan, Y.D. Ren, et al., Chromosomal level reference genome of

ro

Tachypleus tridentatus provides insights into evolution and adaptation of horseshoe

-p

crabs. Mol. Ecol. Resour. 19 (2019) 744-756. DOI: 10.1111/1755-0998.12988.

re

[30] T.D. Wu, C.K. Watanabe, GMAP: a genomic mapping and alignment program

na

10.1093/bioinformatics/bti310.

lP

for mRNA and EST sequences. Bioinformatics 21 (2005) 1859-1875. DOI:

[31] F.A. Simao, R.M. Waterhouse, P. Ioannidis, et al., BUSCO: assessing genome

Jo ur

assembly and annotation completeness with single-copy orthologs. Bioinformatics 19 (2015) 3210-3212. DOI: 10.1093/bioinformatics/btv351. [32] S.E. Abdelghany, M. Hamilton, J.L. Jacobi, et al., A survey of the sorghum transcriptome using single-molecule long reads. Nat. Commun. 7 (2016) 11706. DOI: 10.1038/ncomms11706. [33] S. Foissac, M. Sammeth, ASTALAVISTA: dynamic and flexible analysis of alternative splicing events in custom gene datasets. Nucleic Acids Res. 35 (2007) 297-299. DOI: 10.1093/nar/gkm311. [34] L. Kong, Y. Zhang, Z.Q. Ye, et al., CPC: assess the protein-coding potential of

Journal Pre-proof

transcripts using sequence features and support vector machine. Nucleic Acids Res. 36 (2007) 345-349. DOI: 10.1093/nar/gkm391. [35] L. Sun, H.T. Luo, D.C. Bu, et al., Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41 (2013) e166. DOI: 10.1093/nar/gkt646.

of

[36] L. Wang, H.J. Park, S. Dasari, et al., CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 6 (2013) e74.

ro

DOI: 10.1093/nar/gkt006.

-p

[37] J. Li, W. Ma, P. Zeng, et al., LncTar: a tool for predicting the RNA targets of long

re

noncoding RNAs. Brief. Bioinform. 16 (2015) 806-812. DOI: 10.1093/bib/bbu048.

lP

[38] Y. Zheng, C. Jiao, H. Sun, et al., iTAK: a program for genome-wide prediction

na

and classification of plant transcription factors, transcriptional regulators, and protein kinases. Mol. Plant 9 (2016) 1667-1670. DOI: 10.1016/j.molp.2016.09.014.

Jo ur

[39] H.M. Zhang, T. Liu, C.J. Liu, et al., AnimalTFDB 2.0: a resource for expression, prediction and functional study of animal transcription factors. Nucl. Acids Res. 43 (2015): 76-81. DOI: 10.1093/nar/gku887.[40] M.T. Bassi, M.P. Sperandeo, B. Incerti, et al., SLC7A8, a gene mapping within the lysinuric protein intolerance critical region, encodes a new member of the glycoprotein-associated amino acid transporter family. Genomics 62 (1999): 297-303. DOI: 10.1006/geno.1999.5978. [41] F. Mertens, B. Johansson, T. Fioretos, et al., The emerging complexity of gene fusions in cancer. Nat. Rev. Cancer. 15 (2015):371-81. DOI: 10.1038/nrc3947. [42] Q. Pan, O. Shai, L.J. Lee, et al., Deep surveying of alternative splicing

Journal Pre-proof

complexity in the human transcriptome by high-throughput sequencing. Nat. Genet. 40 (2008): 1413-1415. DOI: 10.1038/ng.259. [43] J. Gao, X.W. Wang, Z.H. Zou, et al., Transcriptome analysis of the differences in gene expression between testis and ovary in green mud crab (Scylla paramamosain). BMC Genomics 15 (2014): 585. DOI: 10.1186/1471-2164-15-585

of

[44] S. Vainio, M. Heikkila, A. Kispert, A., et al., Female development in mammals is regulated by Wnt-4 signalling. Nature 397(1999): 405-409. DOI: 10.1038/17068.

ro

[45] F.R. Lou, T.Y. Yang, Z.Q. Han, et al., Transcriptome analysis for identification of

-p

candidate genes related to sex determination and growth in Charybdis japonica. Gene

Jo ur

na

lP

re

66 (2018): 10-16. DOI: 10.1016/j.gene.2018.07.044.

Journal Pre-proof Author statement

Jo ur

na

lP

re

-p

ro

of

Lou Fangrui: Methodology, Validation, Formal analysis, Writing-Original draft. Song Na: Validation, Supervision, Project administration. Han Zhiqiang: Conceptualization, Methodology, Resources, Writing - Review & Editing, Funding acquisition. Gao Tianxiang: Validation, Supervision, Writing- Reviewing & Editing,