Abstract

Background

   Nuptial pads, a typical sexually dimorphic trait in anurans, are
   located on the first digit of the male forelimb in Rana chensinensis
   and exhibit morphological changes synchronized with breeding cycles.
   However, the genetic mechanisms underlying its formation and seasonal
   changes remain poorly understood.

Results

   To identify genes and biological processes associated with the
   development and seasonal variations of nuptial pads, we conducted a
   comprehensive transcriptome analysis on nuptial pads and hind toe skin
   across both sexes at different breeding periods in R. chensinensis. We
   identified numerous sexually and seasonally differential expression
   genes in nuptial pads. Notably, genes including KRT, TRY, HPDB, AKR1C1,
   and AKR1C3 were identified as potential key regulators of
   keratinization and coloration variation in nuptial pads. We further
   examined gene co-expression modules closely linked to nuptial pad
   development. These modules contained genes involved in signal
   transduction, substance transport, cytoskeletal structure, energy
   metabolism, and protein modification, suggesting that the development
   of nuptial pads is a complex multifaceted regulatory process.
   Furthermore, genes in modules associated with pad development during
   the breeding season were primarily involved in apoptosis, steroid
   hormone synthesis, autophagy, and cytochrome P450 pathways, suggesting
   their pivotal role in pad formation. Additionally, key regulators of
   the cell cycle, such as FOXO4, PIK3C2A, and GSPT2, were implicated in
   influencing nuptial pad development by modulating cell differentiation
   and proliferation.

Conclusions

   Our study provides a valuable reference for investigating the molecular
   basis of sexual dimorphism in R. chensinensis and other amphibian
   species more broadly.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12864-024-11176-3.

   Keywords: Nuptial pads, Rana chensinensis, RNA-seq, Keratin, Steroid
   hormones

Introduction

   The evolution and diversification of sexual dimorphism are of enduring
   interest to evolutionary biologists. Amphibians, which represent a
   pivotal transition from aquatic to terrestrial life among vertebrates,
   display a wide range of sexually dimorphic characteristics.
   Investigating these traits in amphibians is vital for elucidating the
   origins and evolutionary pathways of sexual dimorphism in vertebrates.
   Initial studies in amphibians have focused on sexually dimorphic
   features such as tail fin architecture, pigmentation, vocal sacs, and
   nuptial pads [[30]1]. Of these, nuptial pads have emerged as a focal
   point of research within the anuran clade.

   Nuptial pads, characterized by their extreme morphological diversity
   across various species, serve as a crucial taxonomic criterion for
   amphibians. These unique sexually dimorphic traits play an essential
   role in reproductive behavior, including aiding male competition,
   facilitating grip during amplexus, and releasing protein-based
   pheromones to stimulate reproductive activities [[31]2–[32]7]. Given
   these potential functions, nuptial pads present an ideal model for
   studying sexual dimorphism and the adaptive evolution of reproductive
   behaviors in amphibians. At present, nuptial pads are one of the most
   extensively researched sexually dimorphic traits in anurans [[33]8].

   Nuptial pads consist of thickened layers of epidermal and/or dermal
   tissue layers. They are typically found on the inner margins of the
   thumbs in most male anurans [[34]9–[35]11]. In some species, these
   structures are permanent, while in others, they emerge solely during
   the mating season and recede afterward [[36]8, [37]12]. The seasonal
   transformation of this structure in some species is primarily regulated
   by the synergistic action of androgens and estradiol [[38]13, [39]14].
   However, hormonal regulation is not the sole determinant for nuptial
   pad development; environmental factors such as temperature and
   photoperiod also play a significant role in their seasonal variations
   [[40]15]. Despite this, studies exploring the molecular mechanisms
   underlying the development of amphibian nuptial pads remain scarce. A
   few studies have explored the molecular basis of the nuptial pad glands
   and their secretions [[41]7, [42]16]. These secreted proteins are
   closely associated with the function of nuptial pads in pheromone
   release for reproductive facilitation. However, our understanding of
   the structural and functional characteristics of these proteins, as
   well as the genetic mechanisms underpinning sexually dimorphic gland
   development is still in its early stages.

   The Chinese brown frog (Rana chensinensis) exhibits sexually dimorphic
   traits, particularly the nuptial pads on the male’s forelimbs
   (Fig. [43]1). Only adult males possess these structures, which are
   absent in females. Furthermore, the morphology of male nuptial pads
   exhibits seasonal variations, being more pronounced during the breeding
   season and gradually diminishing thereafter. Located on the first
   digit, these pads feature four tumorous outgrowths, with two proximal
   ones being larger and more pronounced. They bear nuptial spines that
   change color with the seasons. Microscopic examination has shown that
   the pads consist of an epidermal and dermal layer. The epidermis
   comprises a basal layer of densely arranged cuboidal or short columnar
   cells, thin with only two to three cell layers, each with prominent
   nuclei. Granular cells with round to oval nuclei sit above the basal
   layer, enlarging outward and containing keratin granules of varying
   sizes. The dermal layer contains melanocytes and acinar structures of
   the nuptial gland. Seasonal variations affect the epidermal thickness,
   the salience of nuptial spines, and the size and number of nuptial
   gland acini in the nuptial pads of R. chensinensis [[44]17]. The
   development and regression of these nuptial pads are also modulated by
   sex hormones [[45]18], however, the precise molecular mechanisms remain
   to be determined.

