Comparative transcriptomes of three different skin sites for the Asiatic toad (Bufo gargarizans)

Abstract

Toads release toxic dry secretions from glands in their skin. Toxin possesses a wide range of biological effects, but little is known about its specific gene expression pattern and regulatory mechanisms. The Asiatic toad (Bufo gargarizans) is widely used to produce toxin. Here, we explored the gene expression of 30 tissue samples from three different skin sites (parotoid gland, dorsal skin, and abdomen skin) of B. gargarizans. After de novo assembly, 783,130 unigenes with an average length of 489 bp (N50 = 556 bp) were obtained. A total of 9,248 significant differentially expressed genes (DEGs) were detected. There were 8,819 DEGs between the parotoid gland and abdomen skin and 1,299 DEGs between the dorsal skin and abdomen skin, while only 1,283 DEGs were obtained between the parotoid gland and dorsal skin. Through enrichment analysis, it was found that the detected differential gene expressions corresponded to the different functions of different skin sites. Our key findings were the genetic expression of toxin secretion, the protection function of skin, and the related genes such as HSD3B, Cyp2c, and CAT, LGALS9. In conclusion, we provide useful transcript resources to study the gene expression and gene function of B. gargarizans and other amphibians. The detected DEGs between different sites of the skin provided better insights into the genetic mechanisms of toxin secretion and the protection function of skin for amphibians.

Introduction

The Asiatic toad (Bufo gargarizans) belongs to the family of Bufonidae (Anura: Amphibia) and is widely distributed in East Asia. As a true toad, B. gargarizans has a pair of well-developed parotoid glands and numerous granular glands on the back (Lu et al., 2021). These glands secrete potent alkaloid toxins (collectively bufotoxins), which consist of structurally and functionally similar toxic compounds (e.g., bufalin) (Hayes et al., 2009), to defend themselves against predators, parasites, and pathogens (Mina et al., 2015). In addition, the toxin derived from toads and the toxin possesses a wide range of biological effects, having long been used for the treatment of heart failure, tumors, sores, and pains in clinical settings in traditional Chinese medicine (Wu, Gao & Wang, 2010; Ji et al., 2014; State Pharmacopoeia Commission of the PRC, 2020).

In this study, we explored the question of what physiological differences exist in different skin sites in order to deal with the toxin presence. Because Amphibians played a very important role in the process of animal evolution and lived in a special ecological environment, they formed their own defense system (Zhou et al., 2021). The B. gargarizans could produce potent toxins to serve as weapons to fight their predators and to protect themselves (Ujvari et al., 2015). Physiological and biochemical analyses have identified the toxin contains chemical compounds classified into four categories: biogenic amines, bufadienolides, alkaloids, and peptides and proteins (Abdel-Rahman, Ahmed & Nabil, 2010; Mariano et al., 2018; Mariano et al., 2019). Some compounds belonged to the bufadienolides were considered the main active constituents of the toxin, such as bufalin, cinobufagin, arenobufagin, and resibufogenin (Li et al., 2021). Indole alkaloids were considered another type of active constituents of the toxin. Compared with the bufadienolides, there were relatively few reports on antitumor activity of Indole alkaloids (Li et al., 2021). One of the major active constituents of the toxin is bufalin (Wang et al., 2015; Liu et al., 2016). Although previous analyses have identified main compounds (Wu et al., 2020), the specific gene expression pattern and regulatory mechanisms of toxin and its components remain unclear.

RNA-seq is a powerful approach to analysis transcriptomes, which could generate global view of the transcriptional landscapes and correlate protein abundance with mRNA content (Wang, Gerstein & Snyder, 2009; Zhou et al., 2020; Wang et al., 2021). The transcriptomes of some tissues, such as liver, brain, heart, muscle and testicle from B. gargarizans, were sequenced (Yang, Qi & Fu, 2016; Yang et al., 2017; Jin et al., 2018), the gene expression profiles of skin were still missing. Yang et al. found the gene expression variations and genetic signals of high-altitude adaptation in B. gargarizans (Yang, Qi & Fu, 2016; Yang et al., 2017; Jin et al., 2018). Jin et al. (2018) sequenced liver transcriptomes of active and torpid B.gargarizans and gained insights into changes in gene expression patterns in hibernating female and male toads. Besides, although Lu et al. (2021) presented a high-quality genome assembly for B. gargarizans, the gene annotation was unavailable. Therefore, we applied transcriptomic methods to analyze the differentially expressed genes (DEGs) in different skin sites of B. gargarizans to explore the differential gene expressions. Transcriptome sequencing was performed on 30 samples from the parotoid gland, dorsal skin, and abdomen skin of ten individuals using Illumina paired-end sequencing technology. Our main purposes were to: (1) complete the assembly of the transcriptome of B. gargarizans; (2) annotate the skin transcriptome of B. gargarizans; (3) calculate and analyze gene expression; (4) screen out the DEGs among the three skin sites; (5) perform functional enrichment of DEGs to further clarify the gene expression profiles at different skin sites of B. gargarizans and (6) investigate the genetic mechanisms may related to self-protection, toxin secretion and production of anti-tumor components.

Materials and Methods Experimental design and sample collection

Ten wild adult female Asiatic toads (B. gargarizans) were obtained from the wetlands by the lake in Chengdu, Sichuan Province, China in Mid-October 2018. All toads were captured on the same day and lived in similar environments. All animal experiments in this study were approved by the Ethics Committee of the College of Life Sciences, Sichuan University (number 20210309009). The ten toads were either anesthetized and euthanized. After euthanasia, 30 tissue samples were taken from three skin sites; parotoid gland, dorsal skin and abdomen skin. The 30 samples from the three skin sites were collected in cryotubes and immediately stored in liquid nitrogen for later RNA extraction and qPCR.

Library preparation and RNA sequencing

Total RNA was isolated separately from the 30 samples using a TRIzol Kit (Invitrogen, Waltham, MA, USA) according to the manufacturer’s instructions. The construction of a cDNA library and transcriptome sequencing was conducted in Novogene (Beijing, China), using the PolyATract mRNA Isolation Systems Kit (Promega) and Illumina HiSeq 4000 according to the manufacturer’s instructions.

