De novo transcriptome assembly and annotation for gene discovery in the streamside salamander (Ambystoma barbouri)

Tissue collection

Eggs were collected from Wilson County and Rutherford County to represent the southern and central Tennessee ESUs as well as from populations in Preble County, Ohio and Scott and Boone counties within Kentucky (Hubbs et al., 2022). All eggs were brought back to Tennessee Technological University to be staged and separated into treatment groups as part of an incubation common garden experiment exploring differential gene expression in response to chronic temperature stress (Gaupp et al., in prep). Staging consisted of placing each individual egg under a dissecting microscope and characterizing its Harrison stage of development (Duellman and Trueb, 1994). Upon initial staging, eggs ranged from stage 1 to stage 13 (according to descriptions in Schreckenberg and Jackson 1975). Individuals from each clutch were separated randomly between two thermal treatments: 10 °C and 20 °C to simulate both optimal and thermally stressful environments (Thulander and Hall 2025). Eggs were incubated at these temperatures until reaching Harrison developmental stage 40, the stage at which it is possible to hatch for most Ambystoma species. Upon reaching this developmental stage, individuals were euthanized via immersion in 1 g/L solution of Tricaine MS-222 buffered to pH 7 with sodium bicarbonate, as recommended by the American Society of Ichthyologists and Herpetologists (Beaupre et al., 2004) and the American Veterinary Medical Association, until cessation of a visibly detectable heartbeat. Samples were immediately preserved in RNALater, and stored at −80 °C following a 24-hour incubation at 4 °C.

In addition to individuals collected for the incubation common garden experiment, heart, liver, spleen, and tail tissues were collected from individuals of about 1 year of age. These individuals were collected by Thulander and Hall (2025) as eggs and kept at Tennessee Tech for continued observation; however, many individuals developed an unidentified skin condition. Individuals were euthanized via immersion in 1 g/L solution of Tricaine MS-222 buffered to pH 7 with sodium bicarbonate. Whole heart, liver, spleen, and tail tissues were collected from seventeen individuals that were infected and ten healthy individuals. Intact tissue was immediately preserved in RNALater, and stored at −80 °C following a 24-hour incubation at 4 °C. Tissues from six healthy and six infected specimens were included in RNA extraction, library preparation, and transcriptome assembly in order to create a more complete assembly of the A. barbouri transcriptome, including putative immune responses. Individuals were selected based on their population of origin to represent each of the three collection sites within the northern range unit. All individuals were housed in accordance with Tennessee Tech IACUC approved protocol (TTU-IACUC-21-22-008). Individuals were collected under the TWRA Scientific Collection Permits 1323 and 5669, TDEC Scientific Research and Collection Permit 2022–041, and TWRA Importation Permit 37079926.

RNA extraction, RNA-seq library construction and sequencing

Total RNA was extracted from whole larvae and heart, liver, spleen, and tail tissues using Omega Bio-tek E.Z.N.A. Tissue RNA kits (LOT numbers R693443380–51, R683448848–96, and R68345086–103) and bead homogenization using a Mini-BeadBeater-16 and 1.00 mm glass beads. All RNA extractions were treated with DNase following the DNase I digestion protocol detailed in the E.Z.N.A. Total RNA Kit II manual. Extracted RNA samples were assessed for quality and concentration using Agilent Bioanalyzer RNA 6000 Pico reagents. We performed library preparation using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina, following the manufacturer protocol and using PolyA mRNA enrichment. Approximately 200 bp RNA inserts were targeted using NEBNext Sample Purification Beads. The concentrations of the PCR products were quantified using Quantifluor dsDNA kit (Promega) and the quality was assessed using Agilent Bioanalyzer High Sensitivity DNA analysis. Equimolar amounts of individuals’ PCR products were pooled to create sequencing libraries.

Three to four individuals per temperature treatment group for each clutch were included for library preparation and RNA sequencing, totaling 80 sequenced larval individuals. Samples were selected for sequencing based on RNA integrity number (RIN) quantified via Agilent Bioanalyzer RNA 6000 Pico RNA analysis. Samples with the highest RIN per clutch-treatment group were used for sequencing; the average RIN of sequenced samples equaled 7.2. Five libraries, each consisting of 16 larval individuals, were sequenced on an Illumina NextSeq 1000 to produce 150-bp paired-end reads. For heart, spleen, liver, and tail samples, two libraries were sequenced for each tissue-treatment pairing (i.e., two spleen-infected libraries, two spleen-wild-type libraries, etc.) with each library consisting of pooled RNA from 3 individuals for each tissue. For these samples, two pools of 8 libraries each were sequenced on an Illumina NextSeq 1000 to produce 150-bp paired-end reads. Each sequencing pool contained libraries representing each of the analyzed tissues; however, infected and wild-type libraries were sequenced on separate flow cells.

RNA-seq quality assessment and control

We assessed raw-read quality across all samples with FastQC (Andrews 2010) and aggregated the results into a single report using MultiQC (Ewels et al. 2016). Quality filtering and adapter trimming was executed using FastP settings q = 30, length_required = 50, trim_front1 = 10, trim_tail1 = 10, detect_adapter_for_pe, and default settings for all other parameters (Chen et al. 2018). Random sequencing errors within reads were corrected using default settings in Rcorrector (k = 23, maxcorK = 4, wk = 0.95; Song and Florea 2015).

Transcriptome assembly and annotation