Fig. 1.

   [46]Fig. 1
   [47]Open in a new tab

   Different phenotypes of the nuptial pads in male and female R.
   chensinensis during the breeding season. The red circles are nuptial
   pads of the male (left), while there is no nuptial pad on the female’s
   fingers (right)

   To elucidate the molecular mechanisms regulating the development and
   seasonal changes of nuptial pads in R. chensinensis, we performed
   transcriptomic sequencing on nuptial pads and hind toe skin from both
   male and female individuals at the breeding and post-breeding seasons.
   Using differential expression and weighted gene co-expression network
   analysis (WGCNA) analyses, we identified a series of candidate genes
   and biological processes in shaping nuptial pads. The results provide a
   foundation for advancing our understanding of the molecular evolution
   of sexual dimorphism in anurans.

Methods

Tissue collection and RNA isolation

   Adult males (M) and females (F) of R. chensinensis were randomly
   sampled from a breeding site in the Taihang Mountain National Nature
   Reserve, at Yugong Forest Farm, Jiyuan City, Henan Province. The
   collections were conducted during the breeding season in March (stage
   A) and the post-breeding season in May (stage B) of 2021. Given that
   the skin of the hind toe shares a similar epidermal structure with the
   finger skin and is easier to obtain in sufficient quantities, it was
   chosen as the control tissue for the nuptial pad. This selection
   facilitates the identification of genes and pathways associated with
   the development and function of nuptial pads. In accordance with the
   American Veterinary Medical Association (AVMA) Guidelines for the
   Euthanasia of Animals (2020 Edition), our experimental subjects were
   first deeply anesthetized by immersion in a 5 g/L tricaine
   methanesulfonate (MS 222) solution, and subsequently euthanized using
   the method of cerebral spinal destruction [[48]19]. Nuptial pads
   (tissue 1) and hind toe skin (tissue 2) were excised from the subjects,
   immediately placed into sterile tubes, flash-frozen in liquid nitrogen,
   and stored at −80 °C for subsequent analysis. A total of eight groups
   of tissue samples were obtained and labeled as AF1, AF2, AM1, AM2, BF1,
   BF2, BM1, and BM2. Each group contains three biological replicates from
   different individuals, resulting in a total of 24 tissue samples.

   An appropriate 100 mg of each sample was used for RNA extraction by the
   phenol/chloroform method. The RNA degradation and contamination were
   checked via 1% agarose gel electrophoresis. The concentration and
   integrity of the RNA were assessed and confirmed using an RNA Assay
   Kit, a Qubit^® 2.0 Fluorometer, and an RNA Nano 6000 Assay Kit on a
   Bioanalyzer 2100 system, respectively. Only RNA samples with good
   quality were used in further experiments.

Illumina sequencing and gene expression quantification

   The sequencing libraries were prepared using the NEBNext^® UltraTM RNA
   Library Prep Kit from Illumina^® (NEB, USA) and sequenced on the
   Illumina HiSeq 4000 platform. Clean reads were obtained by removing the
   adapter sequences, reads containing poly-N and low-quality reads (reads
   with greater than 5% unknown nucleotides or reads with less than 50 bp)
   from raw reads. Subsequently, the Q20, Q30, and GC content of the clean
   data were calculated. The clean reads were aligned to the R.
   chensinensis reference genome (our unpublished data) using the software
   HISAT2 [[49]20]. The 4.23 Gb of the R. chensinensis reference genome
   was assembled into chromosome-level based on PacBio long reads and HiC
   data, yielding a contig N50 of 2.34 Mb and a scaffold N50 of 484.14 Mb.
   A total of 22,360 protein-coding genes were annotated based on the
   integration of de novo prediction, homology search, and
   transcripts-based assembly methods. Transcripts were assembled using
   StringTie [[50]21]. New transcripts that were not included in the
   reference genome annotation were further verified using the Coding
   Potential Calculator (CPC). HTseq [[51]22] was used to count the read
   numbers mapped to each gene. All steps were performed using default
   parameters. The FPKM of each gene was calculated based on the length of
   the gene and the reads count mapped to this gene. To evaluate the
   reproducibility across samples, we calculated correlation coefficients
   for three biological replicates. A correlation coefficient close to one
   indicates high consistency among replicates.

Differential expression analysis

   To elucidate the genetic and biological pathways underlying the
   formation and seasonal changes of nuptial pads in R. chensinensis, we
   first conducted a comprehensive comparative analysis. We compared gene
   expression patterns related to nuptial pad development among genders
   (males and females), tissues (nuptial pad skin and hind toe skin), and
   breeding periods (breeding and post-breeding seasons). Differential
   expression analysis between two groups was performed using the DESeq2
   package [[52]23]. P values were subsequently adjusted by the
   Benjamini-Hochberg method to control the false discovery rate. Genes
   with an adjusted P value < 0.05 and |log2 fold change (FC)| > 1 were
   identified as differentially expressed genes (DEGs). Gene Ontology (GO)
   enrichment analysis of these DEGs was performed utilizing the
   clusterProfiler R package [[53]24].

