Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Transcriptome analysis reveals the genetic basis underlying the seasonal development of keratinized nuptial spines in Leptobrachium boringii

BMC Genomics volume 17, Article number: 978 (2016) Cite this article

Abstract

Background

The expression of sexually selected traits often varies with populations’ breeding cycles in many animals. The elucidation of mechanisms underlying the expression of such traits is a research topic in evolutionary biology; however, the genetic basis of the seasonal development of their expression remains unknown. Male Leptobrachium boringii develop keratinized nuptial spines on their upper jaw during the breeding season that fall off when the breeding season ends. To illuminate the genetic basis for the expression of this trait and its seasonal development, we assessed the de novo transcriptome for L. boringii using brain, testis and upper jaw skin and compared gene expression profiles of these tissues between two critical periods of the spine growth cycle.

Results

We identified 94,900 unigenes in our transcriptome. Among them, 2,131 genes were differentially expressed between the breeding period when the spines developed and the post-breeding period when the spines were sloughed. An increased number of differentially expressed genes (DEGs) were identified in the upper jaw skin compared with the testis and brain. In the upper jaw skin, DEGs were mainly enriched in cytosolic part, peptidase inhibitor activity and peptidase regulator activity based on GO enrichment analysis and in glycolysis/gluconeogenesis, ribosome biogenesis in eukaryotes and retinol metabolism based on KEGG enrichment analysis. In the other two tissues, DEGs were primarily involved in the cell cycle, DNA replication and melatonin production. Specifically, insulin/insulin-like growth factor and sex steroid hormone-related DEGs were identified in the upper jaw skin. The expression variation of IGF2 and estrogen-related genes may be the main factors regulating the seasonal development of the spines.

Conclusions

Our study provides a list of potential genes involved in the regulation of seasonal development of nuptial spines in L. boringii. This is the first transcriptome survey of seasonally developed sexually selected traits for non-model amphibian species, and candidate genes provided here may provide valuable information for further studies of L. boringii.

Background

Elucidating the mechanisms underlying the expression of sexually selected traits (SSTs) has been an important research topic in evolutionary biology since Darwin’s theory of sexual selection was developed [1, 2]. Intrasexual selection (competition for mates) and intersexual selection (mate choice) are the two principle mechanisms driving the evolution of animal sexual traits [3]. In most species, males tend to compete with each other for acquiring mates, and females tend to be choosy. Thus, SSTs typically evolved in male individuals [1–3]. The best studied SSTs are elaborate ornaments or weapons, such as the peacock’s tail [4], horns of scarab beetles [5], swords of the swordtail fish [6], and antlers of the deer [7]. Among SSTs, weapons have evolved multiple times across the animal kingdom [8]. Specific to male individuals, weapons are often used in physical combat between male rivals competing for mates by establishing dominance hierarchies or defending reproduction resources [1, 3, 8, 9]. However, despite a wealth of information demonstrating the importance of weapons and their evolution by sexual selection, physiological and genetic mechanisms that regulate their development across animal taxa remain unknown [2, 8].

Within species, characteristics of male SSTs in general and the size of weapons in particular strongly depend on individual body size, physiological condition, nutritional condition and genetic quality [5, 10–14]. Moreover, the expression of diverse male SSTs varies with a populations’ breeding cycle [15–18]. Some male SSTs, such as ornamental plumage of birds, are highly expressed during the breeding season [19]. Other SSTs, such as the nuptial pads of amphibians, nuptial coloration of fishes, or nuptial spines of Leptobrachium anurans, only appear during the breeding season [15–18]. Seasonal changes in plasma sex steroid hormones, especially in testosterone, are important factors that regulate the growth of seasonally developed male sexual traits [16, 18, 20, 21]. The increase in testosterone (or related androgen) in males could trigger the expression of such traits, whereas low levels of this hormone would be insufficient to trigger its activation [18, 19, 22]. Recent studies have demonstrated that the insulin/insulin-like growth factor (IGF) pathway is a widely conserved physiological mechanism regulating the development of pronounced sexual traits across animal taxa [23, 24]. In addition, IGF can interact with sex steroid hormones to regulate the development of condition-sensitive sexual traits [25, 26]. However, prior research on the mechanisms controlling the expression of sexual traits mostly used physiological experiments, such as histology and hormone implants experiments, or measured expression levels of a small number of genes, making the genetic basis of these mechanisms largely unknown [15, 16, 20, 21, 27].

In amphibians, most SSTs develop only during the breeding season [16]. In addition, hormonal-based modulatory mechanisms for the expression of some SSTs, such as skin excrescences at the thumb pads and nuptial pads, have been well studied in numerous anurans [16, 28]. Among amphibians, defensive structures such as spines and tusks are developed by males of few anuran species [29, 30]. They have been phenotypically described in diverse publications, but the physiological mechanisms regulating their development remain unexplored [16, 30]. In Leptobrachium boringii, during the breeding season, adult males develop 10 to 16 keratinized nuptial spines on the edge of their upper jaw that are used as weapons during combat with rivals for territory defense and acquiring mates, and these spines fall off once the breeding season ends (Fig. 1) [31, 32]. Therefore, L. boringii could be an excellent model to study the molecular basis underlying the seasonal development of sexually selected traits for seasonal breeding animals.

Fig. 1
figure 1

Different phenotypes of the keratinized nuptial spines in Leptobrachium boringii during different breeding periods. Male L. boringii grow 10 to 16 nuptial spines on the upper jaw during the breeding period (left), and the spines are sloughed at the post-breeding period (right)

Full size image

High-throughput next-generation sequencing technologies allow for the accurate exploration of gene expression profiles of phenotypic traits for living organisms. Differential gene expression analysis allows studying the genetic basis of context-specific tissue development without the need for a reference genome [33]. Here, we sequenced and assembled the transcriptome of the brain, testis and the upper jaw skin of male L. boringii during the breeding and post-breeding periods using Illumina sequencing technology. We then compared the differential gene expression profiles of the three tissues between these two time periods of the spine growth cycle. Our objective was to identify candidate genes involved in the seasonal development of nuptial spines in Leptobrachium anurans.

Results and Discussion

Illumina sequencing and de novo assembly

Due to the absence of a reference genome for L. boringii, we de novo assembled a transcriptome as a reference for read mapping and gene expression profiling in this species. The assembly generated 602,214,290 raw reads for L. boringii. We obtained an average of 31,817,564 clean reads for our 18 samples (Additional file 1: Table S1). In our assembly, 183,601 transcripts with a total length of 258,254,318 bp were obtained, and they were assembled into 94,900 unigenes, which totaled 81,414,717 bp. N50 lengths of the transcripts and unigenes were 2,941 bp and 1,854 bp, respectively (Table 1), and Q20 of all samples were greater than 96% (Additional file 1: Table S1). In addition, the mapping rates of clean reads to the assembled unigenes for all samples ranged from 75.70% to 81.83% (Additional file 1: Table S1). Taken together, our de novo assemblies revealed a high quality compared with previous studies [5, 34, 35]. This high quality of sequence reads and assembly was the foundation of all our subsequent analyses [36].