Sequence preprocessing and transcripts assembly

We obtained the 150-bp paired-end (PE) short reads after Illumina sequencing was completed. The adapters and low-quality reads in raw reads were removed by NGS QC Toolkit v2.3.3 (Patel, Mukesh & Liu, 2012) with stringent criteria (reads with more than 90% bases within ≥ 20 base quality (Q-value)) to obtain clean reads. The transcriptome was de novo assembled using the short-read assembly program Trinity (version 2.8.6) software with the option “-min_contig_length 300” (Grabherr et al., 2011; Haas et al., 2013). After assembly, we used TrinityStats.pl, which comes with Trinity, to counted the number, size and N50 of components and transcripts. The assembled transcripts were clustered by CD-HIT (version 4.8.1) (Fu et al., 2012) with the threshold of 95% similarity and 90% coverage of query sequences to obtain the non-redundant transcript sets as the unigenes. To evaluate the completeness, we employed Benchmarking Universal Single-Copy Orthologs (BUSCO v5.1.2; http://busco.ezlab.org/) to evaluate the gene set of the toad using Vertebrata data (Seppey, Manni & Zdobnov, 2019). The result of BUSCO was shown in Fig. 1A. The protein-coding regions prediction of assembled unigenes was undertaken using TransDecoder-v5.5.0 and the longest open reading frame (ORF) predicted for each contig sequence with a minimum of 100 amino acids long (Haaset al.2013). The result of ORF prediction was shown in Fig. 1B.

Figure 1: Sequencing and de novo assembly. (A) The results of BUSCO estimation. (B) Venn-diagram showing the classification of all identified open reading frames (ORFs). (C) Venn-diagram showing the overlap between results of unigenes aligned to Swiss-Prot, KEGG, COGs, and GO databases. (D) The principal component analysis (PCA) of 29 samples. This analysis used gene expression values, which normalized by Deseq2, and the differences between individuals are eliminated. DOI: 10.7717/peerj.12993/fig-1

Gene annotation

Unigenes from the de novo assembly were aligned into the UniProt database (Rolf et al., 2004) and Clusters of Orthologous Groups (COGs) (Koonin et al., 2004) using DIAMOND (version 2.0.7.145) (Huson & Buchfink, 2015). The unigenes were also aligned into the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Gerlich & Neumann, 2000) database using KAAS (Yuki et al., 2007) using the single-directional best-hit (SBH) method to predict hypothetic function and the background gene set we selected were Homo sapiens, Xenopus laevis, Xenopus tropicalis and Nanorana parkeri. Unigenes aligned to the UniProt database were assigned to Gene Ontology (GO) (Ashburner et al., 2000) terms according to the idmapping_selected.tab file (last modified on 2021-02-10). The WEGO chart was generated by uploading the GO annotation results to http://wego.genomics.org.cn/ (Ye et al., 2018). The results of gene annotation were shown in Fig. 1C.

Differentially expressed gene and gene enrichment analysis