Weighted gene co‑expression network analysis (WGCNA)

   To explore the intricate relationships among numerous DEGs in the
   samples and their relationship with the seasonal development of nuptial
   pads, we constructed gene co-expression networks using the WGCNA R
   package [[54]25]. The mean absolute deviation (MAD) is the average
   distance between each data value and the median, which can be used to
   measure the variability of the dataset. To improve the accuracy of
   network analysis, we selected 75% most variable genes to construct the
   weighted gene co-expression networks. We calculated expression
   correlation coefficients to determine the optimal soft threshold and
   construct the gene network using an unscaled topological model.
   Subsequently, a dynamic tree-cutting method was used to identify gene
   modules with similar expression patterns, with a minimum module size of
   50 and a merge height cut-off of 0.25. Considering differences in
   genders (M: males; F: females), tissues (1: nuptial pads; 2: hind toe
   skin), and periods (A: breeding season; B: post-breeding season), we
   categorized the 24 samples into eight phenotypic traits: AF1, AF2, AM1,
   AM2, BF1, BF2, BM1, and BM2. A module was considered significantly
   correlated with the trait when P < 0.05 and the correlation coefficient
   (r) > 0.5. Furthermore, we conducted GO and KEGG enrichment analysis
   for genes in modules significantly related to the formation of nuptial
   pads using the ClusterProfiler R package. Finally, the 50 most highly
   connected genes in key modules were used to generate a correlation
   network in Cytoscape. The top 10 genes ranked by their degree of
   connectivity were defined as the hub genes of the module.

Results

Transcriptomic data constitution and sample clustering

   During the breeding and post-breeding seasons, we performed
   transcriptomic sequencing on eight samples of nuptial pads and hind toe
   skin from both male and female individuals. Transcriptomic sequencing
   produced approximately 1,185 M raw reads and 1,162.88 M high-quality
   clean reads. The minimum number of raw reads per sample was 40.9 M, and
   the minimum number of clean reads was 40.03 M. The quality of clean
   reads was notably high, with at least 96.23% at Q20 and 90.52% at Q30
   (Additional file 1: Table [55]S1). We further evaluated the
   reproducibility of biological replicates based on gene expression
   levels. Pearson correlation coefficients for the three replicates
   ranged from 93.42 to 96.64% (Fig. [56]2), reflecting the degree of
   linear association between variable pairs. A higher coefficient value
   suggests a stronger correlation.

Fig. 2.

   [57]Fig. 2
   [58]Open in a new tab

   Pearson correlation coefficients between samples. The x and y axes
   correspond to the 24 samples, and the varying colors denote the
   respective correlation coefficient

Sexually different expression genes in nuptial pads at different periods

   Nuptial pads are sexually dimorphic traits that change seasonally, so
   we first compared the gene expression differences between sexes, which
   helps to reveal the molecular basis of sexual dimorphism. We initially
   examined the DEGs between males and females in terms of different
   tissues (1: nuptial pads and 2: hind toe skin) and different periods
   (A: breeding season and B: post-breeding season). The analysis revealed
   that the quantity of DEGs between males and females was consistently
   higher in nuptial pads compared to hind toe skin in two examined
   periods. Furthermore, a significant fluctuation in the number of
   nuptial pads’ DEGs was observed between these periods (Table [59]1,
   Additional file 1: Table [60]S2−5). Specifically, during the breeding
   season (stage A), nuptial pads exhibited 1,956 DEGs, of which 1,053
   were highly expressed in males (1053/1956). The number of sexually DEGs
   decreased markedly at the post-breeding season (stage B), with only 75
   DEGs identified in nuptial pads, including 47 male-biased expression
   genes (47/75).

Table 1.

   Statistics of sex-specific differentially expressed genes
   Compared tissues Compared groups Up   Down Total DEGs
   Nuptial pad      AF1 vs. AM1     1053 903     1956
                    BF1 vs. BM1      47   28      75
   Hind toe skin    AF2 vs. AM2     148   78     226
                    BF2 vs. BM2      27   15      42
   [61]Open in a new tab

   Notably, 1796 sexually differentially expressed genes present only in
   the nuptial pad during the breeding season, and 42 such genes are
   present after the breeding season (Fig. [62]3). We speculate that these
   sexually differentially expressed genes, which are exclusive to the
   nuptial pad at different times, may be associated with the development
   and regression of the nuptial pad.

Fig. 3.

   Fig. 3
   [63]Open in a new tab

   Venn diagram of sexually different expression genes in the nuptial pad
   and skin during different periods

   Subsequently, we performed a GO enrichment analysis on these genes
   (Fig. [64]4, Additional file 1: Table S11-12). Within the nuptial pads,
   sexually DEGs at the breeding season (from AF1 vs. AM1) were
   significantly enriched in biological processes such as protein
   synthesis and modification (translation GO:0006412, structural molecule
   activity GO:0005198), as well as metabolic processes (GO:0008152),
   redox processes (GO:0055114), microtubule-based processes (GO:0007017),
   and positive regulation of apoptotic process (GO:0043065). At the
   post-breeding season, the sexually DEGs (from BF1 vs. BM1)
   significantly enriched in functions including helicase activity
   (GO:0004386), metallopeptidase activity (GO:0008237), and
   metalloendopeptidase activity (GO:0004222) (Fig. [65]4).

Fig. 4.

   [66]Fig. 4
   [67]Open in a new tab

   Significantly enriched GO terms for sex-specific DEGs (q < 0.05). The
   x-axis represents different comparison groups, and the y-axis
   represents GO terms. The circle size represents the number of genes
   included in the item. The circle colors represent adjusted p values (q
   values) from enrichment analyses. BP, Biological Process; CC, Cellular
   Component; MF, Molecular Function