Table 1 Number and length of transcripts and unigenes in the Leptobrachium boringii transcriptome
Full size table

Functional annotation and classification of unigenes

We used the databases NR, Swiss-Prot, PFAM, KOG, GO, and KEGG for unigene annotation. NR and Swiss-Prot databases resulted in reliable protein annotations for approximately 25% of unigenes. The number and percentage of annotated unigenes are presented in Table 2. Among all annotated databases, the largest match of annotation hits (24.51%) was obtained with the GO database. KEGG showed the smallest match (11.19%, Table 2). Of 94,900 unigenes assembled, 29,319 genes (30.89%) exhibited a positive match against at least one database. Thus, many of the unigenes have not been functionally annotated. However, several reasons could explain this annotation failure. Those unsuccessfully annotated unigenes could be partially or misassembled transcripts, sequences from UTR protein regions, or meaningless protein genes. In addition, for non-model amphibians, published genomic data in general and especially genetic data on the expression of sexual traits are still lacking. All of these factors combined with the absence of genomic information for L. boringii would have caused the failure of functional annotation of a part of the identified unigenes. However, considering that previous studies successfully annotated approximately 25,000 genes in other non-model amphibian species, such as Odorrana margaretae, Megophrys sangzhiensis and Rhacophorus omeimontis [34, 37], our gene annotation results (more than 29,000 unigenes annotated) from L. boringii’s transcriptome can be considered as high quality.

Table 2 Summary of unigenes annotated in different databases
Full size table

In our study, a total of 23,260 unigenes (~25%) were assigned to 50 sub-categories of GO terms belonging to the following three main categories: cellular component (CC), molecular function (MF) and biological process (BP). These main categories included 17, 11 and 22 sub-categories, respectively (Fig. 2). The most enriched GO terms were related to cell and organelle parts for the category CC, binding and catalytic activity for MF, and cellular process and metabolic process for BP (Fig. 2). The KOG database mainly describes the phylogenetic classification of proteins encoded in complete genomes [38]. In our annotation, only a small proportion of assembled unigenes were mapped to the KOG database, with 12,357 genes classified functionally into 26 categories (Fig. 3). The most enriched KOG categories were General function prediction (2,746 or 22.2%), Signal transduction mechanisms (2,416 or 19.6%), and Posttranslational modification, protein turnover and chaperones (1,158 or 9.4%). In addition, we further annotated all unigenes into pathways in KEGG for understanding their high-level functions and utilities in the biological system of studied animals [39]. A total of 10,625 unigenes were mapped to 259 signaling pathways in our annotation. Signal transduction was the most enriched pathway, followed by pathways in the endocrine system and immune system (Fig. 4). These annotations may provide a valuable resource for further understanding of specific functions and pathways in L. boringii studies.

Fig. 2
figure 2

Gene Ontology (GO) classification of the assembled unigenes in Leptobrachium boringii

Full size image
Fig. 3
figure 3

Eukaryotic orthologous group (KOG) classification of the assembled unigenes in Leptobrachium boringii

Full size image
Fig. 4
figure 4

KEGG classification of the assembled unigenes in Leptobrachium boringii. A to E represent five categories of the KEGG pathways: a cellular processes; b environmental information processing; c genetic information processing; d metabolism; e organismal systems

Full size image

Analysis of differentially expressed genes in the studied tissues

To explore genes controlling the seasonal development of nuptial spines in L. boringii, we compared the gene expression profiles in the brain, testis and upper jaw skin between the breeding and post-breeding periods and identified the significantly differentially expressed genes (DEGs). The final set of DEGs in each tissue was revealed by the hierarchical clustering of expression patterns (Additional file 2: Figure S1). The hierarchical clustering of DEGs revealed that sample groups were more divergent across tissues than between the studied time periods. Each tissue had a relatively distinct gene clusters and expression pattern, so we performed the following differential expression analysis separately for each studied tissue.

Overall, we finally identified 2,131 DEGs in the three tissues, and 75 DEGs were shared by these tissues (Fig. 5d). For approximately 70% (1,484) of the DEGs that were functionally annotated, we were able to complete an overall comprehensive review of these DEGs. As shown in the Venn diagram (Fig. 5d), more DEGs were found in the upper jaw skin compared with the brain and testis. This finding suggests that the upper jaw skin causes greater changes in gene expression than the sexual glands (testes) and central nervous system (brain) during different breeding periods in L. boringii, indicating that the upper jaw skin might play a relatively important role in the development of nuptial spines compared with the other two tissues. Significant GO enrichment and KEGG enrichment (p < 0.05) in each tissue could provide valuable information to understand the gene function of the DEGs involved in the development of nuptial spines (Table 3; Table 4; Additional file 1: Table S2).

Fig. 5
figure 5

Differential expression analysis of tissue comparisons between the breeding and post-breeding periods. Volcano plot of differentially expressed genes in the brain (a), testis (b) and upper jaw skin (c; spine for short in figure) are separately presented. The X axis represents differential fold-change (log2FC) of the differentially expressed genes (DEGs), and the Y axis represents the –log padj value of the DEGs. padj represents the adjusted P-value of the DEG, and 0.05 was set as the significance level. Significantly up-regulated DEGs in the breeding period compared with the post-breeding period are presented in red, whereas significantly down-regulated DEGs are shown in green. Non-significantly differentially expressed genes are presented in blue. d Venn diagram showing co-differentially expressed genes among different tissue comparisons in Leptobrachium boringii

Full size image
Table 3 The most enriched GO terms of differentially expressed genes in the three tissues
Full size table
Table 4 Significantly enriched KEGG pathways among differentially expressed genes in Leptobrachium boringii
Full size table

To validate the differentially expression profiles obtained from RNA-seq, we selected eight DEGs among the three tissues for qRT-PCR analysis. Six of them closely matched the results detected by RNA-seq, with a correlation coefficient of R = 0.957 (Fig. 6). Although the other two genes (GAG1L, KRT) did not have a perfect match and exhibited approximately 4-fold differences in relative expression, they still exhibited the same expression pattern at the two time periods. The correlation coefficient for all eight genes was R = 0.745. Overall, the qRT-PCR results were consistent with the RNA-sequencing data, indicating that our transcriptome data were credible.

Fig. 6
figure 6

