Significant changes in the lipid composition of B. gargarizans during metamorphosis

To ensure high-quality data acquisition, the performance of the LC–MS platform was monitored before and after analyzing each block based on QC samples. The samples were subjected to qualitative analysis, as shown in Fig. 2A and B. The base peak chromatograms (BPI) of samples at four stages were overlapped and compared. The results showed that the chromatographic peak response intensity and retention time of different samples partially overlapped, indicating that the four groups of samples were comparable.

Fig. 2

Mass spectrum of B. gargarizans at different stages. BPI (base peak intensity chromatogram) of positive and negative ions [A: positive ion; B: negative ion]. Multivariate statistical analysis of lipid molecules [C: PCA score chart; D: OPLS-DA score chart; E: model validation diagram for OPLS-DA]

The unsupervised principal component analysis (PCA) method was applied to obtain analysis results that were more intuitive and visual for the different developmental stages. The gathered results indicated that the samples for each stage could be effectively classified and characterized, and the groups grew further apart as the tadpoles developed. QC samples clustered together displayed satisfactory data quality (Fig. 2C). Orthogonal partial least squares-discriminant analysis (OPLS-DA) was used for variables with significant changes in the micromolecules among each group. The parameters of the OPLS-DA score plot obtained from these four groups (Fig. 2D) were R2X = 0.817, R2Y = 0.967, and Q2 = 0.805, which indicated excellent fitness and reliability. The model was validated by the permutation test (Fig. 2E), and the results showed that R2 and Q2 on the left were both smaller than the initial R2 and Q2 on the right, indicating that the model was reliable and could be used for subsequent lipid metabolism analysis [27].

Comparative research of differential lipids in crucial growth stages was performed to elucidate the regulation mechanism of lipid homeostasis for the development of tadpoles. Progenesis QI selected and identified a total of 89 significantly differential lipid molecules (| log2FC |> 1 and adjusted P-value < 0.05) from the comparison between G31 vs. G38, G38 vs. G42, and G42 vs. G46 (Additional file 2: Table S2). There were 48 differential lipids between the G31 group and G38 group: 16 for fatty acyls (including FA and CAR), 17 for LPC, 13 for ST, and 1 for PE and PI; 27 differential lipids were discovered in the G38 and G42 group: 12 for fatty acyls (including FA and CAR), 8 for ST, 5 for LPC, and 1 for PI and PG; 61 lipids were found to be significantly different between the G42 and G46 groups: 14 for fatty acyls (including FA and CAR), 26 for LPC, 18 for ST, and 1 for PE, PG, and PI. To show the variation of differential lipids in tadpoles over time, a heat map was drawn to visualize the tissue distribution (Fig. 3A). Significant changes in the relative quantities of differential lipids at different growth stages suggested the significance and representativeness of the identified lipid metabolites. The clustering of repeats presented satisfactory biological reproducibility; therefore, differential lipids were classified and analyzed by accumulation (Fig. 3B and Additional file 8: Fig. S1). The most significant proportion of ST was measured at the G31 stage, steadily decreased as development proceeded, and included a large number of bile acids, such as apocholic acid (ST 24:2;O4), 7-ketodeoxycholic acid (ST 24:2;O5), glycochenodeoxycholic acid (ST 24:1;O4;G), and ursocholic acid (ST 24:1;O5). The same trend in ST variation was observed for FA and PI, while PG gradually accumulated from the G31 to G46 stage. PE gradually decreased from G31 to G42, but the expression began to callback at the G46 stage, and unlike the others, there was a significant increase or decrease in the relative amounts of CAR and LPC at the G46 stage.

Fig. 3

Analysis of lipidome differences during four stages of B. gargarizans. A Heat map of differential lipid metabolism. B Comparison of relative amounts in the lipid subgroups. C Comparison of the relative quantities of 7 BDs in each group. ****p < 0.0001, ***p < 0.001, **p < 0.001, *p < 0.01