Comparative gene expression analysis between nuptial pads and Hind toe skin

   To gain a deeper understanding of the biological processes that govern
   nuptial pad development, we then further compared the gene expression
   differences between different tissues to explore the essential genetic
   expression differences between nuptial pads and hind toe skin. This
   comparison allowed us to identify DEGs by contrasting the gene
   expression profiles of two skin tissues from the same sex and
   developmental stages. During the breeding season, we identified 1,067
   DEGs between the hind toe skin and nuptial pads of males, and 48 DEGs
   in the female counterpart comparison. During the post-breeding season,
   the number of DEGs decreased to 119 for males’ comparison and 64 for
   females (Table [68]2, Additional file 1: Table S6-9). Notably, 983
   genes were only expressed in the group of AM2 vs. AM1, while 49 genes
   were unique to the group of BM2 vs. BM1 (Fig. [69]5).

Table 2.

   Comparative statistics of differentially expressed genes in nuptial
   pads and Hind toe skin
   Compared periods Compared group Up  Down Total DEGs
   Breeding period  AF2 vs. AF1    32   16      48
                    AM2 vs. AM1    764 303     1067
   After breeding   BF2 vs. BF1    36   28      64
                    BM2 vs. BM1    60   59     119
   [70]Open in a new tab

Fig. 5.

   Fig. 5
   [71]Open in a new tab

   Venn diagram of differentially expressed genes between nuptial pads and
   hind toe skin at different periods

   We hypothesized that genes only differentially expressed between males’
   nuptial pads and hind toe skins at the breeding season (AM2 vs. AM1)
   may contribute to nuptial pad development. Conversely, DEGs unique to
   the BM2 vs. BM1 comparison could be implicated in the post-breeding
   regression of nuptial pads. According to the GO enrichment analyses,
   DEGs exclusive to the AM2 vs. AM1 group were notably overrepresented in
   biological processes including intermediate filament organization
   (GO:0005882), carbohydrate metabolism processes (GO:0005975),
   transmembrane transport (GO:0055085), protein dephosphorylation process
   (dephosphorylation GO:0016311, protein dephosphorylation GO:0006470),
   and the positive regulation of apoptosis process (GO:0043065)
   (Fig. [72]6 and Additional files: Table S13-14). We identified that
   genes such as keratin 75 (KRT75), dual specificity phosphatase (DUPD1),
   protein phosphatase Mg2+/Mn2+-dependent 1 H (PPM1H), BCL2 interacting
   protein 3 (BNIP3), actinin alpha 3 (ACTN3),
   beta-1,4-galactosyltransferase 4 (B4GALT4), sodium voltage-gated
   channel alpha subunit 2 (SCN2A), ryanodine receptor 1 (RYR1), aquaporin
   2 (AQP2), solute carrier family 15 member 1 (SLC15A1), and potassium
   voltage-gated channel subfamily A member 3 (KCNA3) exhibited
   significant upregulation in males’ nuptial pads compared to the hind
   toe skin. Furthermore, the DEGs exclusive to the comparison of BM2 vs.
   BM1 were significantly enriched in the category of methyltransferase
   activity (GO:0008168) at the molecular function level (q < 0.05).

Fig. 6.

   [73]Fig. 6
   [74]Open in a new tab

   GO significantly enrichment results of DEGs expressed in the comparison
   group between nuptial pad and hind toe skin (q < 0.05). Red, light
   blue, and green represent biological process (BP), cellular component
   (CC), and molecular function (MF), respectively. The first lap
   indicates the GO terms. The second lap indicates the number of genes in
   the genome background and q-values for gene enrichment for the
   specified GO terms. The third lap indicates a bar graph of enriched
   up-regulated genes. The fourth lap indicates a bar graph of enriched
   down-regulated genes. The fifth lap indicates the enrichment factor for
   each GO term

Gene expression differential analysis of nuptial pads across breeding cycles

   To gain a better understanding of the biological mechanisms governing
   seasonal variations in nuptial pad morphology, we further focused on
   the gene expression differences in the nuptial pads between different
   periods. In males’ nuptial pads, a total of 3,426 DEGs were identified
   between the breeding and post-breeding periods, comprising 1,864
   upregulated and 1,562 downregulated genes (Fig. [75]7A, Additional file
   1: Table S9-10). We proposed that the 2,179 DEGs uniquely expressed in
   male nuptial pads were associated with the seasonal transformation
   observed in these structures (Fig. [76]7B). These DEGs were
   significantly enriched in the intermediate filaments and the protein
   tyrosine phosphatase activity (q < 0.05).

Fig. 7.

   [77]Fig. 7
   [78]Open in a new tab

   DEGs exclusively expressed during different breeding periods in pad and
   skin and their GO enrichment results. A: Differentially expressed genes
   in nuptial pads and hind toe skins between breeding and post-breeding
   periods. B: Venn diagram of differentially expressed genes in male
   nuptial pads and skin across different breeding periods. C:
   Significantly enriched GO terms for DEGs in the nuptial pad during
   different breeding periods (BM1 vs. AM1, q < 0.05)

   We separately examined upregulated and downregulated DEGs in nuptial
   pads between breeding and post-breeding season. The upregulated genes
   in the nuptial pads during the breeding season were significantly
   enriched in protein tyrosine phosphatase activity (GO:0004725) and
   transmembrane transporter activity (GO:000520). Conversely, the
   downregulated genes, which showed higher expression post-breeding in
   the nuptial pads, were significantly enriched in processes such as
   protein modification via gamma-glutamyltransferase activity
   (GO:0003810), peptide cross-linking (GO:0018149), protein folding with
   peptidyl-prolyl isomerase activity (GO:0000413, GO:0003755), and
   maintenance of the extracellular matrix structure (GO:0005201;
   Fig. [79]7C).