The quantification of these genes in each sample was performed using RNA-Seq by Expectation Maximization (RSEM) (Dewey & Li, 2011) after mapping to the assembled unigenes using Bowtie2 (version 2.4.2) (Langdon, 2015). We also mapped the samples to the reference genome of Asiatic toads from NCBI (GCA_014858855.1; Lu et al., 2021) using Bowtie2 (version 2.4.2) (Table S8). We eliminated the interference of low expression genes by removing genes from the gene count matrix obtained in the previous step that met all the following conditions: the number of samples with non-zero expression was less than ten. We then used R (version 3.6.2) package DESeq2 (Love, Huber & Anders, 2014) to identify the DEGs, with processed raw read counts as the input. We divided all samples into three groups according to the sampling site: dorsal skin, abdomen skin and parotoid gland. The read counts matrix and the samples’ information were constructed to a DESeqDataSet (dds) by the function, DESeqDataSetFromMatrix. After constructing dds, we performed a minimal pre-filtering to keep only rows that have at least two reads total. The standard differential expression analysis steps were wrapped into a single function, DESeq. Based on the normalized abundances of genes, principal component analysis (PCA) was performed using the function, plotPCA of DESeq2 (version 1.26.0) to visualize the correlation of the 30 RNA-seq samples (Fig. 1D). The Benjamini–Hochberg false discovery rate (FDR) was used to correct P values for multiple testing (FDR ≤ 0.05) and an absolute value of log 2 fold change (log2fc) was used to determine the significant differences in gene expression (log2fc ≥ 2). In this study, we focused on the DEGs that annotated using Swiss-Prot, GO and KEGG. GO and KEGG enrichment analyses were performed by R package clusterProfiler (version 3.14.3) (Yu et al., 2012) with GO and KEGG annotation results of all unigenes of the toad as the background gene set. After we completed the GO and KEGG enrichment of the differential genes, we selected the genes contained in the GO terms and KEGG pathways involved in the toxin secretion and self-protection of the toad. After we got the gene ids related to the toxin secretion and self-protection of the toad, we would get the list of UniProt Accession of these genes according to the Swiss-Prot annotation results (the dna_matches_fmt6.csv in the supplementary files). Finally, we would enter the list of UniProt Accession of these genes in bioDBnet (https://biodbnet-abcc.ncifcrf.gov/db/db2db.php) to get their gene symbols. The 18 genes we selected were shown in Fig. 2.

Figure 2: The results of differences in expression of 18 DEGs in three different skin parts. We calculate the 18 DEGs’ (ABCD3, Acox2, CROT, Pex10, Pecr, Slc27a2, XDH, EHHADH, TPH1, KYAT3, MAOB, Cyp2c, Cyp2d, GPX4, HSD3B, CAT, LGALS9, LTR1) expression level by the gene count matrix. DOI: 10.7717/peerj.12993/fig-2

Results Sampling, sequencing, and de novo assembly

We collected 30 tissue samples of three skin sites from ten wild adult female Asiatic toads (B. gargarizans) in Mid-October 2018. The three skin sites were the parotoid gland, dorsal skin, and abdomen skin. We then constructed 30 cDNA libraries and one of the dorsal skin samples was removed due to the unqualified cDNA library. Finally, a total of 813,227,650 raw reads were generated and 799,684,651 clean reads were obtained from 29 samples. Considering that the genome assembly of the Asiatic toads (B.gargarizans) lacked the corresoonding genome annotation file (.gtf file), we assembled the transcripts by de novo assembly. The clean reads were further assembled into 783,130 unigenes with an average length of 489.2 bp (N50 = 556 bp) after clustering and removing redundancy (Table 1). The average mapping rates of the samples that mapped to the reference genome and assembled unigenes were 62.41% and 88.58% respectively. We found that 45.1% of total complete and single-copy BUSCOs and 46.7% complete and duplicated BUSCOs were identified in this gene set, while 8.2% fragmented BUSCOs and nearly zero missing BUSCOs were identified in this gene set (Fig. 1A). A total of 117,553 ORFs were identified from 103,749 of the 783,130 assembled unigenes by TransDecoder, and 93,919 of the unigenes contained only one ORF (Fig. 1B).

Table 1: Summary of transcripts assembly.

total transcripts	783,130
max length	44,570
min length	201
total length	383106510.0
mean length	489.2
median length	327.0
Contig N10	2,627
Contig N50	556
Contig N90	241
GC content	46.22%

DOI: 10.7717/peerj.12993/table-1

Functional annotation of unigenes

We further annotated the function of the 783,130 assembled unigenes based on sequence similarity. Among these, 111,065 (14.00%), 35,526 (4.54%), 63,490 (8.11%) and 108,439 (13.80%) unigenes were aligned to Swiss-Prot, KEGG, COGs, and GO databases, respectively, and 11,931 unigenes were annotated across all databases (Fig. 1C). The GO annotation showed that the major subcategories were “Cell (GO:0005623)”, “Cell part (GO:0044464)”, “binding (GO:0005488)”, “catalytic activity (GO:0003824)”, “cellular process (GO:0009987)” and “metabolic process (GO:0008152)” (Fig. S1). COG assignments were performed to predict and classify the possible gene functions (Fig. S2). The hits from the COG prediction were classified into 26 functional categories, where the most enriched terms were translation, ribosomal structure and biogenesis (6,075 unigenes, J in Fig. S2), followed by amino acid transport and metabolism (5,860 unigenes, E in Fig. S2) and general function prediction only (5,116 unigenes, R in Fig. S2).

Identification of DEGs and functional enrichments

We performed a PCA analysis on gene expression from the 29 samples with the normalized RNA-seq data and found that the three parts clearly clustered in the first Principal Component (PC1) and PC2, which explained 36% of the variance (Fig. 1D). After that, DEGs were analyzed with pairwise comparison of the three tissue types. A total of 9,248 significant DEGs were detected by DESeq2 (Table S1). DEGs with higher expression levels in one group, when compared to another, were denoted as “upregulated”, while those with lower expression levels were termed as “downregulated”. In general, 8,819 DEGs were identified between parotoid gland and abdomen skin, while only 1,283 DEGs were obtained between parotoid gland and dorsal skin. There were 1299 DEGs between dorsal skin and abdomen skin (Table S1).

We performed GO category and KEGG pathway enrichment analyses to gain insight into the biological roles of DEGs. First, GO and KEGG enrichment analyses were performed with 8819 DEGs between the parotoid gland and abdomen skin. Arachidonic acid metabolic process (GO:0019369), UV protection (GO:0009650), and catalase activity (GO:0004096) were mainly enriched by the upregulated DEGs in the parotoid gland (Table S2), whereas the downregulated DEGs in the parotoid gland were significantly enriched in extracellular matrix organization (GO:0030198), regulation of inflammatory response (GO:0050727), and calcium-dependent cysteine-type endopeptidase activity (GO:0004198) (Table S3). For KEGG enrichment, the upregulated DEGs in the parotoid gland were significantly enriched in peroxisome (ko04146), tryptophan metabolism (ko00380), and FoxO signaling pathway (ko04068) (Table S2), whereas the downregulated DEGs in the parotoid gland were significantly enriched in Calcium signaling pathway (ko04020) (Table S3).

Next, we performed GO and KEGG enrichment with the 1,299 DEGs found between dorsal skin and abdomen skin. The upregulated DEGs of dorsal skin were mainly enriched in steroid biosynthetic process (GO:0006694), arachidonic acid 14,15-epoxygenase activity (GO:0008404), and defense response to gram-negative bacterium (GO:0050829) (Table S4), while the downregulated DEGs of dorsal skin were enriched in melanin biosynthetic process (GO:0042438), embryonic forelimb morphogenesis (GO:0035515) and N-acetylmuramoyl-L-alanine amidase activity (GO:0008745) (Table S5). Glutathione metabolism (ko00480), Glutathione metabolism (ko00480), and Amphetamine addiction (ko00140) were significantly enriched by the upregulated DEGs of dorsal skin in the KEGG enrichment analysis (Table S4), while Neuroactive ligand–receptor interaction (ko04080), and Adrenergic signaling in cardiomyocytes (ko04261) were significantly enriched by the downregulated DEGs of dorsal skin in the KEGG enrichment analysis (Table S5).

Finally, we investigated the 1,283 DEGs identified between the parotoid gland and dorsal skin. The upregulated DEGs in the parotoid gland included enrichment in p38MAPK cascade (GO:0038066), positive regulation of CD4-positive, CD25-positive, alpha-beta regulatory T cell differentiation involved in immune response (GO:0032834), kynurenine metabolic process (GO:0070189), and cysteine-S-conjugate beta-lyase activity (GO:0047804) (Table S6). Whereas the downregulated DEGs in the parotoid gland were mainly enriched in xenobiotic metabolic process (GO:0006805) and calcium-dependent cysteine-type endopeptidase activity (GO:0004198) (Table S7). Glycosaminoglycan biosynthesis - chondroitin sulfate/dermatan sulfate (ko00532), and Amphetamine addiction (ko05031) were significantly enriched by the upregulated DEGs of the parotoid gland in KEGG enrichment (Table S6), while Nitrogen metabolism (ko00910) and Necroptosis (ko04217) were significantly enriched by the downregulated DEGs of the parotoid gland in KEGG enrichment (Table S7).

Combining our enrichment results and annotation results, we discovered 18 DEGs (ABCD3 (p-values: dorsal skin vs parotoid gland: 0.00041; dorsal skin vs abdomen skin: 0.000087; parotoid gland vs abdomen skin: 0.000011), Acox2 (p-values: dorsal skin vs parotoid gland: 0.043; dorsal skin vs abdomen skin: 0.017; parotoid gland vs abdomen skin: 0.00021), CAT (p-values: dorsal skin vs parotoid gland: 0.00041; dorsal skin vs abdomen skin: 0.00041; parotoid gland vs abdomen skin: 0.000011), CROT (p-values: dorsal skin vs parotoid gland: 0.00065; dorsal skin vs abdomen skin: 0.095; parotoid gland vs abdomen skin: 0.000011), Cyp2c (p-values: dorsal skin vs parotoid gland: 0.0057; dorsal skin vs abdomen skin: 0.000087; parotoid gland vs abdomen skin: 0.000011), Cyp2d (p-values: dorsal skin vs parotoid gland: 0.16; dorsal skin vs abdomen skin: 0.013; parotoid gland vs abdomen skin: 0.001), EHHADH (p-values: dorsal skin vs parotoid gland: 0.000043; dorsal skin vs abdomen skin: 0.000043; parotoid gland vs abdomen skin: 0.000011), GPX4 (p-values: dorsal skin vs parotoid gland: 0.0015; dorsal skin vs abdomen skin: 0.00026; parotoid gland vs abdomen skin: 0.000011), HSD3B (p-values: dorsal skin vs parotoid gland: 0.72; dorsal skin vs abdomen skin: 0.00097; parotoid gland vs abdomen skin: 0.00013), KYAT3 (p-values: dorsal skin vs parotoid gland: 0.00065; dorsal skin vs abdomen skin: 0.00041; parotoid gland vs abdomen skin: 0.000011), LGALS9 (p-values: dorsal skin vs parotoid gland: 0.000087; dorsal skin vs abdomen skin: 0.00028; parotoid gland vs abdomen skin: 0.00018), LTR1 (p-values: dorsal skin vs parotoid gland: 0.0062; dorsal skin vs abdomen skin: 0.14; parotoid gland vs abdomen skin: 0.0015), MAOB (p-values: dorsal skin vs parotoid gland: 0.000087; dorsal skin vs abdomen skin: 0.0006; parotoid gland vs abdomen skin: 0.000011), Pecr (p-values: dorsal skin vs parotoid gland: 0.00065; dorsal skin vs abdomen skin: 0.0055; parotoid gland vs abdomen skin: 0.00018), Pex10 (p-values: dorsal skin vs parotoid gland: 0.003; dorsal skin vs abdomen skin: 0.19; parotoid gland vs abdomen skin: 0.000043), Slc27a2 (p-values: dorsal skin vs parotoid gland: 0.013; dorsal skin vs abdomen skin: 0.013; parotoid gland vs abdomen skin: 0.000022), TPH1 (p-values: dorsal skin vs parotoid gland: 0.095; dorsal skin vs abdomen skin: 0.028; parotoid gland vs abdomen skin: 0.0021), XDH (p-values: dorsal skin vs parotoid gland: 0.043; dorsal skin vs abdomen skin: 0.028; parotoid gland vs abdomen skin: 0.000011)) that might relevant to toads’ self-protection and toxin secretion. We also calculate their expression levels in dorsal skin, parotoid glands and abdmen skin (Fig. 2).

Discussion

To date, there have been only a few published complete genome sequences of amphibians mainly due to their large and complex genomes (Lu et al., 2021). Therefore, RNA-seq technology is a potentially useful tool to study their gene annotation, gene discovery, and gene expression (Geng et al., 2019). Recently, the complete genome sequence of the Asiatic toad was released (Lu et al., 2021), which provided a solid foundation to understand its genetic background, but the Asiatic toad genome still lacked the corresponding annotation files. In addition, there were some transcript resources of the Asiatic toads’ liver, brain, heart, muscle, testicles and oothecae (Yang, Qi & Fu, 2016; Yang et al., 2017; Jin et al., 2018), but the skin transcripts, an important tissue of the Asiatic toads, were still missing. In this study, we sequenced and analyzed the 29 skin transcriptomes from three sites on ten Asiatic toads. We constructed high-quality reference gene sets for the skin tissue and further analyzed the specific expression between different sites of the skin.

The toxin from dried secretions of the parotoid glands and skin glands of the Asiatic toad has been used in Traditional Chinese Medicine (TCM) for treating infection and inflammation for hundreds of years (Qi et al., 2014; Zhan et al., 2020). It has also been widely used as animal medicine in other Asian countries, such as Japan and Korea (The Korea Food and Drug Administration, 2016; Editorial board of Japan’s Ministry of Medicine, 2017). Both the parotoid gland and dorsal skin secrete toxin, but the white secretion from the parotoid gland produces the largest amount. Abdomen skin lacks the gland and thus cannot secrete toxin (Xia et al., 2019). We identified the largest number of DEGs (8,819) between the parotoid gland and abdomen skin, whereas there was only 1,283 DEGs between the parotoid gland and dorsal skin. The 1,283 DEGs between the parotoid gland and dorsal skin was also smaller than the 1,299 DEGs between dorsal skin and abdomen skin. The detected differential gene expressions were corresponding to the functions of different sites of the skin tissue. Therefore, the difference of gene expression profiles between the parotoid gland and abdomen skin was considerable, and the 8,819 DEGs may relate to the secretion of toxin.

We identified upregulated DEGs in the parotoid gland that were enriched in several important GO terms associated with toxin synthesis, such as Arachidonic acid metabolic process (GO:0019369). The DEGs enriched in Arachidonic acid metabolic process (GO:0019369) included GPX4, Cyp2c gene family (CYP2C28, Cyp2c29, Cyp2c50, Cyp2c23, Cyp2c18, Cyp2c12, Cyp2c7, Cyp2c19, Cyp2c70 and Cyp2c8) (cytochrome P450 family 2 subfamily C) and Cyp2d gene family (Cyp2d1, CYP2D6, Cyp2d26 and CYP2D15) (Fig. 2) (cytochrome P450 family 2 subfamily D). The protein encoded by GPX4 belongs to the glutathione peroxidase family, members of which catalyze the reduction of hydrogen peroxide, organic hydroperoxides and lipid hydroperoxides, and thereby protect cells against oxidative damage (Huang & Wu, 2020). According to previous research, arachidonic acid was identified as a precursor of toxin (Bruzzone, Buschiazzo & Alonso, 2003). Bufalin, the primary active ingredient of toxin, is a cardiotonic steroid. Therefore, Arachidonic acid metabolic process (GO:0019369), steroid hydroxylase activity (GO:0008395), steroid delta-isomerase activity (GO:0004769) and epoxygenase P450 pathway (GO:0019373) may be associated with toxin secretion. Some of the DEGs are particularly important and need further study. For example, the DEG HSD3B, enriched in steroid biosynthetic process (GO:0006694), encoded an enzyme that catalyzes the oxidative conversion of delta-5-3-beta-hydroxysteroid precursors into delta-4-ketosteroids, which leads to the production of all classes of steroid hormones (Zein et al., 2020). The encoded protein also catalyzes the interconversion of 3-beta-hydroxy- and 3-keto-5-alpha-androstane steroids (Hearn et al., 2020). Furthermore, the Cyp2c and Cyp2d gene family that encode cytochrome P450 family 2 is involved in regulating the biosynthesis and metabolism of steroids (Zhang et al., 2019). The cytochrome P450 proteins are monooxygenases, which catalyze many reactions involved in the synthesis of cholesterol, steroids, and other lipids (Gómez-Tabales et al., 2020). The steroid is a water-soluble component in toxin (Zhanget al.2019). Steroids and bufalin isolated from toxin can induce apoptosis of human tumor cells by regulating the expression of the caspase gene family in cells (Gao et al., 2017).

In addition, several amino acid related terms or pathways were enriched with many upregulated parotoid gland DEGs, such as tryptophan metabolism (ko00380), which contained five DEGs (CAT, MAOB, EHHADH, TPH1 and KYAT3) (Fig. 2). The TPH1 encodes a member of the biopterin-dependent aromatic amino acid hydroxylase family (Choi et al., 2018). The encoded protein is one of two tryptophan hydroxylase enzymes that catalyze the first and rate the limiting step in the biosynthesis of the hormone and neurotransmitter, serotonin (Li & He, 2006). The encoded protein is involved in the development of hypoxia-induced elevations in pulmonary pressures and pulmonary vascular remodeling, and has also been implicated as a regulator of immune tolerance (Flamar et al., 2020). The KYAT3 encodes an aminotransferase that transaminates kynurenine to form kynurenic acid, which is a metabolite of tryptophan (Matysik-Woźniak et al., 2017). The peroxisome (ko04146) was also enriched with many DEGs, such as CAT, Slc27a2, EHHADH, Pex10, Pecr, CROT, Acox2, XDH and ABCD3 (Fig. 2). The protein encoded by Slc27a2 is an isozyme of long-chain fatty-acid-coenzyme A ligase family (Tharp et al., 2016). Although differing in substrate specificity, subcellular localization, and tissue distribution, all isozymes of this family convert free long-chain fatty acids into fatty acyl-CoA esters, and thereby play a key role in lipid biosynthesis and fatty acid degradation (Veglia et al., 2019). The CROT encodes a member of the carnitine/choline acetyltransferase family. The encoded protein converts 4,8-dimethylnonanoyl-CoA to its corresponding carnitine ester. This transesterification occurs in the peroxisome and is necessary for transport of medium- and long- chain acyl- CoA molecules out of the peroxisome to the cytosol and mitochondria. The protein thus plays a role in lipid metabolism and fatty acid beta-oxidation (Okui et al., 2021). The product of Acox2 belongs to the acyl-CoA oxidase family. It encodes the branched-chain acyl-CoA oxidase that is involved in the degradation of long branched fatty acids and bile acid intermediates in peroxisomes (Valentino et al., 2017).

According to previous studies, bufalin can enhance the level of interleukin in mouse lymphocytes with T lymphocyte-dependent immune response (Cho et al., 2018). Modern pharmacological studies show that bufalin has antitumor effects in various cancers, such as leukemia as well as lung, gastric, breast, prostate, ovarian and bladder cancers (Liuet al.2016). LGALS9 is involved in the positive regulation of interleukin synthesis, so it can be inferred that LGALS9 is a key gene involved in regulating the secretion of toxin (Qi et al., 2011).

In conjunction with the genetic expression of toxin secretion, we also explored the protective function from other environmental stressors of the skin. The skin is an important tissue to protect the Asiatic toad, which is the first line of defense against bacteria and other external stimuli (Gabriel & Costa, 2021). UV protection (GO:0009650) was enriched by the upregulated genes in the parotoid gland. The expression levels of CAT in the posterior ear glands and back skin were higher than in the abdominal skin (Fig. 2). The CAT gene is involved in the protection of the Chinese toad from ultraviolet rays (Muramatsu et al., 2005). This gene encodes catalase, a key antioxidant enzyme in the body’s defense against oxidative stress (Yenkoyan, Harutyunyan & Harutyunyan, 2018). The high expression level of these genes could play positive roles in the protection of the Chinese toad from ultraviolet rays. In addition, other genes that might relate to the survival of the Chinese toad were also observed as DEGs. For example, defense response to Gram-negative bacterium (GO:0050829), defense response to Gram-positive bacterium (GO:0050830), and cellular response to virus (GO:0098586) were enriched by the upregulated DEGs in dorsal skin and downregulated in abdomen skin. Positive regulation of NIK/NF-kappaB signaling (GO:1901224) was enriched by the upregulated DEGs in the dorsal skin. A family of Toll-like receptors (TLRs) acts as primary sensors that detect a wide variety of microbial components and elicit innate immune responses. All TLR signaling pathways culminate in activation of the transcription factor nuclear factor-kappaB (NF-κB), which controls the expression of an array of inflammatory cytokine genes (Kawai & Akira, 2007). The protein encoded by LTR1 is a member of the Toll-like receptor (TLR) family which plays a fundamental role in pathogen recognition and activation of innate immunity. LTR1 is ubiquitously expressed, and at higher levels than other TLR genes. They recognize pathogen-associated molecular patterns (PAMPs) that are expressed on infectious agents, and mediate the production of cytokines necessary for the development of effective immunity (Lanki et al., 2019).

In conclusion, we assembled and annotated a comprehensive gene set of the three different skin sites for B. gargarizans. This gene set could provide useful transcript resources to study the gene expression and gene function of B. gargarizans and other amphibians. The detected DEGs between different sites of the skin could provide better insights into the genetic mechanisms of toxin secretion and the protection function of the skin.

Additional Information and Declarations

Competing Interests

The authors declare there are no competing interests.

Author Contributions

Yue Lan and Lewei He analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft.

Xue Dong and Ruixiang Tang performed the experiments, analyzed the data, prepared figures and/or tables, and approved the final draft.

Wanyu Li performed the experiments, authored or reviewed drafts of the paper, and approved the final draft.

Jiao Wang performed the experiments, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft.

Lei Wang performed the experiments, prepared figures and/or tables, and approved the final draft.

Bisong Yue, Megan Price, Tao Guo and Zhenxin Fan conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft.

Animal Ethics

The following information was supplied relating to ethical approvals (i.e., approving body and any reference numbers):

the Ethics Committee of the College of Life Sciences, Sichuan University

Data Availability

The following information was supplied regarding data availability:

The raw RNA-seq data are available at the China National GeneBank DataBase (CDGBdb): CNP0001870.

Funding

This work was supported by the Second Tibetan Plateau Scientific Expedition and Research Program (STEP), Grant No. 2019QZKK05010208. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

Abdel-Rahman MA, Ahmed SH, Nabil ZI. 2010. In vitro cardiotoxicity and mechanism of action of the Egyptian green toad Bufo viridis skin secretions. Toxicology in Vitro 24:480-485

Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT. 2000. Gene ontology: tool for the unification of biology. Nature Genetics 25(1):25-29

Bruzzone A, Buschiazzo J, Alonso TS. 2003. Lipids during Bufo arenarum oogenesis. Zygote 11:95-100

Cho YJ, Lee JE, Park MJ, O’Malley BW, Han SJ. 2018. Bufalin suppresses endometriosis progression by inducing pyroptosis and apoptosis. Journal of Endocrinology 237:255-269

Choi W, Namkung J, Hwang I, Kim H, Lim A, Park HJ, Lee HW, Han K-H, Park S, Jeong J-S, Karsenty G, Ju YS, Choi C, Suh JM, Park JY, Park S, Kim H+7 more. 2018. Serotonin signals through a gut-liver axis to regulate hepatic steatosis. Hepatology 68:169A-70A

Dewey CN, Li B. 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12:323-

Flamar AL, Klose CSN, Moeller JB, Mahlakõiv T, Bessman NJ, Zhang W, Moriyama S, Stokic-Trtica V, Rankin LC, Putzel GG, Rodewald HR, He Z, Chen L, Lira SA, Karsenty G, Artis D+6 more. 2020. Interleukin-33 induces the enzyme tryptophan hydroxylase 1 to promote inflammatory group 2 innate lymphoid cell-mediated immunity. Immunity 52:606-19.e6

Fu L, Niu B, Zhu Z, Wu S, Li W. 2012. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28(23):3150-3152

Gabriel A, Costa S. 2021. Effects of long-term exposure to increased salinity on the amphibian skin bacterium erwinia toletana. 80:779-788

Gao F, Wang X, Li Z, Zhou A, Tiffany-Castiglioni E, Xie L, Qian Y. 2017. Identification of anti-tumor components from toad venom. Oncology Letters 14:15-22

Geng X, Zhang L, Zang X, Guo J, Xu C. 2019. RNA-seq analysis provides insight into molecular adaptations of Andrias davidianus. 229:197-206

Gerlich M, Neumann S. 2000. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Research 28:27-30

Gómez-Tabales J, García-Martin E, Agúndez JAG, Gutierrez-Merino C. 2020. Modulation of CYP2C9 activity and hydrogen peroxide production by cytochrome b(5) Scientific Reports 10(1):1557

Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q. 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology 29:644-652

Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber MD, Ott M, Orvis J, Pochet N, Strozzi F, Weeks N, Henschel R, LeDuc RD, Friedman N, Regev A+9 more. 2013. Transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nature Protocols 8(8):1494-1512

Hayes RA, Crossl MR, Hagman M, Capon RJ, Shine R. 2009. Ontogenetic Variation in the Chemical Defenses of Cane Toads (Bufo marinus): toxin Profiles and Effects on Predators. Journal of Chemical Ecology 35:391-399

Hearn JWD, Sweeney CJ, Almassi N, Reichard CA, Reddy CA, Li H, Hobbs B, Jarrard DF, Chen YH, Dreicer R, Garcia JA, Carducci MA, Di Paola RS, Sharifi N+4 more. 2020. HSD3B1 genotype and clinical outcomes in metastatic castration-sensitive prostate cancer. JAMA Oncology 6:e196496

Huang YY, Wu CH. 2020. Association between the genetic variants of glutathione peroxidase 4 and severity of endometriosis. 17(14):5089

Huson DH, Buchfink B. 2015. Fast and sensitive protein alignment using DIAMOND. Nature Methods 12(1):59-60

Ji Q, Tan CK, Hashimi SM, Zulfiker A, Ming QW. 2014. Toad glandular secretions and skin extractions as anti-inflammatory and anticancer agents. Evidence-Based Complementary and Alternative Medicine 2014:312684

Jin L, Yu JP, Yang ZJ, Merila J, Liao WB. 2018. Modulation of gene expression in liver of hibernating asiatic toads (Bufo gargarizans) International Journal of Molecular Sciences 19(8):2363

Kawai T, Akira S. 2007. Signaling to NF-kappaB by Toll-like receptors. Trends in Molecular Medicine 13:460-469

Koonin EV, Fedorova ND, Jackson JD, Jacobs A. 2004. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biology 5:R7-R

Langdon WB. 2015. Performance of genetic programming optimised Bowtie2 on genome comparison and analytic testing (GCAT) benchmarks. BioData Mining 8(1):1

Lanki M, Seppänen H, Mustonen H, Hagström J, Haglund C. 2019. Toll-like receptor 1 predicts favorable prognosis in pancreatic cancer. PLOS ONE 14:e0219245

Li DW, He L. 2006. Further clarification of the contribution of the tryptophan hydroxylase (TPH) gene to suicidal behavior using systematic allelic and genotypic meta-analyses. Human Genetics 119:233-240

Li FJ, Hu JH, Ren X, Zhou CM, Liu Q, Zhang YQ. 2021. Toad venom: a comprehensive review of chemical constituents, anticancer activities, and mechanisms. Archiv der Pharmazie (Weinheim) 354:e2100060

Liu X, Xiao XY, Shou QY, Yan JF, Chen L, Fu HY, Wang JC. 2016. Bufalin inhibits pancreatic cancer by inducing cell cycle arrest via the c-Myc/NF-κB pathway. Journal of Ethnopharmacology 538-545

Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550

Lu B, Jiang J, Wu H, Chen X, Song X, Liao W, Fu J. 2021. A large genome with chromosome-scale assembly sheds light on the evolutionary success of a true toad (Bufo gargarizans) Molecular Ecology Resources Epub ahead of print 2021 11 January

Mariano D, Di G, Pedro P, Spencer PJ, Pimenta DC. 2018. Biochemical analyses of proteins from duttaphrynus melanostictus (Bufo melanostictus) skin secretion: soluble protein retrieval from a viscous matrix by ion-exchange batch sample preparation. Protein Journal 37(4):380-389

Mariano D, Messias M, Spencer PJ, Pimenta DC. 2019. Protein identification from the parotoid macrogland secretion of Duttaphrynus melanostictus. Journal of Venomous Animals and Toxins Including Tropical Diseases 25:e20190029

Matysik-Woźniak A, Jünemann A, Turski WA, Wnorowski A, Jóźwiak K, Paduch R, Okuno E, Moneta-Wielgoś J, Chorągiewicz T, Maciejewski R, Rejdak R+1 more. 2017. The presence of kynurenine aminotransferases in the human cornea: evidence from bioinformatics analysis of gene expression and immunohistochemical staining. Molecular Vision 23:364-371

Mina AE, Ponti AK, Woodcraft NL, Johnson E, Saporito RA. 2015. Variation in alkaloid-based microbial defenses of the dendrobatid poison frog Oophaga pumilio. Chemoecology 25:169-178

Ministry of Medicine of Japan. 2017. Japanese pharmacopoeia sixteenth edition.

Muramatsu S, Suga Y, Mizuno Y, Hasegawa T, Matsuba S, Hashimoto Y, Volkl A, Seitz J, Ogawa H. 2005. Differentiation-specific localization of catalase and hydrogen peroxide, and their alterations in rat skin exposed to ultraviolet B rays. Journal of Dermatological Science 37:151-158

Okui T, Iwashita M, Rogers MA, Halu A, Atkins SK, Kuraoka S, Abdelhamid I, Higashi H, Ramsaroop A, Aikawa M, Singh SA, Aikawa E+2 more. 2021. CROT (Carnitine O-Octanoyltransferase) is a novel contributing factor in vascular calcification via promoting fatty acid metabolism and mitochondrial dysfunction. Arteriosclerosis, Thrombosis, and Vascular Biology 41:755-768

Patel RK, Mukesh J, Liu Z. 2012. NGS QC toolkit: a toolkit for quality control of next generation sequencing data. PLOS ONE 7:e30619

Qi F, Li A, Inagaki Y, Kokudo N, Tamura S, Nakata M, Tang W. 2011. Antitumor activity of extracts and compounds from the skin of the toad Bufo bufo gargarizans Cantor. International Immunopharmacology 11:342-349

Qi J, Tan CK, Hashimi SM, Zulfiker AH, Good D, Wei MQ. 2014. Toad glandular secretions and skin extractions as anti-inflammatory and anticancer agents. Evidence-Based Complementary and Alternative Medicine 2014:312684

Rolf A, Amos B, Wu CH, Barker WC, Brigitte B, Serenella F, Elisabeth G, Huang H, Rodrigo L, Michele M. 2004. UniProt: the universal protein knowledgebase. Nucleic Acids Research 32:D115-D119

Seppey M, Manni M, Zdobnov EM. 2019. BUSCO: assessing genome assembly and annotation completeness. Methods in Molecular Biology 1962:227-245

State Pharmacopoeia Commission of the PRC. 2020. Pharmacopoeia of the People’s Republic of China. Beijing: People’s Medical Publishing House.

Tharp KM, Khalifeh-Soltani A, Park HM, Yurek DA, Falcon A, Wong L, Feng R, Atabai K, Stahl A. 2016. Prevention of gallbladder hypomotility via FATP2 inhibition protects from lithogenic diet-induced cholelithiasis. American Journal of Physiology-Gastrointestinal and Liver Physiology 310:G855-G64

The Korea Food and Drug Administration,. 2016. Korean pharmacopoeia. Tenth edition monographs.

Ujvari B, Casewell NR, Sunagar K, Arbuckle K, Wüster W, Lo N, O’Meally D, Beckmann C, King GF, Deplazes E, Madsen T+1 more. 2015. Widespread convergence in toxin resistance by predictable molecular evolution. Proceedings of the National Academy of Sciences of the United States of America 112:11911-11916

Valentino A, Calarco A, Di Salle A, Finicelli M, Crispi S, Calogero RA, Riccardo F, Sciarra A, Gentilucci A, Galderisi U, Margarucci S, Peluso G+2 more. 2017. Deregulation of MicroRNAs mediated control of carnitine cycle in prostate cancer: molecular basis and pathophysiological consequences. Oncogene 36:6030-6040

Veglia F, Tyurin VA, Blasi M, De Leo A, Kossenkov AV, Donthireddy L, To TKJ, Schug Z, Basu S, Wang F, Ricciotti E, Di Russo C, Murphy ME, Vonderheide RH, Lieberman PM, Mulligan C, Nam B, Hockstein N, Masters G, Guarino M, Lin C, Nefedova Y, Black P, Kagan VE, Gabrilovich DI+15 more. 2019. Fatty acid transport protein 2 reprograms neutrophils in cancer. Nature 569:73-78

Wang Z, Gerstein M, Snyder M. 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nature Reviews Genetics 10:57-63

Wang YM, Li ZY, Wang JJ, Wu XY, Wang ZM. 2015. Bufadienolides and polyhydroxycholestane derivatives from Bufo bufo gargarizans. Journal of Asian Natural Products Research 17:364-376

Wang J, Lv M, He L, Wang X, Lan Y, Chen J, Chen M, Zhang C, Tang R, Zhou D, Deng X, Li J, Guo T, Price M, Yue B, Fan Z+6 more. 2021. Transcriptomic landscape of persistent diarrhoea in rhesus macaques and comparison with humans and mouse models with inflammatory bowel disease. Gene 800:145837

Wu J-H, Cao Y-T, Pan H-Y, Wang L-H. 2020. Identification of antitumor constituents in toad venom by spectrum-effect relationship analysis and investigation on its pharmacologic mechanism. Molecules 25(18):4269

Wu XY, Gao HM, Wang ZM. 2010. Advances on chemical constituents in crude medicines from toads. Chinese Journal of Experimental Traditional Medical Formulae 16:207-214

Xia F, Gao F, Yao H, Zhang G, Gao B, Lu Y, Wang X, Qian Y. 2019. Identification of angiogenesis-inhibiting peptides from Chan Su. Protein Expression and Purification 163:105445

Yang W, Qi Y, Fu J. 2016. Genetic signals of high-altitude adaptation in amphibians: a comparative transcriptome analysis. BMC Genetics 17:134

Yang W, Qi Y, Lu B, Qiao L, Wu Y, Fu J. 2017. Gene expression variations in high-altitude adaptation: a case study of the Asiatic toad (Bufo gargarizans) BMC Genetics 18:62

Ye J, Zhang Y, Cui H, Liu J, Wu Y, Cheng Y, Xu H, Huang X, Li S, Zhou A. 2018. WEGO 2.0: a web tool for analyzing and plotting GO annotations, 2018 update. Nucleic Acids Research 46(W1):W71-W75

Yenkoyan K, Harutyunyan H, Harutyunyan A. 2018. A certain role of SOD/CAT imbalance in pathogenesis of autism spectrum disorders. Free Radical Biology and Medicine 123:85-95

Yu G, Wang LG, Han Y, He QY. 2012. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics-a Journal of Integrative Biology 16:284-287

Yuki M, Masumi I, Shujiro O, Yoshizawa AC, Minoru K. 2007. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Research 35:W182-W185

Zein J, Gaston B, Bazeley P, De Boer MD, Igo Jr RP, Bleecker ER, Meyers D, Comhair S, Marozkina NV, Cotton C, Patel M, Alyamani M, Xu W, Busse WW, Calhoun WJ, Ortega V, Hawkins GA, Castro M, Chung KF, Fahy JV, Fitzpatrick AM, Israel E, Jarjour NN, Levy B, Mauger DT, Moore WC, Noel P, Peters SP, Teague WG, Wenzel SE, Erzurum SC, Sharifi N+22 more. 2020. HSD3B1 genotype identifies glucocorticoid responsiveness in severe asthma. Proceedings of the National Academy of Sciences of the United States of America 117:2187-2193

Zhan X, Wu H, Wu H, Wang R, Luo C, Gao B, Chen Z, Li Q. 2020. Metabolites from Bufo gargarizans (Cantor, 1842): a review of traditional uses, pharmacological activity, toxicity and quality control. Journal of Ethnopharmacology 246:112178

Zhang Y, Yuan B, Takagi N, Wang H, Zhou Y, Si N, Yang J, Wei X, Zhao H, Bian B. 2019. Comparative analysis of hydrophilic ingredients in toad skin and toad venom using the UHPLC-HR-MS/MS and UPLC-QqQ-MS/MS methods together with the anti-inflammatory evaluation of indolealkylamines. Molecules 24(1):86

Zhou M, Zhang L, Yang Q, Yan C, Jiang P, Lan Y, Wang J, Tang R, He M, Lei G, Sun P, Su N, Price M, Li J, Lin F, Yue B, Fan Z+7 more. 2020. Age-related gene expression and DNA methylation changes in rhesus macaque. Genomics 112:5147-5156

Zhou S-w, Quan J-y, Li Z-w, Ye G, Shang Z, Chen Z-p, Wang L, Li X-y, Zhang X-q, Li J, Liu J-s, Tian H-y+2 more. 2021. Bufadienolides from the eggs of the toad bufo bufo gargarizans and their antimelanoma activities. Journal of Natural Products 84:1425-1433