Accumulation patterns of major BDs in B. gargarizans during metamorphosis

By matching retention times and mass-to-charge ratios of the standards, the relative quantities of the major BDs at four key stages were determined from the mass spectrometry data, including cinobufagin, bufalin, resibufogenin, arenobufagin, desacetylcinobufagin, telocinobufagin, and gambufatalin (Fig. 3C). There was a clear accumulation trend, where BDs significantly increased during the entire metamorphosis process except for cinobufagin. Their relative difference in amounts was evident at different stages. In Fig. 3C, arenobufagin was the highest at the G46 stage, followed by desacetylcinobufagin and telocinobufagin, while the lowest amount was for cinobufagin. The shallow level of cinobufagin did not affect the accumulation trend of the total BDs. The corresponding information is summarized in Additional file 3: Table S3.

Gene regulation process of body remodeling in B. gargarizans

To enhance our understanding of the molecular mechanisms involved in the body remodeling regulation of B. gargarizans, transcriptome sequencing was performed using the Illumina NovaSeq 6000 platform. After raw data filtering and examination of the sequencing error rate, 78.20 GB of clean data were obtained, with 6.03–7.18 Gb of clean data and a Q30 score of at 93.62–94.45% for each sample. The GC content was distributed in the range of 46.42–47.49% (Additional file 4: Table S4). A boxplot showed trends in transcript level changes in the 12 samples that indicated that the data distribution of gene expression levels in the samples was the same basically (Additional file 9: Fig. S2A). After quality control, the cleaned sequences were aligned to the B. gargarizans (Asiatic toad) genome (NCBI_ASM1485885v1, reference genome), with a matching rate of 73.24–77.53%, indicating that the transcriptome data could be further analyzed [1].

A correlation heat map (Additional file 9: Fig. S2B) and the comprehensive PCA (Additional file 9: Fig. S2C) demonstrated satisfactory biological repeatability and significant differences between groups. Moreover, genes were functionally annotated in the GO (Gene Ontology) and KEGG databases. Based on the GO database, 34,129 genes were annotated, and based on the KEGG database, 16.202 genes were functionally annotated. Pairwise comparisons identified four stages of differentially expressed genes (DEGs). As shown in Additional file 9: Fig. S2D and E, 6731 (3626 up and 3105 down), 16,497 (7802 up and 8695 down), and 18,821 (9968 up and 8853 down) DEGs were identified by transcriptome sequencing analysis between different growth groups and the control G31 group. Of these DEGs, 1444 genes were consistently upregulated, and 1385 genes were consistently downregulated, indicating that B. gargarizans activated the expression of numerous genes to promote self-development and metamorphosis [28].

The DEGs of G31 vs. G38, G38 vs. G42, and G42 vs. G46 were focused on, including 6731 DEGs (3626 upregulated and 3105 downregulated) of G31 vs. G38, 8627 DEGs (3616 upregulated and 5011 downregulated) of G38 vs. G42, and 9,150 DEGs (5192 upregulated and 3958 downregulated) of G42 vs. G46. The GO enrichment analysis was carried out to elucidate the specific biological functions of the DEGs in the three comparisons (G31 vs. G38, G38 vs. G42, and G42 vs. G46). The top 20 significantly enriched GO terms concerning cellular component (CC), molecular function (MF), and biological process (BP) are shown in Additional file 5: Table S5 and Fig. S3. GO analysis indicated that B. gargarizans could effectively regulate peptidase activity and the specific combination of vitamin D, iron ion, and heme to facilitate optimal muscle content, digestion, and metabolism mode, allowing its body to more optimally adapt to metamorphosis processes [29,30,31,32]. To confirm the biological pathways of metamorphosis, all DEGs were also assigned to the KEGG database for KEGG pathway analysis (Additional file 6: Table S6 and Fig. S4). The results suggested that the changes in lipid metabolism and energy metabolism provided conditions for B. gargarizans as development proceeded.