Seasonal development of nuptial pads in R. chensinensis: Biological processes
and genes associations

   To elucidate the primary biological processes and gene associations
   implicated in the seasonal development of nuptial pads in R.
   chensinensis, we performed an extensive comparative analysis of DEGs.
   This analysis encompassed contrasts between sexes (AF1 vs. AM1 & BF1
   vs. BM1), between tissues (AM2 vs. AM1 & BM2 vs. BM1), and between
   periods (BM1 vs. AM1).

   We initially focus on the results related to the development of male
   nuptial pads during the breeding season for the groups AF1 vs. AM1, AM2
   vs. AM1, and BM1 vs. AM1 (upregulated genes). We observed that among
   the DEGs in the three comparative groups, 192 genes are consistently
   highly expressed in AM1. These genes are primarily involved in cellular
   metabolism (ACY3, AGPAT4, ARSB, etc.), transmembrane transport
   (transmembrane protein gene families, etc.), and signal transduction
   (ECM1, GRB14, RRAS2, etc.), which are crucial biological functions.
   Specifically, we identified genes associated with sex hormone synthesis
   and signaling, apoptosis, keratinization, and melanin synthesis,
   including cytochrome P450 family 2 subfamily g polypeptide 1 (CYP2G1),
   cytochrome P450 family 46 subfamily A member 1 (CYP46A1),
   sulfotransferase family 2B member 1 (SULT2B1), DAB2 interacting protein
   (DAB2IP), forkhead box A1 (FOXA1), insulin like growth factor binding
   protein 2 (IGFBP2), keratin 17 (KRT17), and tyrosinase (TYR).

   Furthermore, by comparing the results of sex differences (AF1 vs. AM1)
   and tissue differences (AM2 vs. AM1), we found that a series of
   biological processes involving cellular structure and redox activities
   were significantly enriched, including activity related to structural
   molecule (GO:0005198), intermediate filament (GO:0005882),
   oxidoreductase activity (GO:0016491), redox process (GO:0055114), and
   oxidoreductase activity with a single donor that incorporates molecular
   oxygen (GO:0016791). Specifically, five keratin genes, two genes
   associated with tight junction proteins, two genes from the aldo-keto
   reductase family, four from the cytochrome P450 family, and seven
   related to tyrosine metabolism demonstrated substantial overexpression
   in the nuptial pads of males (Additional file 1: Table S17). However,
   no significant enrichment of biological processes was observed in the
   period difference (BM1 vs. AM1 upregulated genes). We speculate that
   these biological processes may be related to the unique development of
   the nuptial pad, a special dermal derivative, compared to other skin
   tissues.

   Subsequently, we performed an integrated comparative analysis of the
   results from the three groups associated with the post-breeding
   regression of male nuptial pads: BF1 vs. BM1, BM2 vs. BM1, and BM1 vs.
   AM1 (downregulated genes). By comparing the DEGs across the three
   groups, the significantly downregulation of glutathione S-transferase
   omega 1 (GSTO1) in the nuptial pad may be associated with its
   degeneration. However, we observed no significant enrichment of DEGs
   when comparing between sexes (BF1 vs. BM1), between tissues (BM2 vs.
   BM1), or across different time points (BM1 vs. AM1 downregulated genes)
   post-breeding season.

Gene co-expression modules related to nuptial pad development in R.
chensinensis

   To further elucidate the regulatory mechanisms underlying the formation
   of nuptial pads in R. chensinensis, we employed WGCNA analysis on 24
   samples. We identified 18,391 with the top 75% MAD values for
   subsequent WGCNA analysis. Based on the correlation coefficients among
   the selected genes, we constructed a gene clustering dendrogram using a
   soft-thresholding power of 22. The gene modules were further
   categorized with a similarity threshold of 0.8 and a minimum module
   size of 50. Finally, we delineated 28 modules with sizes ranging from
   72 to 6,288 genes, among which a “gray” module encompassed genes
   unassignable to any specific modules (Fig. [80]8A-B, Additional file 1:
   Table S19). Gene modules and phenotypic trait correlation analysis
   identified eight modules associated with the development of male
   nuptial pads during the breeding season (AM1; Fig. [81]8C). Among
   these, five modules (darkorange2, tan, lightsteelblue1, lightcyan, and
   midnightblue) exhibited a significant positive correlation with the
   presence of male nuptial pads at the breeding season (r > 0.5,
   p < 0.05). Conversely, three modules (orange, darkgreen, and violet)
   were found to be significantly negatively correlated (r > 0.5,
   p < 0.05). Additionally, the palevioletred3 module demonstrated a
   significant positive correlation with the regression of male nuptial
   pads during the post-breeding season (BM1) (r > 0.5, p < 0.05).

Fig. 8.

   [82]Fig. 8
   [83]Open in a new tab

   WGCNA results of genes. A: The gene cluster dendrogram constructed by
   gene correlation coefficients (r). B: Soft threshold used in WGCNA. C:
   Relationships between modules and traits