Quantitative real-time PCR confirmation of differentially expressed genes identified by RNA-seq. Relative expression levels are fold-changes expressed as the ratio of gene expression at the breeding period to post-breeding period normalized with GADPH gene

Full size image

Differentially expressed genes in the upper jaw skin

A total of 1,181 DEGs were identified in the upper jaw skin between the breeding and post-breeding periods in L. boringii. Among them, 748 genes (63.34%) were down-regulated, and the remaining 433 genes (36.66%) were up-regulated during the breeding period (Fig. 5c). A complete list of the DEGs in the upper jaw skin is provided (Additional file 1: Table S3). To better understand the biological processes of the DEGs in the skin during the two time periods of the development cycle for nuptial spines, we assigned all DEGs to GO and KEGG databases and analyzed the significant enrichments. However, despite the fact that so many DEGs were identified from the upper jaw skin, no significant enrichment of GO terms was identified among all DEGs (Additional file 1: Table S2). In up-regulated DEGs in the upper jaw skin, ncRNA metabolic process, tRNA metabolic process and ether hydrolase activity were found to be significantly enriched (adjust P < 0.05), indicating that the nuptial spines were actively transcribed during the growth period (Additional file 1: Table S2). The most enriched GO terms among DEGs (P < 0.05) in the skin were related to the cytosolic part, peptidase inhibitor activity, peptidase regulator activity and the extracellular region (Table 3a). This finding indicates that the seasonal development of the nuptial spines may be caused by changes in the cytosol cells and their extracellular matrix from the ‘dermis base’ of spines in upper jaw skin. The differential expression of peptidase inhibitor and regulator activities may be involved in the regulation of growth metabolism for the nuptial spines. The significantly enriched KEGG pathway among all DEGs in the upper jaw skin was only glycolysis/gluconeogenesis (Table 4). In addition, both the up-regulated and down-regulated genes were significantly enriched in the pathways related to the category of metabolism and genetic information processing. Yet, different biological processes were enriched between up-regulated and down-regulated genes (Table 4). During the breeding period, DEGs involved in translation and protein processing exhibited high expression levels, whereas during the post-breeding period, DEGs were mainly involved in the metabolism of cofactors and vitamins, amino acids and carbohydrates.

We identified the top 10 up-regulated and 10 down-regulated genes in the upper jaw skin by the order of FDR correction to obtain potential genes that may directly control the seasonal development of the nuptial spines. Our results indicate that the most differentially expressed genes in the upper jaw skin were mainly expressed in the extracellular matrix and were involved in peptidase inhibitor activity, the regulation of transcription and metabolic processes including COL28A1, SERPINB4, OGN and CSTB (Table 5). The top three most down-regulated genes during the breeding period were all lysosomal alpha-gluacosidase-like genes (Table 5). A previous study proposed that lysosome-mediated autophagy may participate in the regulation of cell death by degrading the extracellular matrix [40]. We thus considered that the epidermal cells of the upper jaw skin connected to the keratinized spines may undergo apoptosis when these spines are sloughed. Newt-specific cysteine-rich growth regulator (nsCCN) is a growth factor that was first identified in the newt species Notophthalmus viridescens and has an important role in the regeneration process of heart tissue [41]. In L. boringii, we observed a high expression level of nsCCN in the upper jaw skin after the spines had sloughed. Hence, we hypothesized that nsCCN may play a crucial role in the regeneration of epidermal cells after spines separate from the skin. The surface of the nuptial spines is highly keratinized, and the deciduous part is mainly the keratinized region. Interestingly, many of the keratin-related genes were up-regulated in the upper jaw skin when these spines developed. These genes were down-regulated when these weapons were sloughed (Additional file 1: Table S3). These findings support the involvement of keratin-related genes in the development of nuptial spines.

Table 5 Top 20 most differentially expressed genes in the upper jaw skin
Full size table

Differentially expressed genes in the testis and brain

Using the same method as for the upper jaw skin, we further analyzed DEGs in the testis and brain. In summary, fewer DEGs were identified in these two tissues compared with the upper jaw skin. A total of 843 and 107 genes were differentially expressed in the testis and brain, respectively. In the testis, 499 genes (59.19%) were down-regulated and 344 genes (40.81%) were up-regulated during the breeding period (Fig. 5b; Additional file 1: Table S4). In the brain, the corresponding data were 52 genes (48.60%) and 55 genes (51.40%), respectively (Fig. 5a; Additional file 1: Table S5).

Among the three studied tissues, we found the greatest elevated number of enriched GO terms in the testis with 16 terms significantly enriched among all DEGs and 103 GO terms significantly enriched among down-regulated genes. The enriched GO terms were predominantly related to chromosome, DNA replication metabolism, cell division, cell cycle and spermatogenesis (Table 3b, Additional file 1: Table S2). Consistent with the GO term enrichment results, significant pathways enriched in KEGG were mostly assigned to genetic information processing, cellular processes, and metabolism (Table 4). Cell cycle, DNA replication, Cell cycle-yeast, and Meiosis-yeast were the most enriched pathways among all DEGs. The top most up-regulated and down-regulated genes identified in the testis were mainly related to the above function categories (Additional file 1: Table S6). Overall, the enrichment results overwhelmingly highlighted the activity of genetic processing and cellular processes in the testis. Generally, the seasonal expression of secondary sexual traits and sex accessory structures in most amphibians, especially the nuptial excrescences of anurans, are regulated by gonadal activity and gonadal steroid hormones [16]. However, we identified no significantly enriched terms related to steroid hormones in the testis. In L. boringii, time spent on parental care for offspring by reproducing males represents a significant proportion (greater than 70%) of the breeding season [32, 42]. Such a parental behavior is often correlated with low androgen levels in vertebrate males [21, 43]. Hence, the steroid hormone-related genes were rarely differentially expressed in the testis during the breeding season of L. boringii. We still find some DEGs in the testis that might be indirectly involved in steroid hormone activity. The relaxin and progesterone receptor genes are good examples of such genes (Additional file 1: Table S4). Therefore, we hypothesized that the gonadal activity may not directly participate in the regulation of the seasonal development of the nuptial spines through testosterone but though other factors such as estrogen-related genes.

