- Research article
- Open Access
Transcriptomic analysis of differentially expressed genes in the oviduct of Rhacophorus omeimontis provides insights into foam nest construction
BMC Genomicsvolume 20, Article number: 562 (2019)
The production of foam nests is one of the strategies that has evolved to allow some anuran species to protect their eggs and larvae. Despite considerable knowledge of the biochemical components of and construction behavior leading to anuran foam nests, little is known about the molecular basis of foam nest construction. Rhacophorus omeimontis presents an arboreal foam-nesting strategy during the breeding season. To better understand the molecular mechanism of foam nest production, transcriptome sequencing was performed using the oviduct of female R. omeimontis during the period when foam nest production began and the period when foam nest production was finished.
The transcriptomes of six oviduct samples of R. omeimontis were obtained using Illumina sequencing. A total of 84,917 unigenes were obtained, and 433 genes (270 upregulated and 163 downregulated) were differentially expressed between the two periods. These differentially expressed genes (DEGs) were mainly enriched in extracellular space and extracellular region based on Gene Ontology (GO) enrichment analysis and in the pathways of two-component system, cell adhesion molecules, steroid hormone biosynthesis and neuroactive ligand-receptor interaction based on Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. Specifically, genes encoding lectins, surfactant proteins and immunity components were highly expressed when the foam nest construction began, indicating that the constituents of foam nests in R. omeimontis were likely a mixture of surfactant, lectins and immune defense proteins. During the period when foam nest production was finished, genes related to lipid metabolism, steroid hormone and immune defense were highly expressed, indicating their important roles in regulating the process of foam nesting.
Our study provides a rich list of potential genes involved in the production of foam nests in R. omeimontis. These results provide insights into the molecular mechanisms underlying the process of foam nest construction and will facilitate further studies of R. omeimontis.
Foam nest production is a breeding strategy that has evolved in some homopteran bugs, tunicates, fishes and anuran amphibians to protect their eggs and larvae from environmental challenges [1,2,3,4]. In anuran amphibians, foam nest construction is associated with adaptation to terrestrial life, and 10 out of the 41 different reproductive modes found in anurans involve foam nest construction [5,6,7]. Depending on the anuran species, the foam nests are produced in underground burrows, float on the water surface, or are even suspended from vegetation [4, 8,9,10]. Several functions have been identified for the production of anuran foam nests, such as protecting eggs and embryos against predators, desiccation, microbial degradation and thermal damage, as well as supplying oxygen and serving as a food source for tadpoles [2, 8, 11,12,13].
The production of foam nests has evolved independently multiple times in several lineages of anurans; thus, anurans are a particularly interesting case for which to study foam nests [6, 14]. Recent studies on anuran foam nests have mostly focused on their structure and biochemical composition and the construction behavior leading to them [4, 13, 15,16,17]. Well-studied foam-nesting species include frogs from the genus Leptodactylus, especially the túngara frog [4, 10, 12, 14, 15]. Hissa et al.  and Fleming et al.  reported that the foam nest components of Leptodactylus frogs were a mixture of proteins termed ranaspumins, with surfactant activity, cystatin activity and lectin composition. The foam components were mainly secreted from the posterior convolutions of the oviducts in mature female foam-nesting frogs [18, 19].
Despite considerable knowledge of the biochemical components of and construction behavior leading to anuran foam nests, little is known about the molecular basis of and metabolic pathways involved in the process of foam nest production. In addition, published studies on anuran foam nests typically concern the water-based foams produced by members of Leptodactylidae, and whether the arboreal foam nests of anuran species have similar components remains unclear [10, 20].
The Omei treefrog, Rhacophorus omeimontis (Anura: Rhacophoridae), presents a foam-nesting strategy during breeding seasons, and females commonly mate with multiple males in one spawning group [21,22,23]. During each breeding season, males first arrived at a permanent pond surrounded by vegetation. The female is usually grasped by a male at the time of arrival and then moves to find an appropriate spawning site, usually on a leaf above the water surface [22, 23]. The female releases a clutch of eggs together with foam precursor fluid; at the same time, the paired male expels the sperm and, using a rapid “egg-beater” motion of the legs, whips the mixture into a white foamy mass containing the fertilized eggs. Then, other males gradually join the mating pair, forming a spawning group (Fig. 1A). At the end of spawning, the males leave the mating group, and the female remains until leaf nesting is finished [23,24,25]. After the mating period, the foam nests remain on the vegetation or leaves for almost one month (Fig. 1B). The eggs and embryos in the foam nests commonly exhibit a very low mortality [21, 23]. Therefore, R. omeimontis may be an appropriate model for studying the molecular basis of arboreal foam nest production in anuran species.
Transcriptome sequencing is an efficient way to examine the global molecular response of living organisms to the conditions of certain life history periods [26, 27]. Differential gene expression analysis provides a method for studying the molecular basis of context-specific tissue activity without the need for a reference genome [27, 28]. In this study, we performed transcriptome sequencing of the oviduct of female R. omeimontis from the stage when foam nest production was finished (AFNP_O) to the stage before foam nest production (BENP_O). Then, we compared the differential gene expression profiles of the oviducts between these two periods of foam nest production. Our aim was to investigate genes and pathways that may play important roles during the production of leaf foam nests.
De novo transcriptome sequencing of R. omeimontis
To investigate the expression changes of genes in the oviduct of female R. omeimontis during the production of foam nests, we performed transcriptome sequencing of R. omeimontis oviducts during two periods (AFNP_O and BENP_O) with three biological replicates in each period. Due to the absence of a reference genome for R. omeimontis, we de novo assembled a transcriptome as a reference for read mapping and gene expression profiling in this species. Illumina sequencing of the R. omeimontis oviducts generated a total of 188,814,226 raw reads. After removing adapters, primers and low-quality reads, we obtained a total of 182,618,614 clean reads for the six samples. The sequencing data had an average Q20 of 96.71% and GC content of 46.43% (Additional file 1: Table S1). The clean reads were then assembled into transcripts using Trinity . In the assembly of R. omeimontis, 146,672 transcripts with a total length of 162,215,340 bp were obtained (Table 1). These transcripts were finally assembled into 84,917 unigenes with a total length of 69,985,136 bp. The N50 lengths of the transcripts and unigenes were 2185 bp and 1631 bp, respectively. Most of the unigenes (51,358, 60.48%) were 200–500 bp in length (Additional file 2: Figure S1). A detailed summary of the transcriptome data can be found in Table 1. To verify the assembly quality, we mapped the clean reads to the assembled unigenes, and the mapping rates for the six samples ranged from 77.62 to 81.92% (Additional file 1: Table S1). Taken together, these metrics indicate the high quality of our transcriptome data, which is the foundation for all subsequent analyses . Importantly, this study represents the first transcriptome exploration of oviducts in R. omeimontis and provides valuable gene expression information that can be used to further understand the genetic basis of foam nest production in this species.
Functional annotation and classification of the unigenes
To explore the potential function of the unigenes in the R. omeimontis oviduct during foam nest production, we annotated the unigenes against the NR, Swiss-Prot, Pfam, GO, Clusters of EuKaryotic (KOG) and Kyoto Encyclopedia of Genes and Genomes Orthology (KO) databases. Of the 84,917 unigenes, 28,681 unigenes (33.77%) presented a positive match against at least one database (Table 2). Thus, many of the unigenes have not been annotated, which may be partly due to the limited genetic resources for non-model amphibian species, especially genetic data on expression in the oviduct tissue. Furthermore, unigenes with partial or misassembled transcripts, sequences from untranslated regions (UTRs) of proteins and meaningless proteins could not have a match to any genes. These reasons may explain the annotation failure for those unigenes. The summary annotation information of R. omeimontis unigenes was provided in Table 2.
In GO classification, a total of 22,133 unigenes (26.06%) were assigned to 55 sub-categories of GO terms belonging to the three main categories: Biological Process, Cellular Component and Molecular Function (Fig. 2). The category Biological Process contained 22 sub-categories, and cellular process (64.66%) and metabolic process (53.40%) were the two main enriched GO terms in this category (Fig. 2). The Cellular Component category included 17 sub-categories, and the two most enriched GO terms were cell (42.11%) and cell part (42.05%). In category Molecular Function, 11 sub-categories were found, and the most enriched term was binding (61.78%), followed by the term of catalytic activity (41.51%; Fig. 2).
The COG database mainly contains the phylogenetic classification of proteins encoded in 21 complete genomes of bacteria, archaea and eukaryotes . In the eukaryote-specific COG (KOG) classification, a total of 11,413 unigenes (13.44%) were functionally classified into 25 KOG categories (Fig. 3). Among all KOG classifications, General function prediction only (R) occupied the largest proportion (23.04%), followed by the Signal transduction mechanisms (S; 19.55%) and Posttranslational modification, protein turnover, chaperones (O; 9.66%).
To further understand the gene functions involved in the metabolic pathways, we annotated all unigenes against the KEGG database. A total of 10,735 unigenes (12.64%) were mapped into 263 signaling pathways belonging to 32 categories of pathway hierarchy level 2 and five categories of hierarchy level 1 (Fig. 4). Among the 32 pathways, Signal transduction (13.28%), Endocrine system (6.52%) and Immune system (6.36%) were the most enriched pathways. These annotations will provide valuable resources for further studies of R. omeimontis.
Analysis of differentially expressed genes
To explore DEGs between the two periods of foam nest construction, we estimated the gene expression levels of each sample by RSEM and calculated the FPKM value. The whole transcriptomes of oviduct between the different periods were first analyzed by a principal component analysis. The first principal component (PC1, accounting for 55.48% of the variation) and the second principal component (PC2, accounting for 24.07% of the variation) together explained 79.55% of the variance in the transcriptome data (Fig. 5A). Principal component analysis revealed strong clustering between the two periods of foam nest construction. In addition, the global gene expression profiles of oviduct from the two periods also showed different patterns. The final set of DEGs in oviduct between the two periods was revealed by the hierarchical clustering of expression patterns (Fig. 5B). In total, we identified 433 DEGs in the oviduct of R. omeimontis, among which 270 genes were upregulated and 163 genes were downregulated during the AFNP_O period compared with the BFNP_O period (Fig. 6; Additional file 1: Table S2, Table S3). Among all DEGs, 259 (59.82%) were successfully annotated. To confirm the differential expression profiles obtained with the transcriptome data, we selected a total of 14 genes from DEGs for qRT-PCR analysis. The result revealed that 12 out of 14 genes closely matched the results of RNA-seq, with a correlation coefficient R of 0.92 (Fig. 7).
To better understand the potential biological functions and metabolic pathways that the DEGs were involved in during the two periods of foam nest production, we compared all DEGs to the GO and KEGG databases and applied the significant enrichment analysis. In the upregulated DEGs, extracellular space (GO:0005615) and extracellular region (GO:0005576) were the significantly enriched GO terms (Additional file 1: Table S5; Additional file 2: Figure S2). However, no significant enrichment of GO terms was identified in the downregulated DEGs (Additional file 1: Table S6). Furthermore, all 433 DEGs were classified into 46 KEGG pathways. Among these pathways, the ‘two-component system’, ‘cell adhesion molecules’, ‘steroid hormone biosynthesis’, ‘neuroactive ligand-receptor interaction’, ‘nicotinate and nicotinamide metabolism’ and ‘cytokine-cytokine receptor interaction’ pathways were significantly enriched (Fig. 8). Generally, these significantly enriched pathways were mostly basic metabolism-related pathways. In contrast, in the upregulated DEGs, the significantly enriched KEGG pathways were ‘two-component system’, ‘nicotinate and nicotinamide metabolism’, ‘aminobenzoate degradation’, ‘steroid hormone biosynthesis’ and ‘folate biosynthesis’ (Additional file 1: Table S7; Fig. 8). These DEGs were mostly related to biological processes of reproduction and sex hormone biosynthesis. In the downregulated DEGs, the significantly enriched KEGG pathways were ‘cell adhesion molecules’, ‘phagosome’, ‘cytokine-cytokine receptor interaction’ and ‘lysosome’ (Additional file 1: Table S8; Fig. 8). Therefore, the downregulated DEGs were predominantly participate in pathways with immunity-associated functions.
We further identified the top ten upregulated and top ten downregulated annotated DEGs in the oviduct based on the order of FDRs to obtain potential genes that may directly control the construction or composition of the foam nest (Table 3, Additional file 1: Table S2, Table S3). Among these top 20 annotated DEGs, ten genes were immune-related genes, including REG4, CPAMD9, LYZ2, IL17B, PARP14, CD116, MHC I and C7 (Table 3). Of the top ten upregulated DEGs, three were Regenerating islet-derived protein 4-like genes (REG4, Table 3). And VMO1 gene was the most differentially expressed (log2FC = 10.83, FDR = 4.99E-15, Table 3) of the upregulated genes. Among the top ten downregulated genes, two genes were related to lipid metabolism, including APOC-1 and CYP11A1 gene. And four DEGs were related to immune responses, including MHC I, PARP14, C7 and CD116 gene. In addition, mmp18 gene were also among the top ten downregulated genes.
The production of foam nests is one of the strategies that has evolved to allow some anuran species to protect their eggs and larvae. Despite considerable knowledge of the biochemical components of anuran foam nests, little is known about the molecular basis involved in the construction of foam nest [5,6,7]. In this study, we performed a transcriptome of oviduct in female R. omeimontis from two different periods during foam nest production to explore genes and pathways that may play important roles in this process.
Differential expression analysis revealed 433 DEGs between the two periods during foam nest production. Of these genes, 270 were upregulated and the remaining 163 genes were downregulated. The results from GO and KEGG analysis indicated that various numbers of DEGs were involved in different processes of extracellular metabolism and steroid hormone biosynthesis. The upregulated DEGs were mostly related to biological processes of sex hormone biosynthesis. Before foam nest production began, the females were in the active breeding period; therefore, the significantly enriched pathways in the oviduct tissue were associated with reproduction and sex hormones, such as ‘folate biosynthesis’ and ‘steroid hormone biosynthesis’. Folate, a kind of vitamin B, is an indispensable element of DNA methylation that plays a key role during the process of embryogenesis. A shortage of folic acid may induce neural tube defects [32, 33].
The downregulated DEGs were predominantly participate in pathways with immunity-associated functions. Phagosomes, lysosomes and proteasomes play a part in antigen processing and presentation [34,35,36,37]. Apoptosis is important for the development and normal functioning of the immune system . The TNF signaling pathway, Jak-STAT signaling pathway, NF-kappa B signaling pathway, and PI3K-Akt signaling pathway are activated during the process of immune response, and the cytokine-cytokine receptor interaction contributes to the mediator roles of these signaling pathways [39,40,41,42,43,44,45].
Among the top 20 annotated DEGs, nine genes were immune-related genes, indicating that immune function plays a key role in the production of the leaf foam nest. Of the top ten upregulated DEGs, three were REG4 gene. REG4 is a member of the regenerating (Reg) gene family, which belongs to the calcium-dependent lectin gene superfamily . Lectins are a kind of carbohydrate-binding family proteins that are widely found in viruses, bacteria, plants, animals and humans [47, 48]. Studies on the túngara frog (Engystomops pustulosus) indicated that the foam nest was a mixture of six proteins, named ranaspumins (Rsn-1 to Rsn-6); four of these proteins exhibit lectin activity and may play important roles in forming long-term foam stabilization by contributing to the formation of a cross-linked carbohydrate network at the air-water interface and defending against pathogen attacks . Rsn-3 to Rsn-5 have amino acid sequences similar to those of fucolectins that predominately bind to fucose, whereas Rsn-6 belongs to the C-type lectins that are always associated with galactose binding . Lectins bind to molecular patterns of pathogens, and they hinder microbial dissemination by agglutinating granules with a specific array of carbohydrates . In addition, lectins participate in the process of embryonic development and intracellular glycoprotein transport in some pathways [50,51,52]. Therefore, the lectins may be the constituent proteins of the leaf foam nest of R. omeimontis, and this component was also found in the water-based foam nest of the túngara frog .
The VMO1 gene was the most differentially expressed of the upregulated genes. A previous study revealed that VMO1 reduced the surface tension of tears and produced a more stable tear film in vivo . Thus, we considered that VMO1 in the oviduct of R. omeimontis may have a similar function that could reduce the surface tension of foams and produce a more stable and biocompatible foam nest environment. However, compared to other anuran species with a water-interface foam nest, R. omeimontis did not exhibit similar genes encoding protein surfactants that could reduce surface tension and enable foam nest formation, such as Rsn-1, Rsn-2, ranasmurfin gene and Lv-Rsn-1 [17, 54,55,56,57]. Considering that the composition of the foam nest should be species-specific, the arboreal foam nest constructed by R. omeimontis may have different types of surfactant proteins that differ from those in aquatic foam nests. Moreover, published genetic data for anuran species on the expression in the oviduct are still limited; therefore, many of our unigenes have not been functionally annotated, and these genes may include surfactant protein genes or other important genes associated with foam nest formation in R. omeimontis.
Of the top ten downregulated genes, two genes were related to lipid metabolism, including APOC-1 and CYP11A1 gene. APOC-1 encodes a component of lipoproteins that plays a central role in high density lipoprotein (HDL) and very low-density lipoprotein (VLDL) metabolism of vertebrates [58, 59]. Previous studies suggest that overexpression of APOC-1 in humans and transgenic mice may cause high levels of triglycerides (TGs) as a consequence of APOC-1 inhibition activity of lipoprotein lipase (LPL)-mediated TG lipolysis [60, 61]. LPL hydrolyses the TGs into free fatty acids (FFAs), which are used for energy metabolism or storage in tissues [62,63,64]. The reproductive process in anurans is an energy-consuming activity, especially during the amplexus and spawning period [65, 66]. Additionally, the foams are energetically expensive to make and difficult to maintain, with a tendency to collapse over time unless stabilized mechanically or kinetically by additional processes [56, 67]. Based on our field observation, female R. omeimontis built their foam nest by releasing cloacal mucus fluid and whipping the blend into foam by churning her back legs, and this process commonly takes two to three hours. Thus, we hypothesized that foam nest building is an energy-consuming process for female R. omeimontis, and the highly expressed lipid metabolism-related genes during AFNP_O period may be involved in this process.
CYP11A1 encodes a member of the cytochrome P450 superfamily catalyzing the conversion of cholesterol to pregnenolone . This is the first reaction in the process of steroidogenesis in all mammalian tissues that specialize in the production and metabolism of various steroid hormones [68, 69]. Matrix metalloproteinases (MMPs) play an important role in extracellular matrix remodeling and degradation during organ development and pathological processes . In Xenopus laevis, mmp18 may play a role in larval tissue degeneration and adult organogenesis during metamorphosis [71, 72]. The highly expressed CYP11A1 and mmp18 during AFNP_O period in the oviduct of female R. omeimontis may be important for regulating foam-nesting activity through the steroid hormone pathway and extracellular matrix proteolysis.
Additionally, four immune-related genes were also found in the top ten downregulated genes. Major histocompatibility complex (MHC) molecules vitally participate in immune response by presenting antigens to T cells [73, 74]. The PARP14 gene encodes a protein of poly (ADP-ribose) polymerase family member 14. Poly (ADP-ribose) polymerase is an immediate DNA damage-dependent post-translational modifier of histones and other nuclear proteins that contributes to the survival of injured proliferating cells [75, 76]. The C7 gene encodes a serum glycoprotein that forms a membrane attack complex together with the complement components C5b, C6, C8, and C9 as part of the terminal complement pathway of the innate immune system . CD116 is a receptor for granulocyte-macrophage colony-stimulating factor, which stimulates the production of white blood cells . In anuran amphibians, embryos and larvae are vulnerable to bacterial and fungal pathogen infections [79,80,81]. To protect the offspring, females would maternally transfer antibodies to offspring for immune defense before the larvae are immune at maturity [82, 83]. Therefore, the highly expressed immune-related genes in the oviduct of R. omeimontis could be important in the production of foam nests for protecting the eggs and embryos.
Overall, the above-mentioned candidate genes from the oviduct of R. omeimontis encoded and regulated a range of proteins with a mixture of surfactant, carbohydrate-binding and immune defense activities that together provide a stable, biocompatible, protective foam environment for developing eggs and embryos. The steroid hormone pathway and lipid metabolism may play important roles in regulating the process of foam nest construction.
In this study, transcriptome analysis was performed with the oviduct tissue of R. omeimontis to explore the molecular basis of leaf foam nest construction. Our transcriptome data provided a rich list of genes expressed during the two different periods of foam nesting. A total of 433 DEGs were identified, with 270 genes being upregulated and 163 genes being downregulated. In summary, numerous functional genes and metabolic pathways were likely associated with the process of foam nest construction. During the period when foam nest production began, genes encoding lectin, surfactant protein and immune components were highly expressed, as their products were indispensable components of the foam nest. During the period when the foam nest production was finished, genes related to lipid metabolism, steroid hormone and immune function were highly expressed, and these genes may be important for regulating the process of foam nesting. Our findings will provide useful information for further investigation of anuran foam nest construction at the molecular level.
Sample collection and RNA extraction
The female R. omeimontis individuals were obtained from the population in Badagong Mountain Natural Reserve in Hunan Province, central China (29°47′02″ N, 110°05′27″ E), during the breeding season. The sampling was performed during the period when foam nest production began (BENP_O) and the period after foam nest production was finished (AFNP_O). During each sampling period, three females were caught from groups in amplexus and designated as biological replicates. In BENP_O period, we sampled the oviducts of females immediately when the female that was gripped by a male arrived at the breeding and spawning tree. This period could represent the preparation or beginning stage of foam nest production. In the AFNP_O period, we sampled the oviducts of females immediately when the male left the amplexus. Commonly, females would stay with the foam nest for about half an hour, and continue whipping the foam and wrapping the foam nest with surrounding leaves. Thus, this period could represent the active phase of foam nest construction. In both sampling procedures, we caught females as soon as possible and kept them in individual plastic boxes containing leaves from their natural habitats. The living frogs were euthanized by injecting buffered MS-222 solution and dissected immediately with 0.1% diethyl pyrocarbonate (DEPC)-ddH2O solution-treated scissors. The tissue samples were taken from the oviduct of the frogs and stored in liquid nitrogen until needed.
Total RNA from each oviduct sample was extracted using an RNA extraction kit (Omega Bio-Tek, USA) according to the manufacturer’s instructions. The quality of the total RNA was assessed by 1% gel electrophoresis, and the RNA purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). We measured the RNA concentration using a Qubit® RNA Assay Kit in a Qubit® 2.0 fluorometer (Life Technologies, CA, USA). RNA integrity was checked using a RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Total RNA was extracted separately from each oviduct sample for three replicates in each studied period, and six high-quality RNA samples were obtained.
cDNA library construction and Illumina sequencing
A total of 3 μg of RNA pre-sample was used as input material for cDNA library construction. cDNA libraries were generated using an NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations. In brief, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragments were generated using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H−). Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. To select cDNA fragments of preferentially 150~200 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). PCR was then performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer. PCR products were purified using the AMPure XP system, and library quality was assessed with the Agilent Bioanalyzer 2100 system.
The cDNA samples were then clustered on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina). We sequenced the cDNA libraries on an Illumina HiSeq 2000 platform and paired-end reads were generated, with each read length of 50 bp.
De novo assembly and gene functional annotation
Raw data (raw reads) were firstly processed with in-house Perl scripts (Additional file 3), and clean data (clean reads) were then obtained by removing reads containing adapters, reads containing poly-N and reads of low quality (reads with more than 5% unknown nucleotides or reads of less than 13 bp) from the raw data. The number of raw reads, clean reads and Q20 and Q30 values were calculated at the same time. All subsequent analyses were based on the clean data of high quality. Due to the absence of reference sequences for R. omeimontis, high-quality clean data were de novo assembled into transcripts using Trinity software with min_kmer_cov set to 2 by default and all other parameters set to the default . Unigenes were selected from the longest transcript copy of each gene cluster to avoid redundant transcripts . All unigenes were functionally annotated against the following databases: the NCBI non-redundant protein sequences (NR) database, the Protein family database (Pfam) database (http://pfam.xfam.org/), the manually annotated and reviewed protein sequence database (Swiss-Prot) database, the Clusters of Orthologous Groups of Proteins (COG) database, the Pathway Annotation of the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and the Gene Ontology (GO) database. GO annotation was performed using the Blast2GO method [85, 86].
Differential expression analyses
Gene expression levels were quantified by RSEM for each sample . First, clean data were mapped back onto the assembled transcriptome. 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). Differential expression analysis between the two periods of foam nest construction was performed using the DESeq R package , and the resulting p values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate (FDR) . The genes with an adjusted p value < 0.05 and |log2 (fold change) | > 1 were designated as significant differentially expressed genes (DEGs).
GO enrichment of the DEGs was performed with the goseq R packages based on Wallenius’ non-central hyper-geometric distribution . The KEGG database is a database resource for understanding high-level functions and utilities of biological systems, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/) . KEGG enrichment of the DEGs was applied using KOBAS software and the ‘fdr’ parameter set to ‘BH’ .
To validate the accuracy of our transcriptome data, we selected 14 genes from the significantly differentially expressed genes and performed qRT-PCR using a Bio-Rad CFX96™ Real-Time System. Total RNA was extracted from the oviducts for the two periods using TRIzol reagent. The PCRs were conducted 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 genes were normalized using the housekeeping gene GAPDH as an endogenous control, and the relative gene expression was calculated by the 2-ΔΔCT method. All experiments were performed with three biological replicates, and the data were analysed by one-way analysis of variance (ANOVA). The primers used in this study were designed with Primer Express 3.0 software, and additional information about the primers is provided (Additional file 1: Table S4).
Availability of data and materials
The dataset supporting the conclusions of this study is included within the article and its additional files. The raw data of this study has been deposited in the NCBI Sequence Read Archive (SRA) database under the accession PRJNA529065. The assembled sequence has been deposited at DDBJ/EMBL/GenBank under the accession GHMZ00000000. The version described in this paper is the first version, GHMZ01000000.
Apolipoprotein C-I-like gene
Complement component 7
DNA complementary to RNA
Cytochrome P450 family 11 subfamily A polypeptide 1
Differentially expressed genes
False discovery rate adjusted p value
Gene Ontology database
Kyoto Encyclopedia of Genes and Genomes database
Clusters of Orthologous Groups of proteins database
Major histocompatibility complex
non-redundant protein sequences database
Poly [ADP-ribose] polymerase 14 gene
Protein family database
Quantitative reverse transcription real-time PCR
Regenerating islet-derived protein gene
the manually annotated and reviewed protein sequence database
Vitelline membrane outer layer protein 1 homolog precursor
Andrade DV, Abe AS. Foam nest production in the armoured catfish. J Fish Biol. 1997;50:665–7.
Cooper A, Kennedy MW, Fleming RI, Wilson EH, Videler H, Wokosin DL, et al. Adsorption of frog foam Nest proteins at the air-water Interface. Biophys J. 2005;88:2114–25.
Castilla JC, Manriquez PH, Delgado AP, Gargallo L, Leiva A, Radic D. Bio-foam enhances larval retention in a free-spawning marine tunicate. Proc Natl Acad Sci U S A. 2007;104:18120–2.
Dalgetty L, Kennedy MW. Building a home from foam--tungara frog foam nest architecture and three-phase construction process. Biol Lett. 2010;6:293–6.
Haddad C, Prado C. Reproductive modes in frogs and their unexpected diversity in the Atlantic Forest of Brazil. Bioscience. 2005;55:207–17.
Wells KD. The ecology and behavior of amphibians: University of Chicago Press; 2010.
Iskandar DT, Evans BJ, McGuire JA. A novel reproductive mode in frogs: a new species of fanged frog with internal fertilization and birth of tadpoles. PLoS One. 2014;9:e115884.
Prado C, Toledo LF, Zina J, Haddad C. Trophic eggs in the foam nests of Leptodactylus labyrinthicus (Anura, Leptodactylidae): an experimental approach. Herpetol J. 2005;15:279–84.
Shepard DB, Caldwell JP. From foam to free-living: ecology of larval Leptodactylus labyrinthicus. Copeia. 2005;2005:803–11.
Hissa DC, Bezerra WM, Freitas CD, Ramos MV, Lopes JL, Beltramini LM, et al. Frog foam Nest protein diversity and synthesis. J Exp Zool. 2016;325:425–33.
Downie JR. Functions of the foam in the foam-nesting leptodactylid physalaemus-pustulosus. Herpetol J. 1998;1:302–7.
Menin M, Giaretta AA. Predation on foam nests of leptodactyline frogs (Anura: Leptodactylidae) by larvae of Beckeriella Niger (Diptera: Ephydridae). J Zool. 2003;261:239–43.
Bastos RP, Haddad CFB, Pombal JP Jr. Foam nest in Scinax rizibilis (Amphibia: Anura: Hylidae). Zoologia (Curitiba). 2010;27:881–6.
Faivovich J, Ferraro DP, Basso NG, Haddad C, Rodrigues MT, Wheeler WC, et al. A phylogenetic analysis of Pleurodema (Anura: Leptodactylidae: Leiuperinae) based on mitochondrial and nuclear gene sequences, with comments on the evolution of anuran foam nests. Cladistics. 2012;28:460–82.
Heyer WR, Rand AS. Foam nest construction in the leptodactylid frogs Leptodactylus pentadactylus and Physalaemus pustulosus (Amphibia, Anura, Leptodactylidae). J Herpetol. 1977;11:225–8.
Hissa DC, Vasconcelos IM, Carvalho AF, Nogueira VL, Cascon P, Antunes AS, et al. Novel surfactant proteins are involved in the structure and stability of foam nests from the frog Leptodactylus vastus. J Exp Biol. 2008;211:2707–11.
Fleming RI, Mackenzie CD, Cooper A, Kennedy MW. Foam nest components of the tungara frog: a cocktail of proteins conferring physical and biological resilience. Proc Biol Sci. 2009;276:1787–95.
Alcaide MF, Lavilla EO, Alcaide AP. Histology and Histochemistry of the albumin glands in some foam-nesting anurans. South American Journal of Herpetology. 2009;4:151–63.
Furness AI, McDiarmid RW, Heyer WR, Zug GR. Oviduct modifications in foam-nesting frogs, with emphasis on the GenusLeptodactylus (Amphibia, Leptodactylidae). S Am J Herpetol. 2010;5:13–29.
Cooper A, Vance SJ, Smith BO, Kennedy MW. Frog foams and natural protein surfactants. Colloids Surf A Physicochem Eng Asp. 2017;534:120–9.
Liao WB, Lu X. Breeding behaviour of the Omei tree frogRhacophorus omeimontis (Anura: Rachophoridae) in a subtropical montane region. J Nat Hist. 2010;44:2929–40.
Luo Z, Li C, Wang H, Shen H, Zhao M, Gu Q, et al. Male-male competition drives sexual selection and group spawning in the Omei treefrog, Rhacophorus omeimontis Behav. Ecol Sociobiol. 2016;70:593–605.
Zhao M, Li C, Zhang W, Wang H, Luo Z, Gu Q, et al. Male pursuit of higher reproductive success drives female polyandry in the Omei treefrog. Anim Behav. 2016;111:101–10.
Liao WB, Lu X. A comparison of reproductive output of the Omei Treefrog (Rhacophorus omeimontis) between high and low elevations. Anim Biol. 2011;61:263–76.
Liao WB, Lu X. Variation in body size, age and growth in the Omei treefrog (Rhacophorus omeimontis) along an altitudinal gradient in western China. Ethol Ecol Evol. 2011;23:248–61.
Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57–63.
Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet. 2011;12:87–98.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bioinformatics. 2011;12:323.
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.
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.
Tatusov RL, Fedorova ND, Jackson JD, Jacobs AR, Kiryutin B, Koonin EV. The COG database: an update version includes eukaryotes. BMC Bioinformatics. 2003;4:41.
Christensen B, Rosenblatt DS. 9 effects of folate deficiency on embryonic development. Baillière's clinical haematology. 1995;8:617–37.
McKay JA, Williams EA, Mathers JC. Folate and DNA methylation during in utero development and aging. Biochem Soc Trans. 2004;32:1006–7.
Garin J, Diez R, Kieffer S, Dermine J-F, Duclos S, Gagnon E, et al. The phagosome proteome. J Cell Biol. 2001;152:165–80.
Houde M, Bertholet S, Gagnon E, Brunet S, Goyette G, Laplante A, et al. Phagosomes are competent organelles for antigen cross-presentation. Nature. 2003;425:402–6.
Kloetzel P-M. The proteasome and MHC class I antigen processing. BBA-Mol Cell Res. 1695;2004:225–33.
Sijts EJ, Kloetzel PM. The role of the proteasome in the generation of MHC class I ligands and immune responses. Cell Mol Life Sci. 2011;68:1491–502.
Opferman JT, Korsmeyer SJ. Apoptosis in the development and maintenance of the immune system. Nat Immunol. 2003;4:410–5.
Li Q, Verma IM. NF-kappaB regulation in the immune system. Nat Rev Immunol. 2002;2:725–34.
Silverman N, Maniatis T. NF-κB signaling pathways in mammalian and insect innate immunity. Genes Dev. 2001;15:2321–42.
Chen G, Goeddel DV. TNF-R1 signaling: a beautiful pathway. Science. 2002;296:1634–5.
Agaisse H, Perrimon N. The roles of JAK/STAT signaling in Drosophila immune responses. Immunol Rev. 2004;198:72–82.
Hazeki K, Nigorikawa K, Hazeki O. Role of phosphoinositide 3-kinase in innate immunity. Biol Pharm Bull. 2007;30:1617–23.
Sauer S, Bruno L, Hertweck A, Finlay D, Leleu M, Spivakov M. T cell receptor signaling controls Foxp3 expression via PI3K, Akt, and mTOR. P Natl Acad Sci. 2008;105:7797–802.
Juntilla MM, Koretzky GA. Critical roles of the PI3K/Akt signaling pathway in T cell development. Immunol Lett. 2008;116:104–10.
Chen S, Gou WF, Zhao S, Niu ZF, Zhao Y, Takano Y, et al. The role of the REG4 gene and its encoding product in ovarian epithelial carcinoma. BMC Cancer. 2015;15:471.
Zhu-Salzman K, Shade RE, Koiwa H, Salzman RA, Narasimhan M, Bressan RA. Carbohydrate binding and resistance to proteolysis control insecticidal activity of Griffonia simplicifolia lectin II. P Natl Acad Sci. 1998;95:15123–8.
Murdock LL, Shade RE. Lectins and protease inhibitors as plant defenses against insects. J Agr Food Chem. 2002;50:6605–11.
Gupta G, Surolia A. Collectins: sentinels of innate immunity. Bioessays. 2007;29:452–64.
Hauri HP, Appenzeller C, Kuhn F, Nufer O. Lectins and traffic in the secretory pathway. FEBS Lett. 2000;476:32–7.
Li Y, Chia JM, Bartfai R, Christoffels A, Yue GH, Ding K, et al. Comparative analysis of the testis and ovary transcriptomes in zebrafish by combining experimental and computational tools. Comp Funct Genom. 2004;5:403–18.
Masse K, Baldwin R, Barnett MW, Jones EA. X-epilectin: a novel epidermal fucolectin regulated by BMP signalling. Int J Dev Biol. 2004;48:1119–29.
Wang Z, Chen ZY, Yang Q, Jiang YB, Lin LP, Liu XL. Vitelline membrane outer layer 1 homolog interacts with lysozyme C and promotes the stabilization of tear film. Invest Ophth Vis Sci. 2014;55:6722–7.
Oke M, Ching RT, Carter LG, Johnson KA, Liu H, McMahon SA, et al. Unusual chromophore and cross-links in ranasmurfin: a blue protein from the foam nests of a tropical frog. Angew Chem Int Ed Engl. 2008;47:7853–6.
Mackenzie CD, Smith BO, Meister A, Blume A, Zhao X, Lu JR, et al. Ranaspumin-2: structure and function of a surfactant protein from the foam nests of a tropical frog. Biophys J. 2009;96:4984–92.
Cooper A, Kennedy MW. Biofoams and natural protein surfactants. Biophys Chem. 2010;151:96–104.
Cavalcante Hissa D, Arruda Bezerra G, Birner-Gruenberger R, Paulino Silva L, Uson I, Gruber K, et al. Unique crystal structure of a novel surfactant protein from the foam nest of the frog Leptodactylus vastus. Chembiochem. 2014;15:393–8.
Li WH, Tanimura M, Luo CC. The apolipoprotein multigene family: biosynthesis, structure, structure-function relationships, and evolution. J Lipid Res. 1988;29:245–71.
Kim KY, Cho YS, Bang IC, Nam YK. Isolation and characterization of the apolipoprotein multigene family in Hemibarbus mylodon (Teleostei: Cypriniformes). Comp Biochem Physiol B Biochem Mol Biol. 2009;152:38–46.
Berbee JF, Van der Hoogt CC, Sundararaman D, Havekes LM, Rensen PC. Severe hypertriglyceridemia in human APOC1 transgenic mice is caused by apoC-I-induced inhibition of LPL. J Lipid Res. 2005;46:297–306.
Westerterp M, Berbee JF, Delsing DJ, Jong MC, Gijbels MJ, Dahlmans VE, et al. Apolipoprotein C-I binds free fatty acids and reduces their intracellular esterification. J Lipid Res. 2007;48:1353–61.
Peterson J, Bihain BE, Bengtsson-Olivecrona G, Deckelbaum RJ, Carpentier YA, Olivecrona T. Fatty acid control of lipoprotein lipase: a link between energy metabolism and lipid transport. P NATL ACAD SCI. 1990;87:909–13.
Andersson Y, Majd Z, Lefebvre A-M, Gv M, Sechkin AV, Kosykh V, et al. Developmental and pharmacological regulation of apolipoprotein C-II gene expression. Arterioscl Throm Vas. 1999;19:115–21.
Murata M, Tojo S. Utilization of lipid for flight and reproduction in Spodoptera litura (Lepidoptera: Noctuidae). Eur J Entomol. 2002;99:221–4.
Ryser J. Weight loss, reproductive output, and the cost of reproduction in the common frog. Rana temporaria Oecologia. 1989;78:264–8.
Zhou W, Li X, Wu JP, Li LF, Li MH. Changes in reproductive and accompanying organs of Rana chaochiaoensis in the Kunming area. Zool Res. 2008;29:89–94.
Clarkson JR, Cui ZF, Darton RC. Protein denaturation in foam: II. Surface activity and conformational change. J Colloid Interf Sci. 1999;215:333–8.
Hanukoglu I. Steroidogenic enzymes: structure, function, and role in regulation of steroid hormone biosynthesis. J Steroid Biochem Mol Biol. 1992;43:779–804.
He G, Xu W, Chen Y, Liu X, Xi M. Abnormal apoptosis of trophoblastic cells is related to the up-regulation of CYP11A gene in placenta of preeclampsia patients. PLoS One. 2013;8:e59609.
Snoek-van Beurden PAM, Von den Hoff JW. Zymographic techniques for the analysis of matrix metalloproteinases and their inhibitors. Biotechniques. 2005;38:73–83.
Stolow MA, Bauzon DD, Li J, Sedgwick T, Ling V, Sang QA. Identification and characterization of a novel collagenase in Xenopus laevis: possible roles during frog development. Mol Biol Cell. 1996;7:1471–83.
Ishizuya-Oka A, Hasebe T, Shi YB. Apoptosis in amphibian organs during metamorphosis. Apoptosis. 2010;15:350–64.
Benacerraf B. Role of MHC gene products in immune regulation. Science. 1981;212:1229–38.
Messaoudi I, Guevara Patino JA, Dyall R, LeMaoult J, Nikolich-Zugich J. Direct link between mhc polymorphism, T cell avidity, and diversity in immune defense. Science. 2002;298:1797–800.
Aguiar RC, Takeyama K, He C, Kreinbrink K, Shipp MA. B-aggressive lymphoma family proteins have unique domains that modulate transcription and exhibit poly (ADP-ribose) polymerase activity. J Biol Chem. 2005;280:33756–65.
Iwata H, Goettsch C, Sharma A, Ricchiuto P, Goh WW, Halu A, et al. PARP9 and PARP14 cross-regulate macrophage activation via STAT1 ADP-ribosylation. Nat Commun. 2016;7:12849.
Ying LS, Zhang FR, Pan XD, Chen KY, Zhang N, Jin JY. Complement component 7 (C7), a potential tumor suppressor, is correlated with tumor progression and prognosis. Oncotarget. 2016;7:86536.
Metcalf D. Binding of 125I-labeled granulocyte colony-stimulating factor to normal murine hemopoietic cells. J Cell Physiol. 1985;124:313–21.
Gomez-Mestre I, Touchon JC, Warkentin KM. Amphibian embryo and parental defenses and a larval predator reduce egg mortality from water mold. Ecology. 2006;87:2570–81.
Hamdoun A, Epel D. Embryo stability and vulnerability in an always changing world. Proc Natl Acad Sci U S A. 2007;104:1745–50.
Gomez-Mestre I, Touchon JC, Saccoccio VL, Warkentin KM. Genetic variation in pathogen-induced early hatching of toad embryos. J Evol Biol. 2008;21:791–800.
Poorten TJ, Kuhn RE. Maternal transfer of antibodies to eggs in Xenopus laevis. Dev Comp Immunol. 2009;33:171–5.
Walke JB, Harris RN, Reinert LK, Rollins-Smith LA, Woodhams DC. Social immunity in amphibians: evidence for vertical transmission of innate defenses. Biotropica. 2011;43:396–400.
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.
Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:1.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.
Storey JD. The positive false discovery rate: a Bayesian interpretation and the q-value. Ann Stat. 2003;31:2013–35.
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. 2007;36:480–4.
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.
We would like to thank Chenling Li for help with the sample collection for this study. We also thank the staff of the Badagongshan National Nature Reserve for assistance in field investigation and sample collection.
This study was supported by the National Natural Science Foundation of China (No. 31470442). The funder provided financial support in samples collection and the transcriptomic sequencing. It did not directly participate to the design of the study, analysis, interpretation of data and writing the manuscript.
Ethics approval and consent to participate
All animals used in the current study were approved by the Institutional Animal Care and Use Committee in Central China Normal University. Care and handling of experimental animals were carried out in accordance with guidelines of Hubei Province Animal Management Regulations (011043145–029–2013-000009).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1-S9. Table S1: The statistics and quality of reads for the six samples of the Rhacophorus omeimontis transcriptome. Table S2: The upregulated differentially expressed genes in the oviduct of Rhacophorus omeimontis. Table S3: The downregulated differentially expressed genes in the oviduct of Rhacophorus omeimontis. Table S4: The primers used for quantitative RT-PCR analysis in this study. Table S5: GO enrichment of the upregualted DEGs. Table S6: GO enrichment of the downregualted DEGs. Table S7: KEGG enrichment of the upregualted DEGs. Table S8: KEGG enrichment of the downregualted DEGs. Table S9: The annotation information of all unigenes. (XLSX, 14,044 kb) (XLSX 14036 kb)
Figure S1-S2. Figure S1. The unigene length distribution of the Rhacophorus omeimontis transcriptome. Figure S2. Significant enrichment analysis of GO terms in the upregulated differentially expressed genes. (DOCX, 119 kb) (DOCX 118 kb)
The Perl scripts for transcriptome analysis in this study. (ZIP 5 kb)