Functional enrichment of genes in target modules

   GO enrichment analysis revealed distinct biological processes
   associated with various gene modules. The darkorange2 module (r = 0.60,
   p < 0.01) was predominantly enriched in the cell surface receptor
   signaling pathway (GO:0007166). In contrast, the lightcyan module
   (r = 0.60, p < 0.01) showed significant enrichment in complement
   activation (GO:0006956) and anion transmembrane transporter activity
   (GO:0008509). The lightsteelblue1 module (r = 0.72, p < 0.01) primarily
   involved genes associated with cytoskeleton structure and function,
   particularly those related to the structural constituent of the
   cytoskeleton (GO:0005200), microtubule-based process (GO:0007017), and
   microtubule (GO:0005874). The midnightblue module (r = 0.61, p < 0.01)
   was enriched in genes involved in glycolytic process (GO:0006096),
   acetylcholine-gated cation-selective channel activity (GO:0022848),
   motor activity (GO:0003774), and the myosin complex (GO:0016459).
   Meanwhile, the tan module (r = 0.61, p < 0.01) was enriched in protein
   modification processes, including protein ubiquitination (GO:0016567),
   protein deubiquitination (GO:0016579), phosphatidylinositol
   dephosphorylation (GO:0046856), and autophagy (GO:0006914). The
   darkgreen module (r=−0.70, p < 0.01) mainly involved genes in energy
   metabolism, such as ATP synthesis coupled proton transport
   (GO:0015986), ATP hydrolysis coupled proton transport (GO:0015991),
   cytochrome-c oxidase activity (GO:0004129), and translation
   (GO:0006412). Moreover, the orange module (r=−0.58, p < 0.01) was
   associated with methyltransferase activity (GO:0008168), and flavin
   adenine dinucleotide binding (GO:0050660). Lastly, the violet module
   (r=−0.66, p < 0.01) was mainly involved in rRNA processing
   (GO:0006364). These gene modules are implicated in biological processes
   relevant to the development of male nuptial pads during the breeding
   season, as depicted in Additional file 2: Fig. [84]S1. The genes in the
   palevioletred3 module (r = 0.62, p < 0.01) were significantly enriched
   in the intermediate-filament GO:0005882 category. These genes are
   associated with the regression of male nuptial pads following
   post-breeding season (Additional file 2: Fig. [85]S1).

   KEGG enrichment analysis indicated that the genes within the lightcyan
   module were mainly implicated in apoptosis (ko04215), cytochrome P450
   activities (ko00199), and steroid hormone biosynthesis (ko00140).
   Conversely, genes in the tan module were mainly associated with
   apoptosis (fly ko04214) and autophagy (other ko04136). The violet
   module’s genes were primarily engaged in intracellular protein
   synthesis and mitochondrial biogenesis, incorporating ribosome
   biogenesis (ko03009), mitochondrial biogenesis (ko03029), messenger RNA
   biogenesis (ko03019), and translation factor (ko03012). Genes from the
   darkgreen module showed a strong association with ribosomes (ko03010,
   ko03011), oxidative phosphorylation (ko00190), protein export
   (ko03060), and propanoate metabolism (ko00640) (Fig. [86]9). These
   pathways have a profound correlation with the development of male
   nuptial pads during the breeding season (Additional file 2: Fig.
   [87]S2−5).

Fig. 9.

   [88]Fig. 9
   [89]Open in a new tab

   KEGG enrichment of key modules. The figure displays only the terms with
   q < 0.05 for each module. The x-axis represents each module, and the
   y-axis represents KEGG terms

   During the post-breeding season, genes in the palevioletred3 module
   related to male nuptial pads were mainly implicated in the
   Staphylococcus aureus infection (ko05150). Amphibian skin secretes a
   variety of antimicrobial peptides which likely contribute to resistance
   against S. aureus infection. Moreover, the native habitats of
   amphibians might provide favorable conditions for the proliferation of
   S. aureus, thus elevating the infection risk.

Identification of hub genes and networks

   The top 50 genes, selected for their high connectivity within the focal
   module, were used to construct a co-expression network via Cytoscape.
   The ten genes exhibiting the highest degree of connectivity were
   designated as hub genes of the module. The gene interaction networks in
   each key module and the hub genes are displayed in Fig. [90]10 and
   Additional file 1: Table S20.

Fig. 10.

   [91]Fig. 10
   [92]Open in a new tab

   Gene interaction network diagram of the top 50 genes ranked by kME
   values in each key module. The red genes represent the hub genes of
   each module