Probably due to the limited number of DEGs in the brain, no significant enrichment of GO terms and KEGG pathways was identified. Among the enriched GO terms (p < 0.05), “aralkylamine N-acetyltransferase activity” was the most enriched term (Table 3c). Aralkylamine N-acetyltransferase is involved in the production of melatonin; it functions in the modulation of circadian rhythm and is essential for seasonal reproduction of vertebrates [44]. The top DEGs in the brain were mainly related to signal transduction (Additional file 1: Table S6). Among the top DEGs, kisspeptin precursor gene was the most down-regulated during the breeding period (Table 3). Kisspeptin is a G-protein coupled receptor ligand that can stimulate the secretion of aldosterone and the release of insulin [45]. Furthermore, we found a consistent expression pattern of ncCCN in the brain as well as in the upper jaw skin with ncCCN exhibiting high expression levels during the post-breeding period. Previous studies reported that the phenotypic variation of sexual traits would lead to alterations of gene expression in the brain [46]. Thus, we proposed that the brain may indirectly modulate the seasonal expression of nuptial spines through the genes specified above but also through other still unidentified genes.

Potential genes regulating the seasonal development of the nuptial spines

Previous studies have proposed that the development and maintenance of male sexually selected traits in vertebrates often rely on the regulation of androgens, including testosterone [16, 21, 47]. In addition, the expression of SSTs tends to be more sensitive to the body condition than to other phenotypic traits of the holders, and recent discoveries implicated the widely conserved insulin/insulin-like growth factor pathway as the common physiological mechanism regulating the growth of these sexual traits [23, 24]. Based on our results and results of previous studies, we further selected DEGs with functional categories related to insulin/insulin-like growth factor and sex steroid hormone-related genes from all DEGs between the two studied time periods to identify potential genes regulating the seasonal development of the nuptial spines.

Insulin/insulin-like growth factor-related genes

As a result, we selected three significant DEGs between the two studied periods associated with insulin/insulin-like growth factor (IGFBP2, IGF2-B, nsCCN) in the upper jaw skin. Consistent with previous studies, IGFBP2 and IGF2-B revealed up-regulated expression when the nuptial spines were developed (Table 6). The insulin-like growth factors (IGF-I and IGF-II) are produced by most tissues of vertebrates and play roles in the regulation of proliferation and differentiation of cells in most tissues [26]. IGF-binding proteins (IGFBP) alter the interaction of IGFs with their cell surface receptors and modulate the function of the IGFs in cell culture [48]. The up-regulated IGF2-B and IGFBP2 genes during the breeding period may increase the metabolism and growth of the cells of the upper jaw skin connected with the nuptial spines, allowing the growth of these weapons. However, the exact molecular pathway involved must be verified by further experiments. Additionally, the up-regulated nsCCN in the upper jaw skin and brain when nuptial spines are sloughed may play an important role given their vital role in the tissue regeneration process.

Table 6 Significantly differentially expressed genes involved in the two candidate pathways in Leptobrachium boringii
Full size table

The expression of striking sexual traits is always associated with the holder’s body condition, as their expression is commonly sensitive to the nutrition history [24, 49]. However, prior studies on condition-dependent sexual traits rarely focused on amphibians, and such sexual traits are often tightly coupled with individual variation within a population [25]. In our study, for the first time, we found that the expression of a striking male sexual trait was likely associated with the changes in physiological condition during a reproductive cycle of an amphibian species. Insulin-like growth factor 2 (IGF2) may be the main peptide affecting the development of nuptial spines in L. boringii, and this is different from the regulation pattern for exaggerated traits in most vertebrates through IGF1 [23, 27]. In conclusion, the seasonal development of the nuptial spines in L. boringii may operate in a condition-dependent manner, and the IGF2 gene may be involved in the regulation of this development.

Sex steroid hormone-related genes

Five significantly DEGs involved in sex steroid hormone-related pathways were identified in the upper jaw skin, and two related DEGs were identified in the testis (Table 6). Seasonally breeding amphibians often undergo dramatic changes in the concentration of sex hormones and in reproductive characteristics such as sexual morphological traits and reproductive behaviors triggered by those hormones at the onset and termination of each breeding cycle [16, 50]. Specifically, 17β-hydroxysteroid dehydrogenases (17β-HSD) are a group of alcohol oxidoreductases that play a key role in the regulation of intracellular concentrations of steroid hormones, including androgens and estrogens [51, 52]. CYP3A5 encodes a member of the cytochrome P450 enzyme superfamily catalyzing the metabolism of steroid hormone into less biologically active metabolites [53, 54]. During the post-breeding period, the high expression of 17β-HSD genes (HSD17B2, HSD17B6 and HSD17B11) and CYP3A5 may increase the degradation and metabolism of steroid hormones and thus decrease the concentration of active androgen and estrogen in the upper jaw skin. Generally, testosterone regulates the expression of male secondary sexual traits directly by stimulating intracellular androgen receptors in the target tissues or indirectly by aromatizing androgens into estrogen and binding to estrogen receptors [21]. Interestingly, in the upper jaw skin, no androgen receptor but an estrogen receptor gene (NR3A1) was differentially expressed between the breeding and the post-breeding periods (Table 6). Once activated, NR3A1 (or ERα) can translocate into the nucleus and regulate the activity of eukaryotic gene expression, cellular proliferation and differentiation in target tissues [55, 56]. However, to date, only a few amphibian species have been reported to have α isoform of the estrogen receptor (ER) [57, 58]. Our study is the first to report the expression of ERα in the skin of an amphibian species, indicating that the upper jaw skin where the spines develop may be a target tissue whose growth and replacement are regulated by estrogen.

Commonly, the presence of male parental care behavior will decrease the expression of secondary sexual traits [28]. However, male L. boringii express and maintain nuptial spines, a secondary sexual trait, during the entire breeding season while simultaneously taking care of offspring for a relatively long time period (at least two months) [31, 42]. Together with the low androgen-dependent parental care behavior in L. boringii, the seasonal acquisition and maintenance of nuptial spines may be regulated by estrogen. Moreover, the up-regulated 7β-HSD genes and ERα would have an inhibitory effect on androgen output during the post-breeding period [59, 60]. Thus, estrogens may play an important role in the seasonal regulation of the nuptial spines in L. boringii. However, elucidating whether the seasonal development of the nuptial spines is regulated by androgen directly, by aromatization of androgen into estrogens exerting related effects, or by both still needs further investigation.

Relaxin is a small peptide hormone that can enhance male sperm motility and fertilization capacity in vertebrates [61, 62]. The expression of fRLX was altered with the annual reproductive cycle of a frog species Pelophylax esculentus and regulated by testosterone [63, 64]. The up-regulated fRLX in the testis of L. boringii during the breeding period may enhance male reproductive activity. The progesterone receptor (PGR) gene was also up-regulated in the testis during the breeding period. This gene is activated by progesterone and can inhibit androgen-dependent reproductive behaviors in males [65, 66]. Thus, we hypothesize that the high expression of PGR during the breeding season may increase the propensity of males for providing parental care at the expense of courtship and mating activities. Finally, there may be other unknown important genes participating in the regulation of seasonal expression for this trait in L. boringii; thus, further study is needed.