In order to create a comprehensive transcriptome for the species, we de novo assembled the trimmed and corrected reads for all of the sequenced individuals and tissues into a single transcriptome using a multi-assembly approach (Pérez-Moreno et al. 2023). In addition to three-stage de novo assembler Trinity (Grabherr et al. 2011), we employed SPAdes using the rnaSPAdes algorithm (Bushmanova et al. 2019) with the default k-mer settings as a second assembly method. The assemblies resulting from Trinity and rnaSPAdes were merged into a single non-redundant assembly using the EvidentialGene tr2aacds pipeline (Gilbert 2013). We used Kraken2 (Wood et al. 2019) to identify and remove sequences from potential sources of contamination with the applied confidence setting of 0.1. To remove redundant transcripts and fragmented assemblies we ran the hierarchical clustering tool CD-HIT-est (Li and Godzik, 2006) using default parameters and a similarity threshold of 95%, which produced unigenes. This assembly was further filtered for adaptor, vector, and foreign contamination using FCS-adaptor and FCS-GX (Tvedte 2024). Following this step, we used TransDecoder (Haas et al. 2013) with default parameter settings to predict long open reading frames (ORFs) before annotation. Annotation was performed using the software Trinotate v 4.0.1(Bryant et al. 2017) which integrates multiple tools and databases to provide a thorough annotation of the transcriptome assembly. The methods employed by Trinotate v 4.0.1 include homology search, protein domain and family identification, and protein signal peptide and transmembrane domain prediction. The databases employed by Trinotate include EggNOG-Mapper (EggNM), Gene Ontology (GO), Infernal, Kyoto Encyclopedia of Genes and Genomes (KEGG), Pfam, SignalP, and SwissProt. The assembly was annotated separately with KEGG Orthology Family (Ko Fam) information downloaded from genome.jp/keg/ko.html on February 17^th 2025.

Quality assessment

To calculate summary statistics and assess the completeness and contiguity of the transcriptome assembly, we used SeqKit2 v2.8.2 (Shen et al. 2024), BUSCO v5.0.0 (Benchmarking Universal Single-Copy Orthologs; Simão et al. 2015), and Bowtie2 v2.5.4 (Langmead and Salzberg 2012). SeqKit2 was used to calculate the number of sequences, minimum and maximum sequence length, total number of bases, GC content, N50, L50, the number of sequences under 200 base pairs in length, average sequence length, and median sequence length ( Table 1). The BUSCO databases tetrapod_obd10 and vertebrata_odb10 were queried against our assembly to obtain the number of complete, fragmented, and missing highly conserved orthologous genes found in our dataset. Lastly, concatenated, cleaned reads were mapped against the final de novo transcriptome assembly using default settings in Bowtie2 (Langmead and Salzberg 2012) to obtain the percentage of reads that mapped to the transcriptome assembly. Higher percentages of mapped reads indicate a greater quality assembly.

Dataset validation

The final assembled transcriptome had a total of 850,204 unigenes (of which the largest size was 71,527 bp) a N50 of 1,628 bp, L50 of 130,416, and E90N50 of 3,651 ( Table 1). These assembly statistics were comparable to those of other published amphibian transcriptome assemblies ( Table 2), suggesting an accurate and contiguous assembly. The high number of total unigenes reported here is not uncommon in salamander species (as demonstrated in Table 2) and could be attributed to similar isoforms within the assembly, paralogous genes, potentially fragmented unigenes, and a high number of genes expressed at low levels (Bryant et al., 2017; Grabherr et al., 2011; Smith et al., 2009). These factors may be particularly pronounced in salamanders because of their exceptionally large and repetitive genomes, which contain abundant duplicated and repetitive sequences that can increase transcript diversity and assembly complexity (Sun et al., 2012; Nowoshilow et al., 2018). It is likely that isoform diversity is a large contributor to the high total number of unigenes since reducing the CD-HIT-est similarity threshold from 99% to 95% lowered the number of unigenes from 1,046,871 to 855,398. The assembly contained complete sequences for 98.1% of the 3,354 genes corresponding to BUSCO’s database vertebrata_odb10 (Single copy: 37.6%, Duplicate: 60.5%, Fragmented: 0.6%, Missing: 1.3%) and 95.4% of the 5,310 genes corresponding to BUSCO’s database tetrapoda_odb10 (Single copy: 38.4%, Duplicate: 57.0%, Fragmented: 1.0%, Missing: 3.6%; Figure 1). Similarly to the high number of unigenes, the high duplicate BUSCO percentage observed may be attributed to similar isoforms of BUSCO markers and/or high allelic variation which could be resolved by applying a lower similarity threshold in CD-HIT-est. However, further collapsing similar isoforms and allelic variation may remove genetic diversity that is essential for downstream differential expression analyses. Read mapping rates to the assembled transcriptome were high (Bowtie2 overall alignment rate 97.34%), indicating strong representation of the sequenced reads in the assembly; the majority of mapped reads were multi-mapped (80.99%), consistent with isoform redundancy and shared sequence among alternatively spliced transcripts. Of 534,476 TransDecoder identified protein coding sequences, 276,770 were annotated with information from at least one of the applied annotation databases (EggNM, GO, Infernal, KEGG, Ko Fam, Pfam, SignalP, SwissProt; Figure 2).

Figure 1. Ambystoma barbouri transcriptome assembly BUSCO results, calculated with both vertebrate_odb10 and tetrapoda_odb10 lineages show a high percentage of complete BUSCO markers found within the assembly.

Complete BUSCOs found in the assembly are represented by different shades of blue (light blue for complete, single-copy; dark blue for complete, duplicate copy). Markers missing from the assembly are represented by red and fragmented BUSCOs are represented by yellow.

Figure 2. Ambystoma barbouri transcriptome annotation breakdown of the number of unigenes annotated across eight annotation databases: EggNOG-mapper (EggNM), gene ontology (GO), infernal, kyoto encyclopedia of genes and genomes (KEGG), KEGG orthology family (Ko Fam), Pfam, SignalP, and SwissProt.