Discussion

   This study utilized differential transcriptome expression analysis and
   WGCNA methods to investigate the molecular mechanisms of the formation
   and seasonal development of nuptial pads in R. chensinensis. We
   pinpointed a number of pivotal pathways and genes implicated in nuptial
   pad formation and regression. Although the lack of experimental samples
   hinders further validation of these genes, this study does enhance our
   understanding of the molecular determinants of the nuptial pads in
   amphibians.

   Despite the near-identical gene sequences shared by both sexes within a
   species, distinctive sexual dimorphisms persist [[93]26]. Generally,
   the emergence of such characteristics stems from differential gene
   expression between the sexes at various developmental stages [[94]27].
   Our study revealed a marked discrepancy in the quantity of sex-specific
   DEGs present in the nuptial pads versus the hind toe skin of R.
   chensinensis, both of which are skin tissues. Specifically, during the
   breeding season, there is a nearly ninefold discrepancy in the number
   of sex-differentially expressed genes between the nuptial pad and the
   skin of the hind toe, with counts of 1956 and 226 genes, respectively.
   This suggests a substantial set of genes uniquely expressed in nuptial
   pads that are essential for their structural development, unlike those
   in typical skin tissues.

   In analyzing differential gene expression between the nuptial pads and
   hind toe skin of R. chensinensis, we identified a suite of genes
   prominently expressed during nuptial pad development. These genes are
   primarily involved in cellular structure maintenance and redox
   biological processes. Specifically, throughout the breeding season, the
   activities related to structural molecule and intermediate filaments
   were significantly upregulated in nuptial pads. Keratin genes such as
   KRT75, KRT13, KRT5, KRT14, KRT12.4.L, and KRT17 emerged as key genes.
   Particularly, KRT17 showed significant and specific high expression
   across all comparison groups during the breeding period. Keratins, a
   type of intermediate filament proteins, serve as the primary structural
   components of epithelial cells, the stratum corneum of the skin, and
   structural proteins of ectoderm cells [[95]28]. Their link to sexual
   dimorphism in amphibian skin structures is well-documented [[96]29,
   [97]30]. Nuptial pads, as specialized keratinized skin tissue, exhibit
   a close association with keratin genes. During the breeding season,
   high expression of keratin genes in nuptial pads may lead to
   intensified keratin synthesis and skin keratinization in R.
   chensinensis. After the breeding season, the expression of keratin
   genes wanes, culminating in the gradual recession of the keratinized
   pads. Therefore, the dynamic expression of keratin genes is
   instrumental in the seasonal development of nuptial pads in R.
   chensinensis.

   Beyond keratins, genes implicated in melanin synthesis and metabolism
   were also identified in nuptial pads. TYR, the tyrosinase gene, is
   pivotal in vertebrates’ melanin synthesis, catalyzing the conversion of
   tyrosine to dopamine (DOPA) and subsequently to dopamine quinone (DOPA
   quinone), which ultimately leads to melanin formation [[98]31].
   Furthermore, we detected the gene encoding 4-hydroxyphenylpyruvate
   dioxygenase b (HPDB). HPDB is a component of the tyrosine metabolism
   pathway, transforming 4-hydroxyphenylpyruvate acid (HPPA) into
   homogentisic acid (HGA), a melanin synthesis precursor [[99]32]. The
   nuptial pads of R. chensinensis exhibit a darker hue during the
   breeding season, which subsequently fades during the post-breeding
   season. Both TYR and HPDB show significant upregulation in nuptial pads
   during the breeding season. We propose that the color variations of
   nuptial pads across different periods are not only attributable to
   keratinization levels but also importantly influenced by TYR and HPDB,
   which play crucial roles in the color development and maintenance of
   nuptial pads.

   Our research has shown that within aldo-keto reductase family, the
   AKR1C1 and AKR1C3 genes are specifically and highly expressed in
   nuptial pads. These genes encode key enzymes that catalyze the
   reduction of aldehydes and ketones, playing a crucial role in the
   biosynthesis and metabolism of steroids and bile acids [[100]33]. AKR1
   enzymes are capable of facilitating the reduction reaction of various
   steroid hormones, such as androgens, estrogens, and progestogens,
   typically by converting the ketone group to an alcohol group. This
   alteration modulates the hormones’ activity and bioavailability. For
   example, the transformation of androgens into more bioactive forms,
   like dihydrotestosterone (DHT), augments their effects [[101]34]. Given
   that the development and regression of nuptial pads are
   androgen-dependent [[102]18], we hypothesize that the elevated
   expression of AKR1C1 and AKR1C3 during the breeding season may
   intensify androgen activity, fostering maximal nuptial pad development.
   It is noteworthy that, in our comparative analysis across different
   periods of the nuptial pad, we did not observe significant changes in
   the expression of these genes. This may suggest that these genes are
   highly expressed in the male nuptial pad and, compared to other skin
   tissues, may be associated with the formation of this specialized
   dermal derivative.

   Based on differential gene expression analyses, we identified genes
   that exhibited significant differences in expression across sexes,
   tissues, and seasons. These genes potentially played a direct role in
   the development and seasonal changes of nuptial pads. Given that genes
   often interact with other genes when performing their biological
   functions, they tend to show a trend of co-expression. Consequently, we
   further utilized WGCNA analysis to identify gene modules that were
   highly co-expressed and related to the development of the nuptial pad.
   We believed that by integrating these two complementary sets of
   results, a more comprehensive understanding of the molecular mechanisms
   underlying the development of the nuptial pad could be achieved.

   Through WGCNA analysis, we identified several gene modules associated
   with the seasonal development of nuptial pads. These modules comprise
   genes involved in various complex biological processes, such as signal
   transduction, substance transport, cytoskeleton and structure, energy
   metabolism, and protein modification. This indicates that the
   development of nuptial pads is regulated by a comprehensive network.
   Notably, gene modules linked to nuptial pad development during the
   breeding season were implicated in essential processes like cell
   apoptosis, steroid hormone biosynthesis, cell autophagy, and cytochrome
   P450 pathways. These processes may be central to the formation of
   nuptial pads.

   As a specialized skin derivative, the development of nuptial pads
   requires programmed cell death. The stratum corneum, representing the
   epidermis’s outermost layer, is composed of anuclear and deceased
   keratinocytes that arise from apoptosis [[103]35]. Genes associated
   with apoptosis, including classic apoptotic factors CASP7, CASP8, and
   APAF1, were significantly enriched within the modules correlated with
   male R. chensinensis nuptial pad development during the breeding
   season. APAF1 facilitates the mitochondrial apoptosis pathway,
   contributing to the apoptotic bodies’ formation, activation of Caspase
   9, and the subsequent triggering Caspase 3, thus promoting apoptosis 36
   [[104]36]. Following the breeding season, no apoptotic genes or
   pathways were discerned in the related modules. We propose that the
   persistence of keratinized nuptial pads after the breeding season may
   be resulted from a diminished reception of apoptotic signals by the
   skin’s epidermal cells, leading to reduced keratinocyte formation and
   its gradual loss over time due to natural wear.

   Autophagy, a pervasive cellular mechanism, allows eukaryotic cells to
   degrade and recycle cellular components, addressing damaged or toxic
   substances during development or stress [[105]37]. Typically, a basal
   level of autophagy is maintained to preserve cellular equilibrium.
   Under unfavorable growth conditions, autophagy can be induced to
   eliminate superfluous or damaged materials and to recycle components,
   thereby providing substrates and products necessary for cellular
   metabolism [[106]38]. The core autophagy process is mediated by
   autophagy-related (ATG) genes [[107]39]. The gene module associated
   with the development of R. chensinensis nuptial pads, designated as the
   ‘tan module’, includes essential autophagy genes such as ATG16, ATG101,
   ATG2, ATG4, ATG7, and ATG9. The critical role of autophagy in
   restructuring internal cellular architecture is evident during the
   nuptial pads’ rapid proliferation and differentiation, a process
   necessitated by the complex environments and competitive pressures
   encountered by R. chensinensis during the breeding season.

   In vertebrates, sex hormones play a pivotal role in the induction of
   sexually dimorphic traits. Seasonal fluctuations in hormone
   concentrations and reproductive features are particularly pronounced in
   amphibians that breed cyclically [[108]40]. Our research indicates that
   genes in the lightcyan module, which correlate with the development of
   nuptial pads during the breeding season, are significantly
   overrepresented in the steroid hormone synthesis pathway. This group
   includes CYP11B, CYP21A2, and SULT2B1, all of which are integral to
   steroid hormone production. The enzyme 11β-hydroxylase, encoded by
   CYP11B, is instrumental in androgen synthesis [[109]41]. CYP21A2 is
   crucial for the generation of adrenal cortical hormones, catalyzing the
   transformation of 17OHP into 11-deoxycortisol, which is a precursor of
   aldosterone. These conversions are key to the production of
   glucocorticoids and mineralocorticoids, respectively [[110]42]. SULT2B1
   contributes to the sulfation of hormones, including sex hormones, and
   can sulfate cholesterol, a sex hormone synthesis precursor [[111]43].
   The modulation of SULT2B1 activity may influence signal transduction of
   sex hormones; sulfated hormones might fail to bind to their receptors,
   diminishing their biological activity. Such regulatory mechanisms could
   be crucial in the differentiation and preservation of sexual dimorphism
   in R. chensinensis nuptial pads. For many amphibians, sex hormones
   trigger seasonal morphological shifts, including the emergence of
   secondary sexual features like keratinized spines. The seasonal
   variations in R. chensinensis nuptial pads might be associated with the
   steroid hormone synthesis-related genes highlighted in our study.

   Moreover, we identified several genes involved in cell cycle control
   and proliferation within pivotal hub genes of specific modules,
   including GON4L, MTMR14, ARHGEF2, FOXO4, and GSPT2 [[112]44–[113]48].
   These genes significantly influence the development of sexually
   dimorphic traits by coordinating cellular differentiation,
   intracellular signal, and proliferation, thus playing a vital role in
   the complex regulatory mechanisms governing nuptial pad formation.

   Through comprehensive comparative transcriptome analysis and WGCNA, we
   have identified a suite of biological processes and core genes
   potentially associated with the seasonal development of R. chensinensis
   nuptial pads. Moreover, it is necessary to acknowledge the limitations
   of this work. We only investigated the the nuptial pad during the
   breeding and post-breeding periods, omitting the pre-breeding phase,
   which may affect our understanding of the developmental context of the
   nuptial pad. Furthermore, due to the lack of further experimental
   validation of the analysis results, there are still some genes whose
   functions in the nuptial pad we cannot fully elucidate. Nonetheless,
   our extensive analysis, coupled with insights from existing research,
   provides a significant compilation of biological processes, pathways,
   and candidate genes that enhance our understanding of the molecular
   mechanisms of nuptial pad development in R. chensinensis. This study
   offers a valuable reference for investigating sexual dimorphism in
   amphibian skin specializations. In our future work, we intend to
   investigate samples from a more complete breeding cycle to more
   thoroughly explore the molecular mechanisms of the nuptial pad and
   other similar sexually dimorphic structures.

