Rna seq data

The library sequencing was performed using the Pacbio Sequel platform. The raw data is processed by the software package SMRTlink to obtain the subreads sequence. The circular consensus sequence (CCS) is generated through correction of subreads. According to whether the sequence contains 5′- or 3′-end primer and polyA tail, the sequences are divided into full-length sequences and non-full-length sequences. The full-length sequences were clustered using isoform-level clustering (ICE) to obtain the cluster consensus sequences. The non-full-length sequences were used to correct the consensus sequences to generate high-quality sequences for subsequent analysis. Additional nucleotide errors in consensus reads were corrected using the Illumina RNA-seq data with the software proovread. The GMAP was used to align consensus reads to the cabbage reference genome.

High-quality genomic DNA was extracted from whole blood of a 70-day old male Zhedong goose (A. cygnoides) reared in Xianshan County, Zhejiang Province, China. We constructed 12 paired-end sequencing libraries for whole genome sequencing (WGS) using a WGS kit (Illumina), according to the manufacturer’s recommended protocol. We next sequenced the DNA on a HiSeq 2000 sequencing platform and assembled the short sequences using SOAPdenovo software [44 (link)]. The genome size was calculated from the total length of the sequence reads divided by the sequencing depth. To estimate the sequencing depth, we counted the frequency of each 17-mer from the WGS sequencing reads and plotted the copy-number distribution. The peak value of the frequency curve represented the overall sequencing depth. We used the algorithm, where k_num is the k-mer number, k_depth is the K-mer depth, b_num is the base number, and b_depth is the base depth. G denotes the genome size, and k_depth is the overall depth estimated from the K-mer distribution. To assess the completeness of the assembly, we aligned the unigenes from Illumina RNA-Seq data to the assembled sequence using the BLAT algorithm with default parameters.

A genome browser was prepared using the assembled genome sequences using the JavaScript-based Genome Browser (JBrowse) 1.11.6 [35 (link)]. The assembled sequence and gene models can be found at http://marinegenomics.oist.jp/habu/.
Alignment of amino acid and nucleotide sequences of each protein family was first performed using ClustalW online at http://clustalw.ddbj.nig.ac.jp. Then the alignment of nucleotide sequences was made based upon the amino acid sequence alignment. We carefully and repeatedly confirmed the correspondence between nucleotide sequences of the genome and mRNAs and amino acid sequences of the predicted proteins.
Phylogenetic trees were constructed by the maximum-likelihood method using IQ-TREE (http://www.iqtree.org) [36 (link)], based on aligned nucleotide sequences. Numbers on branches are bootstrap values with 1000 × resampling. The optimal evolutionary model for each phylogenetic tree was selected using ModelFinder [37 (link)] implemented in IQ-TREE.
The expression level of transcript variants in venom gland was analyzed by Pertea’s protocols using open source software tools [38 (link)], HISAT (hierarchical indexing for spliced alignment of transcripts. Ver. 2.1.0) [39 (link)], StringTie (Ver. 1.3.4d) [40 (link)], and Ballgown (Ver. 2.16.0) [41 ] with Illumina RNA-Seq data.

SMRTlink 6.0 software was used to process raw sequence data, and circular consensus sequences (CCS) were obtained by subread BAM files. The CCS BAM files were then classified as full-length non-chimera (FLNC) or non-full length (nFL) according to 5′-primer, 3′-primer and poly-A. Consensus sequences were obtained by clustering full-length sequences with isoform-level clustering (ICE), and additional nucleotide errors in consensus reads were corrected using the Illumina RNA-seq data using the LoRDEC software (Salmela and Rivals, 2014 (link)). The consensus sequences were corrected with non-full-length and non-chimeric sequences, resulting in polished consensus sequences (Fu et al., 2012 (link)). Finally, any redundancy in corrected consensus reads was removed by CD-HIT to obtain transcripts for subsequent analyses and as reference sequences (Fu et al., 2012 (link)). The RSEM software was used to estimate gene expression levels for each sample by mapping the transcript sequences to obtain the read count of each transcript (Li and Dewey, 2011 (link)). Considering the sequence depth and gene length of fragments, the gene expression level was calculated with the fragments per kilobase million (FPKM) method.

For analysis of the P. dicentrarchi genome and transcriptome, trophonts (10⁷) were concentrated by centrifugation, frozen in liquid nitrogen and sent on dry ice to Future Genomic Technologies (Leiden, Netherlands). For sequencing of the complete genome of the ciliate, a combination of short reading sequencing (Illumina technology) and long reading sequencing (Nanopore technology) was used (Oxford Nanopore Technologies). For de novo assembly of the parasite genome, the data sets were combined using the TULIP program v0.4^{40 (link)}. Transcript sequences from Illumina RNA- Seq data (fragments of approximately 100 bp), obtained by amplification by SBS, were assembled using Trinity software (v2.6.5)^{41 (link)}, included within the Galaxy application (https://usegalaxy.org/). The assembled sequences were analyzed by homology, with Blastgo 5.0 software (Biobam, Spain), and annotated. The sequences that encode proteins that are potentially related to the ciliate AOX were then selected from the Tetrahymena thermophila gene and protein sequences database by using the BLASTx tool Wiki TGD (https://www.ciliate.org/blast/blast_link.cgi).

Additional nucleotide errors in consensus reads were corrected using Illumina RNA-seq data. Consensus reads were aligned to the reference genome using a genomic mapping and alignment program (GMAP) [38 (link)]. Reference genome and gene model annotation files were downloaded directly from the genome website (http://fishomics.gooalgene.com/#/download, accessed on 23 September 2021). We selected Hisat2 (v2.1.0) as the mapping tool. Gene structure analysis was performed using the Transcriptome Analysis Pipeline from Isoform Sequencing (TAPIS) [39 (link)]. Gene transcript functions were annotated based on the following databases: NR (NCBI non-redundant protein sequences), NT (NCBI non-redundant nucleotide sequences), Pfam (Protein family), KOG/COG, Swiss-Prot (a manually annotated and reviewed protein sequence database), KO (KEGG Ortholog database), and GO. We used the BLAST software and set the e-value to “1 × 10⁻¹⁰” in NT database analysis [40 (link)]. For NR, KOG, Swiss-Prot, and KEGG database analyses, we used the Diamond BLASTX software and set the e-value to “1 × 10⁻¹⁰.” The Hmmscan software was used for Pfam database analysis [41 (link)].

Salmon (v0.9.1)^{48 (link)} generated transcript expression estimates for each library by mapping the raw Illumina RNA-seq data to the D. rerio transcriptome (vGRCz11)^{49 (link)}. Gene expression estimates were produced by combining all transcript estimates from the same gene⁵⁰ . Genes with low expression (an average log2 transformed count per million less or equal to one) were removed from further analysis^{51 (link)}. The gene expression estimates were normalized based on the sequencing depth of each library. EdgeR (v3.12.1)^{52 (link)} was used to find differentially expressed genes. Generalized linear models were used to test for differential expression based on contrasting the acrylamide treatment against the control. Genes with an FDR adjusted p-value ≤ 0.05 were considered differentially expressed. The sequencing data have been archived in the NCBI Short Read Archive with the BioProject accession number PRJNA515927.