Moreover, in the present study, the focus was on the detailed functions of the common DEGs (1444 upregulated and 1385 downregulated) from three comparison groups (G31 vs. G38, G31 vs. G42, and G31 vs. G46), as common DEGs were more related to metamorphosis development of B. gargarizans. For upregulated DEGs, GO terms were mainly enriched in sarcomere organization, methylation, and muscle contraction in BP; myosin filament, organelle membrane, and troponin complex in CC; and monooxygenase activity in MF (Additional file 12: Fig. S5A). For downregulated DEGs, GO terms were mainly enriched in digestion, cholesterol catabolic process, and peptide catabolic process in BP; brush border membrane and apical plasma membrane in CC; and serine-type endopeptidase activity and iron ion binding in MF (Additional file 12: Fig. S5B).

To further understand the pathways involved in differentially expressed genes, the KEGG database was used to identify DEG-enriched pathways. The results showed that DEGs were significantly enriched in 74 KEGG pathways (36 for upregulated DEGs, 22 for downregulated DEGs, and 16 for up- and downregulated DEGs). These pathways were divided into five levels: metabolism, human diseases, organic systems, environmental information processing, and cellular processes. Of these, nicotinate and nicotinamide metabolism, the estrogen signaling pathway, Staphylococcus aureus infection, retinol metabolism, fat digestion, and absorption were the top five enriched pathways in upregulated DEGs. In contrast, protein digestion and absorption, pancreatic secretion, primary bile acid biosynthesis, influenza A, and neuroactive ligand-receptor interaction rank in the top five pathways in downregulated DEGs (Additional file 12: Fig. S5C and D).

Gene regulation process of lipid changes in B. gargarizans

To understand the differences in lipid synthesis during metamorphosis (G31–G46), transcriptomics and lipidomics data were integrated for analysis. Pearson analysis was performed to assess the correlation of 89 differential lipids with DEGs from pairwise comparisons (G31 vs. G38, G31 vs. G42, and G31 vs. G46). Differential lipid-related DEGs were shown in Additional file 7: Table S7. The related DEGs with Pearson's correlation coefficient (PCC) greater than 0.8 were selected for KEGG enrichment analysis (Additional file 13: Fig. S6); they were mainly enriched in “primary bile acid biosynthesis (P-value = 4.40E−24),” “PPAR signaling pathway (P-value = 1.10E−23),” “cholesterol metabolism (P-value = 8.20E−17),” “steroid hormone biosynthesis (P-value = 9.49E−17),” “vitamin digestion and absorption (P-value = 9.41E−16)” and other pathways, among which most DEGs were involved in the “PPAR signaling pathway” (ListHit = 37).

Peroxisome proliferator-activated receptors (PPARs) are nuclear hormone receptors that regulate lipid metabolism, cell differentiation, energy metabolism, and inflammatory response [33]. PPARγ can upregulate the expression of the muscle isoform of carnitine palmitoyltransferase I (CPT1) and FA β-oxidation. It can also regulate the metabolic pathway of cholesterol and polyunsaturated fatty acids by changing the expression of CYP7A1/CYP8B1 [34,35,36]. As a critical transcription factor in energy metabolism, PPARγ can provide energy for metamorphosis. In the process of G31-G46, the expression of PPARγ, α, and β/δ significantly decreased; the downregulated retinoid X receptor (RXR) also formed heterodimers PPAR-RXR that can combine with the promoter or heterodimers PPAR-RXR through heterodimerization with PPARα, β/δ, or γ to regulate lipid metabolism, fat formation, and maintain metabolic homeostasis [37]. Fig. 4 demonstrates the primary regulatory process of the PPAR signaling pathway during metamorphosis.

Fig. 4