Conclusions

   Our transcriptome analysis revealed a significant upregulation in the
   expression of keratin genes TRY and HPDB in the nuptial pads of R.
   chensinensis during the breeding season, underscoring their crucial
   roles in pad morphology. Biological processes such as signal
   transduction, cytoskeletal structure, energy metabolism, and protein
   modification are intricately connected to nuptial pad development in
   this stage. KEGG pathway enrichment analysis has disclosed that
   pathways associated with apoptosis, autophagy, and steroid hormone
   synthesis and their respective genes are intricately linked to the
   developmental stages of nuptial pads. A decline in the expression of
   genes associated with sex hormone synthesis from the breeding to the
   post-breeding season suggests that the seasonal development of R.
   chensinensis nuptial pads may be influenced by hormonal fluctuations,
   supporting the findings of previous studies on hormonal control of
   nuptial pads. Our findings demonstrate that molecular biological
   processes, pathways, and key genes are intimately associated with the
   development of the R. chensinensis nuptial pads, providing valuable
   insights for the study of sexual dimorphism in amphibian skin
   derivatives. Future research will necessitate further experimental
   validation of these identified genes and pathways.

Supplementary Information

   [114]Supplementary Material 1.^ (1,010.1KB, xlsx)
   [115]Supplementary Material 2.^ (906.3KB, docx)

Acknowledgements