Taken together, candidate genes involved in insulin/insulin-like growth factor and sex steroid hormone pathways indicate that the nuptial spines in L. boringii exhibit high estrogen sensitivity and low insulin/insulin-like growth factor sensitivity during the post-breeding period when the spines are sloughed. Thus, the seasonal development and maintenance of the nuptial spines in male L. boringii are likely modulated by the circulating sex steroid hormones and the nutritional status. Candidate genes in the testis and brain may indirectly participate in this process by regulating the synthesis and metabolism activity of steroid hormones and insulin growth factor. Further study should focus on these two candidate pathways and use testable experiments to validate the molecular mechanism.

Conclusions

In this study, we conducted a transcriptome analysis of an anuran species to explore the genetic basis of the seasonal growth of a sexual weapon trait. Our overall transcriptome provided a rich list of unigenes expressed in adult males of L. boringii during the breeding and the post-breeding periods. DEGs between different growth periods of the keratinized nuptial spines in L. boringii were identified. In total, 94,900 unigenes and 2,131 DEGs were generated in our transcriptome. These genes greatly contribute to current genetic resources for the anuran species. The number of DEGs identified from the upper jaw skin was larger than the testis and the brain, indicating that the upper jaw skin might contribute more to the seasonal changes in the spine activity. Candidate genes related to steroid hormone signaling pathways and the insulin/insulin-like growth factor pathway are likely involved in regulating the seasonal growth of the spines. The candidate genes and pathways are promising for future studies on the seasonal growth of the nuptial spines in Leptobrachium species and may be helpful for other vertebrates.

Methods

Sample collection and RNA extraction

L. boringii male individuals were caught from the population of the Badagong Mountain Natural Reserve (Hunan, China, 29°39′-29°49′ N, 109°41′-110°09′ E) during the breeding season of 2014. Sampling was performed during the breeding period (in late March when the nuptial spines were developed) and during the subsequent post-breeding period (in late May once the spines had sloughed, Fig. 1). At each sampling time, three individuals were caught as biological replicates for each study period. Once caught, these toads were kept in single plastic boxes containing the same plant leaves and water from their natural environment and transported to the laboratory. The living toads were euthanized by injecting buffered MS-222 solution (3 g/l) through the dorsal skin. Each toad was dissected with 0.1% diethyl pyrocarbonate (DEPC)-ddH2O solution-treated scissors, and a tissue sample (each approximately 100 mg in scale) was taken from the brain, testis and the upper jaw skin and stored in liquid nitrogen until needed. Specifically, the tissue upper jaw skin was taken from the skin area that contained a ‘papillary dermis base’ of the spine, and the tissue brain was taken from the whole brain region of each sampled individual. The testis sample was randomly taken from the whole tissue with the certain portion of each individual.

Total RNA of each sample was extracted using an RNA extraction kit (Omega Bio-Tek, USA) with TRIzol® Reagent (Invitrogen, CA, USA) following the manufacturer’s instructions. RNA degradation and contamination was tested using 1% agarose gels. We measured RNA purity and RNA concentration using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA), respectively. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Total RNA was extracted from each tissue sample separately for three replicates of each breeding period, and 18 high-quality RNA samples were finally acquired.

Illumina sequencing and de novo assembly

For each tissue, using 3 μg total RNA as input material, a cDNA library was constructed using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s instructions. In brief, mRNA was purified from total RNA using poly-T oligo (dT) magnetic beads, and fragments were generated using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNase H−), and second strand cDNA was subsequently synthesized with DNA Polymerase I and RNase H. The library fragments with a preferential length of 150 to 200 bp were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then, PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer. We purified PCR products using the AMPure XP system and assessed library quality using the Agilent Bioanalyzer 2100 system. The above cDNA samples were then clustered using the TruSeq PE Cluster Kit (Illumina, USA) on a cBot Cluster Generation System. We sequenced all libraries on an Illumina Hiseq 2000 platform (San Diego, CA, USA), and 100 bp paired-end reads were then generated for each of them.

Clean reads were obtained by removing the trimming adapter sequences, reads containing poly-N and low quality reads (reads with greater than 5% unknown nucleotides or reads with less than 13 bp) from raw reads. The number of raw reads, clean reads and Q20 were calculated at the same time. All subsequent analyses were based on the clean reads. Due to the absence of reference sequences for L. boringii, high quality clean reads were de novo assembled using Trinity software (r20140717) with the min_kmer_cov set to 2 and with default values for all remaining parameters [67]. Unigenes was selected from the longest transcript copy of each gene clusters to avoid redundant transcripts [68]. All unigenes from the 18 tissue samples were combined in a unigene database for the studied species, and all subsequent analyses were performed on this database.

Functional annotation of unigenes