Differential lipid-related DEGs regulate lipid metabolism in the PPAR signaling pathway. The colors of each cell from left to right in the heatmaps represent the gene expression trends among G31, G38, G42, and G46. Gene expression levels are expressed in FPKM values. Fatty acid binding protein (FABP), retinoid x receptor (RXR), peroxisome proliferator-activated receptor (PPAR), 3-hydroxy-3-methylglutaryl-coa synthase 2 (HMGCS2), apolipoprotein (Apo), phospholipid transfer protein (PLTP), malic enzyme 1 (ME1), stearoyl-CoA desaturase (SCD), cytochrome p450 family 27 (CYP27), cytochrome p450 family 7 subfamily a member 1 (CYP7A1), cytochrome p450 family 8 subfamily b member 1 (CYP8B1), acyl-CoA synthetase (ACS), cd36 antigen (FATCD36), enoyl-CoA hydratase (BIEN), sterol carrier protein 2 (SCP-X), acyl-CoA oxidase (ACO), carnitine O-palmitoyltransferase 1 (CPT1), adipocyte (AP2), adiponectin (ADIPO), sorbin and SH3 domain-containing protein 1 (CAP), matrix metalloproteinase-1 (MMP), ubiquitin C (UBC), and phosphoenolpyruvate carboxykinase (PEPCK)

The DEGs related to differential lipid and BD accumulation were significantly enriched in the PPAR signaling pathway. To validate the reliability of RNA-Seq data, 6 DEGs with differential lipid correlation coefficients greater than 0.9 (| PCC |≥ 0.9) were selected and examined by qRT-PCR. Except for FABPL in G38 and G42, the expression trends of the selected DEGs were consistent with the Illumina sequencing results (Additional file 14: Fig. S7), indicating the reliability and dependability of the RNA-Seq data.

Identification of key pathways and genes for BD accumulation/synthesis in B. gargarizans

A relatively low level of cinobufagin was not considered based on the assay results of BDs. There were significant differences in the relative quantities of other BDs at the G46 stage compared to other critical stages in the metamorphosis. Therefore, G46 was compared with G31, G38, and G42, and identification was achieved for 3260 common DEGs associated with BD accumulation (Additional file 15: Fig. S8A). GO and KEGG enrichment analyses were performed for these DEGs. The top 30 significantly enriched GO terms concerned CC, MF, and BP (Additional file 15: Fig. S8B). The 3260 DEGs with enriched GO terms were mainly involved in digestion, positive regulation of cell division, negative regulation of endopeptidase activity, polysaccharide digestion, polysaccharide catabolic process, maintenance of gastrointestinal epithelium (BP), extracellular space, extracellular region, myosin filament, myofibril, organelle membrane (CC), serine-type endopeptidase inhibitor activity, serine-type endopeptidase activity, heme binding, aromatase activity, metalloendopeptidase inhibitor activity, iron ion binding, motor activity, aspartic-type endopeptidase activity, and peptidase inhibitor activity (MF). The top 30 significantly enriched KEGG terms are shown in Additional file 15: Fig. S8C, with 3260 DEGs mainly involved in protein digestion and absorption, pancreatic secretion, IL-17 signaling pathway, influenza A, neuroactive ligand-receptor interaction, tight junctions, pathogenic Escherichia coli infection, salmonella infection, fat digestion and absorption, and retinol metabolism. The results of these enrichment analyses were too broad, and the range of genes associated with BDs still required narrowing.

To determine the regulatory genes of BD accumulation, Pearson analysis was used to identify DEGs strongly correlated with BDs (|PCC| of ≥ 0.9) to narrow the range. 689 DEGs were identified and subjected to GO and KEGG enrichment analysis. The results showed that DEGs were significantly enriched in pathways such as the “pentose phosphate pathway,” “retinol metabolism,” and “metabolism of xenobiotics by cytochrome p450.” Furthermore, some of the DEGs were enriched in GO terms such as “aromatase activity,” “IL-17 signaling pathway,” “carbohydrate binding,” “iron ion binding,” and “monosaccharide binding” (Fig. 5A and B), which provided a molecular function for further study on the genes related to bufadienolides. “Pentose phosphate pathway” and “aromatase activity” are critical biosynthetic pathways [38, 39], and the DEGs involved in the pathways (GO term) included CYP3A29, CYP2K1, CYP2B4, CYP2D15, CYP2K4, CYP2F3, CYP2C54, CYP2A10, CYP2G1, CYP2H2, PGDH, PFK-M, TKT, ALDA, TALDO1, and GPI.