To functionally annotate the unigenes in L. boringii, we searched our unigene set against the datasets of non-redundant protein sequences (NR) in the National Center for Biotechnology Information (NCBI; http://www.ncbi.nlm.nih.gov), the Protein family database (PFAM; http://pfam.xfam.org/) and the manually annotated and reviewed protein sequence database (Swiss-Prot database; http://www.ebi.ac.uk/swissprot/) using NCBI-BLAST (v 2.2.30) [69] with an E-value cut-off of 1 × 10−5. Annotation was extended to the Functional classifications of Gene Ontology (GO; http://www.geneontology.org/), the Clusters of Orthologous Groups of proteins (KOG/COG; http://www.ncbi.nlm.nih.gov/COG), and Pathway Annotation of the Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/pathway.html). GO annotation was performed using Blast2GO pipelines with an E-value ≤ 1E-3 [70, 71]. KOG annotation was searched with an E-value ≤1E-5 [38]. We searched the unigenes against the KEGG database using KOBAS (v 2.0) [72].

Differential expression analyses

Gene expression levels in each tissue were estimated using the software package RSEM (v 1.2.0) based on the expectation-maximization (EM) algorithm [73]. Firstly, clean reads of each sample were mapped back onto the assembled transcripts. Then, the read count for each gene was obtained from the mapping results. The gene expression level was estimated by calculating the fragments per kb per million reads (FPKM) [74]. The profiles of gene expression level generated for each tissue were combined and used to detect DEGs between the breeding and the post-breeding periods. We used DESeq R package (v2.15.3) [75] to perform differential gene expression analyses in the three tissues between these two time periods. The resulting p-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate [76]. Genes with an adjusted p-value (FDR) < 0.05 and |log2 (fold-change)| > 1 were assigned as significant level.

Enrichment analysis of differentially expressed genes

GO enrichment of the DEGs was analyzed using GOseq R packages based on Wallenius non-central hyper-geometric distribution [77], which can adjust for gene length bias in DEGs. The smallest P-value indicates the highest degree of GO enrichment. KEGG pathway enrichment of the DEGs was tested using KOBAS software [72]. The hypergeometric test was applied in the enrichment analysis. The minimum threshold of statistical significance was fixed at p = 0.05.

Data confirmation by quantitative PCR analysis

To confirm the accuracy of our transcriptome data, we selected eight genes among the significantly DEGs in the three tissues and performed quantitative PCR experiments. PCRs were performed with a CFX96 Touch Real-Time PCR Detection System (BIO-RAD, Hercules, CA, USA) using SYBR Green I Dye with a 20-μl total volume mixture containing 2 μl of cDNA as the template. The expression levels of the tested genes were normalized with a housekeeping gene (glyceraldehyde-3-phosphate dehydrogenase, GAPDH) as an internal control, and relative gene expression was analyzed with the 2-△△Ct method [78]. Mean values and standard errors were determined by three biological replicates. The specific primers for these genes were designed using PRIMER 5, and information on these primers is additionally provided (Additional file 1: Table S7).

Abbreviations

CYP3A5:

Cytochrome P450 Family 3 subfamily a member 5

DEG:

Differentially expressed genes

FDR:

False discovery rate adjust p value

GO:

Gene ontology

HSD17B:

17-beta-hydroxysteroid dehydrogenase

IGF:

Insulin/insulin-like growth factor

IGFBP:

Insulin-like growth factor-binding protein

KEGG:

Kyoto encyclopedia of genes and genomes

KOG/COG:

Clusters of orthologous groups of proteins

NR:

Non-redundant protein sequences

NR3A1:

Estrogen receptor 1

nsCCN:

Newt-specific cysteine-rich growth regulator

PGR:

Progesterone receptor

PFAM:

Protein family database

qRT-PCR:

Quantitative reverse transcription real-time PCR

SSTs:

Sexually selected traits

Swiss-Prot:

The manually annotated and reviewed protein sequence database.

References

  1. Darwin C. The descent of man, 2 Vols. London. 1871;81:130–1.

    Google Scholar 

  2. Cornwallis CK, Uller T. Towards an evolutionary ecology of sexual traits. Trends Ecol Evol. 2010;25:145–52.

    Article  PubMed  Google Scholar 

  3. Andersson MB. Sexual selection. Princeton: Princeton University P ress; 1994.

  4. Ohlsson T, Smith HG, Råberg L, Hasselquist D. Pheasant sexual ornaments reflect nutritional conditions during early growth. Proc R Soc of B. 2002;269:21–7.

    Article  Google Scholar 

  5. Warren IA, Vera JC, Johns A, Zinna R, Marden JH, Emlen DJ, et al. Insights into the development and evolution of exaggerated traits using de novo transcriptomes of two species of horned scarab beetles. PloS One. 2014;9:e88364.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Kang JH, Manousaki T, Franchini P, Kneitz S, Schartl M, Meyer A. Transcriptomics of two evolutionary novelties: how to make a sperm-transfer organ out of an anal fin and a sexually selected "sword" out of a caudal fin. Ecol Evol. 2015;5:848–64.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Malo AF, Roldan ER, Garde J, Soler AJ, Gomendio M. Antlers honestly advertise sperm production and quality. Proc R Soc of B. 2005;272:149–57.

    Article  Google Scholar 

  8. Emlen DJ. The Evolution of Animal Weapons. Annu Rev Ecol Evol Syst. 2008;39:387–413.

    Article  Google Scholar 

  9. Emlen DJ, Marangelo J, Ball B, Cunningham CW. Diversity in the weapons of sexual selection: horn evolution in the beetle genus Onthophagus (Coleoptera: Scarabaeidae). Evolution. 2005;59:1060–84.

    Article  CAS  PubMed  Google Scholar 

  10. David P, Bjorksten T, Fowler K, Pomiankowski A. Condition-dependent signalling of genetic variation in stalk-eyed flies. Nature. 2000;406:186–8.

    Article  CAS  PubMed  Google Scholar 

  11. Cotton S, Fowler K, Pomiankowski A. Condition dependence of sexual ornament size and variation in the stalk‐eyed fly Cyrtodiopsis dalmanni (Diptera: Diopsidae). Evolution. 2004;58:1038–46.

    Article  PubMed  Google Scholar 

  12. Kijimoto T, Snell-Rood EC, Pespeni MH, Rocha G, Kafadar K, Moczek AP. The nutritionally responsive transcriptome of the polyphenic beetle Onthophagus taurus and the importance of sexual dimorphism and body region. Proc R Soc B. 2014;281:20142084.

  13. Faivre B, Grégoire A, Préault M, Cézilly F, Sorci G. Immune activation rapidly mirrored in a secondary sexual trait. Science. 2003;300:103.

    Article  CAS  PubMed  Google Scholar 

  14. Thompson CW, Hillgarth N, Leu M, McClure HE. High parasite load in house finches (Carpodacus mexicanus) is correlated with reduced expression of a sexually selected trait. Am Nat. 1997;149:270–94.

    Article  Google Scholar 

  15. Svensson PA, Pelabon C, Blount JD, Forsgren E, Bjerkeng B, Amundsen T. Temporal variability in a multicomponent trait: nuptial coloration of female two-spotted gobies. Behav Ecol. 2008;20:346–53.

    Article  Google Scholar 

  16. Sever DM, Staub NL. Hormones, sex accessory structures, and secondary sexual characteristics in amphibians.In: Norris DO editior.Hormones and Reproduction in Vertebrates Volume 2: Amphibians. Maryland Heights: Elsevier Inc; 2011.

  17. Zheng Y, Li S, Fu J. A phylogenetic analysis of the frog genera Vibrissaphora and Leptobrachium, and the correlated evolution of nuptial spine and reversed sexual size dimorphism. Mol Phylogenet Evol. 2008;46:695–707.

    Article  CAS  PubMed  Google Scholar 

  18. Loree’A H, Propper CR. Effects of androgens on male sexual behavior and secondary sex characters in the explosively breeding spadefoot toad, Scaphiopus couchii. Horm Behav. 1997;31:89–96.

    Article  Google Scholar 

  19. Lindsay WR, Webster MS, Schwabl H. Sexually selected male plumage color is testosterone dependent in a tropical passerine bird, the red-backed fairy-wren (Malurus melanocephalus). PloS One. 2011;6:e26067.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Peters A, Astheimer LB, Boland CR, Cockburn A. Testosterone is involved in acquisition and maintenance of sexually selected male plumage in superb fairy-wrens, Malurus cyaneus. Behav Ecol Sociobiol. 2000;47:438–45.

    Article  Google Scholar 

  21. Hau M. Regulation of male traits by testosterone: implications for the evolution of vertebrate life histories. BioEssays. 2007;29:133–44.

    Article  CAS  PubMed  Google Scholar 

  22. Gonzalez G, Sorci G, Smith LC, Lope F. Testosterone and sexual signalling in male house sparrows (Passer domesticus). Behav Ecol Sociobiol. 2001;50:557–62.

    Google Scholar 

  23. Emlen DJ, Warren IA, Johns A, Dworkin I, Lavine LC. A mechanism of extreme growth and reliable signaling in sexually selected ornaments and weapons. Science. 2012;337:860–4.

    Article  CAS  PubMed  Google Scholar 

  24. Warren IA, Gotoh H, Dworkin IM, Emlen DJ, Lavine LC. A general mechanism for conditional expression of exaggerated sexually-selected traits. BioEssays. 2013;35:889–99.

    PubMed  Google Scholar 

  25. Dantzer B, Swanson EM. Mediation of vertebrate life histories via insulin-like growth factor-1. Biol Rev Camb Philos Soc. 2012;87:414–29.

    Article  PubMed  Google Scholar 

  26. Cohick W, Clemmons DR. The insulin-like growth factors. Annu Rev Physiolo. 1993;55:131–53.

    Article  CAS  Google Scholar 

  27. Shingleton AW, Frankino WA. New perspectives on the evolution of exaggerated traits. BioEssays. 2013;35:100–7.

    Article  PubMed  Google Scholar 

  28. Emerson SB, Carroll L, Hess DL. Hormonal induction of thumb pads and the evolution of secondary sexual characteristics of the Southeast Asian fanged frog, Rana blythii. J Exp Biol. 1997;279:587–96.

    CAS  Google Scholar 

  29. Shine R. Sexual selection and sexual dimorphism in the Amphibia. Copeia. 1979;2:297–306.

  30. Duellman WE, Trueb L. Biology of amphibians. Baltimore: John Hopkins University press; 1994.

  31. Cameron MH, Xianjin HE, Jinzhong FU. Keratinized Nuptial Spines Are Used for Male Combat in the Emei Moustache Toad (Leptobrachium boringii). Asian Herpetol Res. 2011;2:142–8.

    Article  Google Scholar 

  32. Hudson CM, Fu J. Male-biased sexual size dimorphism, resource defense polygyny, and multiple paternity in the Emei moustache toad (Leptobrachium boringii). PloS One. 2013;8:e67502.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet. 2011;12:87–98.

    Article  CAS  PubMed  Google Scholar 

  34. Qiao L, Yang W, Fu J, Song Z. Transcriptome profile of the green odorous frog (Odorrana margaretae). PloS One. 2013;8:e75211.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Shao Y, Wang L-J, Zhong L, Hong M-L, Chen H-M, Murphy RW, et al. Transcriptomes reveal the genetic mechanisms underlying ionic regulatory adaptations to salt in the crab-eating frog. Sci Rep. 2015;5:17551.

  36. Robertson G, Schein J, Chiu R, Corbett R, Field M, Jackman SD, et al. De novo assembly and analysis of RNA-seq data. Nat Methods. 2010;7:909–12.

    Article  CAS  PubMed  Google Scholar 

  37. Huang L, Li J, Anboukaria H, Luo Z, Zhao M, Wu H. Comparative transcriptome analyses of seven anurans reveal functions and adaptations of amphibian skin. Sci Rep. 2016;6:24069.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Tatusov RL, Fedorova ND, Jackson JD, Jacobs AR, Kiryutin B, Koonin EV, et al. The COG database: an updated version includes eukaryotes. BMC Bioinformatics. 2003;4:1.

    Article  Google Scholar 

  39. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:D480–4.

    Article  CAS  PubMed  Google Scholar 

  40. Levine B, Klionsky DJ. Development by self-digestion: molecular mechanisms and biological functions of autophagy. Dev Cell. 2004;6:463–77.

    Article  CAS  PubMed  Google Scholar 

  41. Looso M, Michel CS, Konzer A, Bruckskotten M, Borchardt T, Krüger M, et al. Spiked-in pulsed in vivo labeling identifies a new member of the CCN family in regenerating newt hearts. J Proteome Res. 2012;11:4693–704.

    Article  CAS  PubMed  Google Scholar 

  42. Zheng Y, Deng D, Li S, Fu J. Aspects of the breeding biology of the Omei mustache toad (Leptobrachium boringii): polygamy and paternal care. Amphibia-Reptilia. 2010;31:183–94.

    Article  Google Scholar 

  43. Wingfield JC, Hegner RE, Dufty Jr AM, Ball GF. The "challenge hypothesis": theoretical implications for patterns of testosterone secretion, mating systems, and breeding strategies. Am Nat. 1990;136:829–46.

  44. Singh M, Jadhav HR. Melatonin: functions and ligands. Drug Discov Today. 2014;19:1410–8.

    Article  CAS  PubMed  Google Scholar 

  45. Pasquier J, Kamech N, Lafont A-G, Vaudry H, Rousseau K, Dufour S. Molecular evolution of GPCRs: Kisspeptin/kisspeptin receptors. J Mol Endocrinol. 2014;52:T101–17.

    Article  CAS  PubMed  Google Scholar 

  46. Jašarević E, Geary DC, Rosenfeld CS. Sexually selected traits: a fundamental framework for studies on behavioral epigenetics. ILAR J. 2012;53:253–69.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Dufty AM, Clobert J, Møller AP. Hormones, developmental plasticity and adaptation. Trends Ecol Evol. 2002;17:190–6.

    Article  Google Scholar 

  48. Jones JI, Clemmons DR. Insulin-Like Growth Factors and Their Binding Proteins: Biological Actions. Endocr Rev. 1995;16:3–34.

    CAS  PubMed  Google Scholar 

  49. Andersson M. Evolution of condition-dependent sex ornaments and mating preferences: sexual selection based on viability differences. Evolution. 1986;40:804–16.

  50. Wilczynski W, Lynch KS, O'Bryant EL. Current research in amphibians: studies integrating endocrinology, behavior, and neurobiology. Horm Behav. 2005;48:440–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Labrie F, Luu-The V, Lin S-X, Claude L, Simard J, Breton R, et al. The key role of 17β-hydroxysteroid dehydrogenases in sex steroid biology. Steroids. 1997;62:148–58.

    Article  CAS  PubMed  Google Scholar 

  52. Baker ME. Evolution of 17β-hydroxysteroid dehydrogenases and their role in androgen, estrogen and retinoid action. Mol Cell Endocrino. 2001;171:211–5.

    Article  CAS  Google Scholar 

  53. Waxman DJ, Lapenson DP, Aoyama T, Gelboin HV, Gonzalez FJ, Korzekwa K. Steroid hormone hydroxylase specificities of eleven cDNA-expressed human cytochrome P450s. Arch Biochem Biophys. 1991;290:160–6.

    Article  CAS  PubMed  Google Scholar 

  54. Tsuchiya Y, Nakajima M, Yokoi T. Cytochrome P450-mediated metabolism of estrogens and its regulation in human. Cancer Lett. 2005;227:115–24.

    Article  CAS  PubMed  Google Scholar 

  55. Hess R, Bunick D, Lubahn DB, Zhou Q, Bouma J. Morphologic changes in efferent ductules and epididymis in estrogen receptor-a knockout mice. J Androl. 2000;21:107–21.

    CAS  PubMed  Google Scholar 

  56. Bjornstrom L, Sjoberg M. Mechanisms of estrogen receptor signaling: convergence of genomic and nongenomic actions on target genes. Mol Endocrinol. 2005;19:833–42.

    Article  PubMed  Google Scholar 

  57. Wu KH, Tobias ML, Thornton JW, Kelley DB. Estrogen receptors in Xenopus: duplicate genes, splice variants, and tissue-specific expression. Gen Comp Endocr. 2003;133:38–49.

    Article  CAS  PubMed  Google Scholar 

  58. Takase M, Iguchi T. Molecular cloning of two isoforms of Xenopus (Silurana) tropicalis estrogen receptor mRNA and their expression during development. Biochim Biophys Acta Gene Struct Expression. 2007;1769:172–81.

    Article  CAS  Google Scholar 

  59. Fasano S, Pierantoni R, Minucci S, Di Matteo L, Chieffi G. Seasonal fluctuations of estrogen-binding activity in the testis of the frog, Rana esculenta. Gen Comp Endocr. 1989;75:157–61.

    Article  CAS  PubMed  Google Scholar 

  60. Soma KK, Tramontin AD, Featherstone J, Brenowitz EA. Estrogen contributes to seasonal plasticity of the adult avian song control system. J Neurobiol. 2004;58:413–22.

    Article  CAS  PubMed  Google Scholar 

  61. Weiss G. Relaxin in the male. Biol Reprod. 1989;40:197–200.

    Article  CAS  PubMed  Google Scholar 

  62. Samuel CS, Tian H, Zhao L, Amento EP. Relaxin is a key mediator of prostate growth and male reproductive tract development. Lab Investigat. 2003;83:1055–67.

    Article  CAS  Google Scholar 

  63. de Rienzo G, Aniello F, Branno M, Minucci S. Isolation and Characterization of a Novel Member of the Relaxin/Insulin Family from the Testis of the Frog Rana esculenta 1. Endocrinology. 2001;142:3231–8.

    PubMed  Google Scholar 

  64. De Rienzo G, Aniello F, Branno M, Izzo G, Minucci S. The expression level of frog relaxin mRNA (fRLX), in the testis of Rana esculenta, is influenced by testosterone. J Exp Biol. 2006;209:3806–11.

    Article  PubMed  Google Scholar 

  65. Bottoni L, Lucini V, Massa R. Effect of progesterone on the sexual behavior of the male Japanese quail. Gen Comp Endocr. 1985;57:345–51.

    Article  CAS  PubMed  Google Scholar 

  66. Schneider JS, Burgess C, Sleiter NC, DonCarlos LL, Lydon JP, O'Malley B, et al. Enhanced sexual behaviors and androgen receptor immunoreactivity in the male progesterone receptor knockout mouse. Endocrinology. 2005;146:4340–8.

    Article  CAS  PubMed  Google Scholar 

  67. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.

    Article  CAS  PubMed  Google Scholar 

  69. Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25:3389–402.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.

    Article  CAS  PubMed  Google Scholar 

  72. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21:3787–93.

    Article  CAS  PubMed  Google Scholar 

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

    Google Scholar 

  74. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:1.

    Article  Google Scholar 

  76. Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Nat Acad Sci. 2003;100:9440–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:1.

    Article  Google Scholar 

  78. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2− ΔΔCT method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We would like to thank Chenliang Li and Chunlin Liao for their assistance in the collection of the samples used in this study. We would also like to thank the molecular and behavioral ecology research group for their kind help and useful advice for data analysis and manuscript writing.

Funding

This study has been supported by the National Natural Science Foundation of China (No.31270425 and No. 31470442).

Availability of data and materials

The dataset supporting the conclusions of this article is included within the article and its Additional files 1 and 2. The assembling sequence and expression profiling data have been deposited in NCBI’s Gene Expression Omnibus (GEO) repository through GEO accession number GSE89016 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE89016).