Fig. 5

Regulatory analysis of BD accumulation. A and B The top 30 GO terms and KEGG pathways of significant enrichment for the BD-related DEGs are shown. C The top 300 DEGs with the highest correlation coefficients were plotted as co-occurrence networks with bufadienolides. The red diamonds represent BDs, and the circles represent BD-related DEGs. Darker colors and larger circles indicate stronger correlations

Pearson analysis was used to evaluate the strength of the correlation between DEGs and BDs, and the top 150 DEGs with the highest correlation coefficients were plotted as co-occurrence networks with BDs (Fig. 5C). In Fig. 5C, the gene correlation coefficient related to resibufagin was higher; desacetylcinobufagin, arenobufagin, and telocinobufagin were more closely related and possessed a more significant number of potential regulatory genes, with bufalin and gamabufatalin having fewer related genes. The accumulation pattern of cinobufagin was different from that of the other BDs. Thus, cinobufagin was not selected for preparing the co-occurrence networks of the common DEGs.

Visualization of the distribution of BDs in B. gargarizans

To understand how BDs are synthesized and accumulated in toads during metamorphosis, DESI-MSI was used to create for the first time a tissue distribution map of BDs in whole-body cryosections of juvenile toads (G46 stage). The workflow of the tissue distribution study of BDs in toads is summarized in Additional file 16: Fig. S9. Two sagittal sections and one coronal section were examined, and light microscopy analysis was used to compare adjacent sections of the sample to determine the BD-enriched organization. Additional file 17: Fig. S10 shows that the organs were well differentiated, indicating that the cryosection method and mass spectrometry imaging conditions for this experiment were suitable. They are feasible for studying the tissue distribution of BDs in this model.

The results of mass spectrometry imaging are shown in Fig. 6, and include desacetylcinobufagin, resibufogenin, bufalin, arenobufagin, and cinobufagin. Although the response strengths of the five BDs were different, they were all enriched in the same tissues. Fig. 6A(1–5) shows that these BDs were concentrated in the liver of toads but not in the immature parotid gland, suggesting that the juveniles, unlike the adults, do not accumulate BDs in the immature parotid gland but produce a large number of BDs in the liver. Not only did Fig. 6B(1–5) support this result, but it was also found that BDs collected in the gallbladder. Cinobufagin was the most obvious in Fig. 6B5. The function of the gallbladder is to receive bile secreted by hepatocytes. The finding that a large number of BDs were enriched in the liver suggested that BDs were likely to be produced by hepatocytes and accumulate in the gallbladder. A horizontal section was also examined, and BDs were still mainly distributed in liver tissue (Fig. 6C). Optical microscopy analysis using a Zeiss microscope (bright field, top-illumination) to observe the position of each organ in the whole-body sections confirmed that these tissues were derived from the liver (Additional file 17: Fig. S10a–c). Additionally, most of the metabolic pathways involved in BD-related DEGs in Fig. 5B are related to the liver. Thus, the liver plays a crucial role in the biosynthesis and accumulation of BDs.

Fig. 6

The results of DESI-MSI. The mass spectrometry imaging of desacetylcinobufagin, resibufogenin, bufalin, arenobufagin, and cinobufagin in the three frozen sections in positive ion mode. A represents the left sagittal sections (sagittal section A), B represents the mid-sagittal sections (sagittal section B), and C represents the horizontal sections (horizontal section C). The corresponding information is summarized in Table 1

Table 1 Information of 5 mainly BDs identified in DESI-MSI