Authors’ contributions

HW conceived the experiments, participated in its design and corrected the manuscript. WZ performed the experiments and wrote the manuscript. WZ, YG and LH conducted the data analysis. YG, JL and EK collected the samples and corrected the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

The sampling and carrying procedures of living toads used in this study had no significant effects on their body condition and survival. The sampling of living toads was permitted by the management bureau of the Badagongshan National Nature Reserve. These procedures and the protocol used to obtain tissue samples from these animals have been previously approved by the Central China Normal University Animal Welfare Committee.

Author information

Authors and Affiliations

  1. Institute of Evolution and Ecology, International Research Centre of Ecology and Environment, School of Life Sciences, Central China Normal University, 152 Luoyulu, Hongshan District, Wuhan, 430079, China

    Wei Zhang, Yue Guo, Jun Li, Li Huang, Eric Gilbert Kazitsa & Hua Wu

Authors
  1. Wei Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yue Guo
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jun Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Li Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Eric Gilbert Kazitsa
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Hua Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Hua Wu.

Additional files

Additional file 1: Table S1.

The statistics and quality of reads throughout all samples in the Leptobrachium boringii transcriptome. Table S2. Significant enrichment of GO terms in the three tissues. Table S3. Differentially expressed genes in the upper jaw skin. Table S4. Differentially expressed genes in the testis. Table S5. Differentially expressed genes in the brain. Table S6. Top 10 most differentially expressed genes in testis and brain of Leptobrachium boringii. Table S7. Specific primers used for quantitative PCR validation. (XLSX 769 kb)

Additional file 2: Figure S1.

Heatmap comparing differentially expressed genes of the three tissues between different breeding periods. (DOCX 265 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, W., Guo, Y., Li, J. et al. Transcriptome analysis reveals the genetic basis underlying the seasonal development of keratinized nuptial spines in Leptobrachium boringii . BMC Genomics 17, 978 (2016). https://doi.org/10.1186/s12864-016-3295-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-016-3295-9

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us