Skip to main content
Download PDF

Lipidomics combined with transcriptomic and mass spectrometry imaging analysis of the Asiatic toad (Bufo gargarizans) during metamorphosis and bufadienolide accumulation

Chinese Medicine volume 17, Article number: 123 (2022) Cite this article

Abstract

Background

To adapt to life on land, Asiatic toads (Bufo gargarizans) must remodel their bodies and refine their chemical defenses in water. The full scope of the mechanisms underlying these processes has yet to be revealed. Bufadienolides (BDs) are chemical defense substances secreted by toads when they are in danger, and they have high medicinal value in treating heart failure, cancer, and hepatitis. However, the artificial breeding of toads to increase BDs has been unsuccessful due to the high mortality of toad larvae during metamorphosis.

Method

Toad larvae at different growth stages were selected to study the changes in the metamorphosis process under the same growth conditions. The differences of tadpoles were explored, including body remodeling, energy metabolism, synthesis and regulation of BDs, through lipidomic technology, transcriptomic technology, and mass spectrometry imaging technology during metamorphosis.

Results

During metamorphosis, tadpoles underwent significant changes in lipid metabolism due to body remodeling to adapt to terrestrial life, which involved ketosis, lipogenesis, cholesterol metabolism, and fatty acid oxidation. The accumulation trend of BDs was observed. “Pentose phosphate pathway” and “Aromatase activity” may be the critical pathway and GO term in BD synthesis, involving 16 genes predominantly expressed in the liver. The involved genes were mainly expressed in the liver, consistent with the synthetic site observed by mass spectrometry imaging.

Conclusion

Together, our findings presented the changes in the toad larvae during metamorphosis and highlighted the accumulation process of BDs as well as the regulatory pathways and synthetic site, providing research and theoretical basis for future development of the toad resources.

Introduction

Toads are one of the most widely distributed amphibians. There are over 300 species in the world, of which the Asiatic toad (Bufo gargarizans) is an ecological generalist, widely distributed in East Asia, and has extremely high medicinal value [1]. Bufadienolides (BDs) enriched in the parotid venom glands of the toads have exhibited strong anti-tumor and anti-inflammatory activities [2, 3]. Chansu (unrefined BDs) is a traditional Chinese animal-derived Chinese medicine in China. However, with the continuous advancement of urbanization in various places, and the increasing environmental degradation, the number of these toads has gradually decreased [4, 5]. To protect species diversity, the Chinese government included B. gargarizans on the protected animals list called “Three Animals.” Chansu also was included on the list of prohibited export commodities, which reflects the scarcity of toads and BDs. As early as the 1960s, the chemical synthesis of BDs had been attempted, but it failed to be applied to industrial production [6, 7].

Obviously, artificial breeding and biosynthesis are shortcuts that will solve the current contradiction between resource conservation and utilization. The implementation is easy for obtaining high-purity medicinal bufadienolides. However, there are considerable obstacles in this solution pathway because of the high mortality of toads during the breeding process. The toads are susceptible to external interference during metamorphosis [8,9,10]. In addition, the regulatory genes governing BD biosynthesis remain unclear, and the enzymatic metabolic pathways and transformational relationships between BDs and synthetic precursors cannot be elucidated.

BDs are important chemical defense substances produced during the growth process of toads and already exist in toad eggs [11, 12]. Parotid venom glands for secreting and storing BDs begin to develop after landing. It is possible that the BDs produced will gradually become more abundant with the growth and development of B. gargarizans. However, the body fat of toads will be gradually decreasing. On the one hand, the tadpoles eat less during the metamorphosis and consume a lot of lipids to maintain energy supply [13]; on the other hand, cholesterol is the synthetic precursor of BDs, and the decrease in lipids may also be related to this. Therefore, conducting a correlation study between lipidomics and transcriptomics of toad larvae at different growth stages is very important. Thus far, there have been few omics studies on the growth and development of B. gargarizans, mainly because the basic information of toads collected in the wild was unclear, such as age, living environment, and genetic background. To ensure high conservation of the samples, we started hatching from the egg belt and tested tadpoles or toads at different stages. Thus, the results will reflect the differences in BDs, lipids, and gene expression at different growth stages.

Transcriptomics and lipidomics are potent tools for biological research because they can accurately reflect the changes in the organism and reveal the relationship between lipids and physiological processes from a unique perspective [14,15,16,17,18,19]. Xu et al. revealed the significant changes in the composition of glycerolipids, glyphospholipids, and sphingolipids in the response of mammalian white adipose tissue to short-term cold exposure, as well as the potential gene regulation process through the conjoint analysis of lipidome and transcriptome [20]. In addition, mass spectrometry imaging technology is a new molecular method that can directly obtain information such as molecular mass, relative abundance, and tissue distribution of small molecular substances in animals from the surface of tissue slices and display them in the form of color images. It is a method widely used to study the distribution of drugs and their metabolites in the body [21,22,23].

The following research was carried out to provide a scientific and technological basis for artificial breeding and biosynthesis and clarify the changes in the body and BDs during metamorphosis. The pre-metamorphosis (Gosner stage 31 and 38), metamorphosis climax (Gosner stage 42), and metamorphosis complete (Gosner stage 46) of B. gargarizans were selected from the same egg mass for examination. Ultra-high performance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS) was used to conduct comprehensive lipidomic analysis, identification, and semi-quantitative analysis of differential lipids and BDs, to reveal the lipid changes in body composition and accumulation of BDs. Combined with transcriptomic analysis, the key regulatory pathways and genes associated with lipids and BDs during metamorphosis were analyzed. Furthermore, desorption electrospray ionization mass spectrometry imaging (DESI-MSI) was used to study the synthesis position of BDs, and the important role of the liver in BD biosynthesis was determined. The results of this study will provide a theoretical basis for metamorphosis research, BD biosynthesis, and artificial breeding of B. gargarizans.

Materials and methods

Animal materials

To reduce the influence of age and environment on sample quality, tadpole eggs from one egg mass were collected in Sikai Town, Liangshan Yi Autonomous Prefecture, Sichuan Province, China. Tadpole eggs were hatched by Xichang Fuhua Jinchan Bio-Tech Co., Ltd. (Sichuan, China). After hatching, these tadpoles were placed in a 4.0 m by 10.0 m by 0.7 m outdoor rearing pond. The temperature of the outdoor rearing tank was maintained between 13 and 17 °C, and the humidity was greater than 70%. The developmental stages of tadpoles were identified according to Gosner, and four stages in metamorphosis were detected. Namely, G31 (Gosner stage 31, pre-metamorphosis), G38 (Gosner stage 38, pre-metamorphosis), G42 (Gosner stage 42, metamorphosis climax), and G46 (Gosner stage 46), metamorphosis complete), as shown in Fig. 1. All of the tadpoles were grouped according to their stage and euthanized by flash freezing in dry ice after fasting for 12 h.

Fig. 1
figure 1

The four life stages of Bufo gargarizans

Full size image

Sample processing

With four tadpoles as one sample, three biological replicates were applied to the later transcriptomic analysis and six biological replicates for lipidomic analysis. The samples were ground into powder using liquid nitrogen and mixed evenly for subsequent analysis. Precisely, 100 mg sample was weighed and put into a centrifugal tube containing 1 mL methanol. After that, the filled centrifugal tubes were vortexed for 30 s and sonicated for 30 min in an ice-cold sonication bath. The suspension was clarified by centrifugation at 12,000 rpm for 10 min at 4 ℃. The 700 µL of supernatant was moved to a clean centrifuge tube, dried under a nitrogen blow dryer, and then dissolved in 150 µL methanol. Exactly 10 μL of each treated sample was taken and mixed into the QC sample for testing on the machine. All reagents used were chromatographically pure.

Chromatographic conditions

Chromatographic separation was performed on a Thermo Ultimate 3000 hyperbaric liquid chromatography (LC) system. The chromatography column used was an Accucore C18 (2.6 μM, 2.1 × 150 mm) purchased from ThermoFisher Scientific. Positive/negative ion mode was used. The mobile phases were 0.1% formic acid in water (A) and acetonitrile (B). The gradient program used was as follows: 0–1 min, 25–40% B; 1–2 min, 25–40% B; 2–12 min, 40–60% B; and 12–22 min, 60–100% B. The flow rate was 0.4 mL/min, injection volume was 10 µL, and the column temperature was 35 °C. The sample was placed in a 4 °C automatic sampler for analysis.

Mass spectrometry conditions

Primary and secondary mass spectrometric data were collected under software (Xcalibur V3.0, ThermoFisher Scientific, San Jose, CA, USA)-controlled Base Peak mode using a hybrid linear ion trap Orbitrap Velos mass spectrometer (LTQ-Orbitrap Velos, Thermo). The ESI source parameters were as follows: ion spray voltage, 3.4 kV; sheath gas flow rate, 35 arb; aux gas flow rate, 10 arb; capillary temperature, 350 °C; and S-lens RF level, 60%. The mass resolution of Fourier transform (FT) was 30,000 with a full scan in the range of m/z 100–1200. The MS/MS experiments were set as data-dependent scans. Data processing was performed using Xcalibur 3.0 and Chromeleon Xpress (version 7.2; ThermoFisher Scientific, CA, USA).

Lipidomic data processing and statistical analysis

For non-targeted lipidomics, all the data obtained were processed using Progenesis QI software for imputing raw data, aligning and selecting peaks, and normalization to produce peak intensities for retention time (tR) and m/z data pairs. The resulting data were exported to SIMCA software (version 14.1.0, Sweden) for principal component analysis (PCA) and orthogonal partial least-squares discrimination analysis (OPLS-DA). Identification was performed using Progenesis QI (version 2.4, Waters, USA) software’s online Lipid Maps Structure Database (LMSD) database, Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and Human Metabolome Database (HMDB) database, and theoretical fragment recognition was conducted. The precursor ion mass number deviation was 10 ppm in each case, and the fragment ion mass number deviation was 5 ppm in each case. Relevant statistical analysis of BDs was performed using GraphPad Prism 8.0 (GraphPad Software, USA).

Sequencing of mRNA and data processing

RNA sequencing (RNA-seq) analysis was outsourced to OE Biotech (Shanghai, China). Briefly, total RNA was extracted using the mirVana miRNA Isolation Kit (Ambion, USA) following the manufacturer’s protocol. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, USA). The untreated samples with RNA Integrity Number (RIN) ≥ 7 were subjected to the subsequent analysis. Libraries were constructed using a TruSeq Stranded mRNA LTSample Prep Kit (Illumina, USA) according to the manufacturer’s instructions. Then, these libraries were sequenced on the Illumina NovaSeq6000, and 125 bp/150 bp paired-end sequences were generated.

The clean sequences were mapped to the B. gargarizans genome (ASM1485885v1) using HISAT2. FPKM of each gene was calculated using Cufflinks, and the sequence counts for each gene were obtained by HTSeqcount. Differential expression analysis was performed using the DESeq2 (2012) R package. P-value < 0.05 and fold-change > 2 or fold-change < 0.5 was set as the threshold for differential expression genes (DEGs). GO enrichment (http://www.geneontology.org/) and KEGG pathway analysis (https://www.genome.jp/kegg/) of DEGs were performed to demonstrate the significantly enriched GO terms and metabolic pathways of DEGs in different groups and samples. All of the charts (bar charts, boxplots, Venn diagrams, bubble charts, and heat maps) were drawn using Oebiotech tools, a free online platform for data analysis (https://cloud.oebiotech.cn/task).

RT-qPCR

RNA was extracted from tadpoles of the four stages and used for RT-qPCR. The primer sequences were designed in the laboratory and synthesized by TsingKe Biotech based on the mRNA sequences obtained from the NCBI database (Additional file 1: Table S1). Three biological repeats and three technical repeats were performed. Quantification was performed by a two-step reaction process: RT (reverse transcription) and PCR. Each RT reaction consisted of 0.5 μg RNA, 0.5 μL of gDNA Remover, and 2 μL of 5 × TransScript All-in-One SuperMix for qPCR, in a total volume of 10 μL. Reactions were performed using the GeneAmp® PCR System 9700 (Applied Biosystems, USA) for 15 min at 42 °C, and 5 s at 85 °C. A 10 × dilution in nuclease-free water was performed for the 10 μL RT reaction mix, which was then stored at − 20 °C.

Real-time PCR was performed using a LightCycler® 480 II Real-time PCR Instrument (Roche, Swiss) with 10 μL PCR reaction mixture that included 1 μL of cDNA, 0.2 μL of forward primer, 5 μL of 2 × PerfectStart™ Green qPCR SuperMix, 3.6 μL of nuclease-free water, and 0.2 μL of reverse primer. Reactions were incubated in a 384-well optical plate (Roche, Swiss) at 94 °C for 30 s, followed by 45 cycles of 94 °C for 5 s, and 60 °C for 30 s. The expression levels of mRNAs were normalized to ACTB and were calculated using the 2−ΔΔCt method [24].

DESI-MSI

The toads at the G46 stage were frozen to death with dry ice, and then their limbs were excised and embedded in 10% gelatin until they could be demolded. Whole-body cryosections were obtained using a Leica CM 3500 (Leica Microsystems, Nussloch, Germany) microtome at 100-μm increments. Each section was adhered to a positively charged coated slide (Citotest, Haimen, China). Serial tissue sections for optical microscopy and DESI-MSI were alternately collected to allow correlation and organ identification.

Mass spectrometric data were collected under software (HDImaging V1.4, Waters, Milford, MA, USA)-controlled Base Peak mode using a mass analyzer Xevo G2-XS QToF mass spectrometer (Waters, Milford, MA, USA). The DESI source parameters were as follows: collision voltage, off; spatial resolution, 200 μm × 200 μm; scan rate, 500 μm/s; spray solvent, 95% methanol in water with 200 pg/μL LE; spraying rate, 5 μL/min; gas nebulizer (N2) pressure, 0.5 Mpa; spray angle, 60°. The mass analysis mode of the Xevo G2 XS GToF was sensitivity mode with a full scan in the range of m/z 300–500. The ion mode was positive, and the instrument was tuned to a signal strength above 1e6 using a red Sharpie pen (rhodamine 6G, m/z 443). Data processing was performed using HDImaging 1.4. To observe the parotid gland and additional tissues and organs, three regions of the body in the left sagittal, mid-parasagittal, and horizontal planes were analyzed [25, 26]

Results

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
figure 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]

Full size image

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
figure 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

Full size image

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
figure 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)

Full size image

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
figure 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

Full size image

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
figure 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

Full size image
Table 1 Information of 5 mainly BDs identified in DESI-MSI
Full size table

Discussion

To adapt to life on land, toads must complete morphological changes and remodeling of organs and cell types in water [40]. In this study, the changes in the lipidome and transcriptome during toad development and metamorphosis, as well as the accumulation characteristics of BDs, were reported. The results showed that the lipid and gene expression in larvae were markedly up- or downregulated in response to body changes during metamorphosis, and BDs continued to increase from G38 and peaked at G46. Additionally, the distribution of BDs in juveniles was observed using whole-body cryosections and DESI-MSI. The results indicated that the liver plays a vital role in the biosynthesis and accumulation of BDs.

Lipid and transcriptional changes in B. gargarizans during metamorphosis

89 differential lipids were identified in this process, and the total relative content of differential lipids gradually decreased (Additional file 8: Fig. S1). The most apparent decrease in differential lipids was in sterol lipids (ST), which mainly include a large quantity 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 enrichment analysis results for DEGs (G31 vs. G38, G38 vs. G42) also suggested a decrease in bile secretion and tyrosine metabolism, indicating that loss of appetite and decreased food intake occured during metamorphosis (Additional file 11: Fig. S4B and D).

Before metamorphosis, tadpoles actively ate to provide adequate nutritional preparation and induce the process of metamorphosis [9, 41]. After metamorphosis began, the hormones leptin and stress neuropeptide corticotropin-releasing factor acted on feeding control centers in the hypothalamus to inhibit feeding [42]. Decreased eating not only assisted with intestinal remodeling, but also reduced the body fat of toads, providing favorable conditions for adapting to terrestrial life. However, large amounts of energy were still required for the morphological changes during metamorphosis. The tail of the tadpole began to be absorbed at the G38 stage, and the relative quantities of fatty acids (FAs) were also decreasing, indicating that the tail changed from an energy-consuming organ to an energy-supplying organ, replacing the fat depot as an energy source for metamorphosis [9]. However, the low appetite phase was likely to last longer than the tail absorption period. Amino acids or proteins that flowed from the tail to the liver were adequately mobilized as extrahepatic nutrients due to the inconvenience of bulk storage. KEGG enrichment analysis showed that the gluconeogenesis function of tadpoles in the metamorphosis climax stage was continuously enhanced (Additional file 11: Fig. S4C and E). Enhanced gluconeogenesis provided energy for the low appetite phase following complete absorption of the tail, and the increased CAR and decreased FA have suggested that starvation-induced enhanced lipid catabolism in the liver and muscle [43].

Clinically, PG and PI are auxiliary indicators for the diagnosis of fetal lung maturity. In Fig. 3B, the opposite trend was observed for them. The ratio of PG to PI in the amniotic fluid increases gradually as the fetal lung develops to maturity. In this study, to facilitate respiration in air, the primary larval gills had transformed into lung breathing through morphological innovations, which was also a process of lung maturation. The relative quantities of lyso-phosphatidyl cholines (LPC) decreased at G46. The main component, LPC 20:4(5Z,8Z,11Z,14Z), is derived from arachidonic acid, and its level is linearly related to the number of preserved beta cells [44]. The consumption of LPC 20:4(5Z,8Z,11Z,14Z) from G42 to G46 may be caused by the remodeling of the endocrine pancreas, while the results of the KEGG enrichment analysis showed that pancreatic secretion was reduced since the G31 stage, indicating that beta cell death may be simply the final process of pancreatic remodeling. In the adult Xenopus pancreas, islets contain only insulin-producing beta cells, and it has been proposed that the maturation of the pancreatic system occurs at the metamorphosis climax [41]. However, beta cells do not replicate and are involved in limited cell death and changes in insulin expression during this process. LPC22:6 confers protective effects on the heart and brain, and supplementing it in cardiac arrest rats can aid in cardiac resuscitation. Thus, the high consumption after G31 may promote cardiac development. The KEGG enrichment analysis also showed that cardiac muscle contraction was continuously enhanced during metamorphosis (Additional file 12: Fig. S5B) [45].

Pinelli et al. performed a comparative analysis of apoptotic and mitotic indices in the developing eye of toad larvae, and they found that the peaks of apoptosis and proliferation in retinal cells appeared in G31 and G38 [46], which was consistent with our enrichment analysis results (Additional file 11: Fig. S4A and C). PE(P-20:0/0:0), which is a marker of apoptotic cells [47, 48], was highly expressed when the organ began to remodel, gradually decreased with redevelopment, and finally increased again when the tail was apoptotic (Fig. 3B).

Mechanism of regulating differential lipids and BDs

The metabolism and composition of lipids were comprehensively analyzed in 4 different developmental stages for tadpoles, and it was found that the genes relevant to differential lipids were mainly enriched in the PPAR signaling pathway. There are three different isoforms of peroxisome proliferator-activated receptors (PPARs) (i.e., PPARα, PPARβ/δ, and PPARγ), and they function as ligand-activated transcription factors with DNA-binding domains that recognize response elements in the promoter region of specific target genes linked to inflammation, apoptosis, cell proliferation, and differentiation [49]. The PPAR signaling pathway was important in regulating lipid metabolism during metamorphosis.

In the present study, the sharp decrease in the expression of genes related to fatty acid oxidation after G31 was due to the priority consumption of amino acids and glycogen rather than liver lipids. The expression trend of cholesterol metabolism-related genes supported this point. Starvation will lead to ketone production, and the increase in ketone bodies induced proliferation, differentiation, and apoptosis in intestinal crypts, which accelerated the remodeling process.

Extrahepatic tissues also used ketone bodies to generate large amounts of acetyl CoA, which can inhibit pyruvate dehydrogenase activity, limit glucose consumption, activate pyruvate carboxylase, and promote gluconeogenesis. This can reduce glucose uptake by extrahepatic tissues to ensure that sufficient glucose reaches the brain tissue and red blood cells. In addition, brain tissues cannot use long-chain fatty acids, but ketone bodies can be used for energy [50]. The expression trend of related genes of ketogenesis (HMGCS2) and gluconeogenesis (UBC) also reflected that. The continuous decrease in SCD-1, an essential gene of adipogenesis, inhibited adipogenesis so that the limbs could better control the body after landing.

The PPAR signaling pathway identified by combined lipidomic and transcriptomic analysis reveals the regulatory process of lipids and genes during metamorphosis, which is an excellent method for analyzing the biosynthetic mechanism of BDs. According to the accumulation trend of BDs at the four stages, 3260 DEGs were identified, and 689 DEGs strongly associated with BDs were identified by Pearson correlation analysis. KEGG and GO enrichment analyses were subsequently performed on the strongly related DEGs. The enrichment results prompted two important synthetic pathways for us. Two important synthetic pathways, “pentose phosphate pathway” and “aromatase activity,” were led according to the enrichment results. The primary physiological function of the pentose phosphate pathway (PPP) is the production of NADPH, which provides reducing power for biosynthesis. It is likely that the NADPH produced by PPP assists in converting the parent nucleus of cholesterol into the cis-fused ring of BDs. In addition, “aromatase activity” was related to catalyzing the reduction of aliphatic rings to generate aromatic rings. The DEGs in “aromatase activity” may be critical regulatory genes for the cyclization of cholesterol side chains into rings.

The pathways involved 16 DEGs, of which 10 were derived from cytochrome P450s, including CYP3A29, CYP2K1, CYP2B4, CYP2D15, CYP2K4, CYP2F3, CYP2C54, CYP2A10, CYP2G1, CYP2H2. Cytochrome P450s are monooxygenases that catalyze reactions involving drug metabolism and synthesis of lipids (cholesterol, steroids, etc.). The regulations of cytochrome P450s were required for both the synthesis and accumulation of BDs and the self-protection caused by endogenous BDs. The cytochrome P450s in DEGs were almost all from the CYP2 family. The annotation information and literature indicate that these CYP2s have the following functions in synthesising BDs. Firstly, CYP2F3 has a clear role in dehydrogenation [51]. It may be involved in the dehydrogenation of 3-OH, which is a key step in the isomerization of the steroidal nucleus of cholesterol to cis-fused rings. Secondly, the conversion of cholesterol to BDs also involves hydroxylation at the C-1β, 3β, 5β, 11α, 12α/β, 14β, and 16β sites. CYP2K1, CYP2D15, CYP2A10 and CYP2H2 may be related to that, which can hydroxylate steroids or lipids [52,53,54,55]. Lastly, CYP2B4, CYP2C54, CYP2G1, and CYP2H2 were involved in the cyclooxygenase P450 pathway [56], which may play a key role in the generation of unsaturated six-membered lactone ring from bufadienolide. It can be inferred that CYP2s are very important in the biosynthesis of BDs.

Synthetic sites and distribution of BDs

DESI-MSI combined with whole-body cryosectioning was an effective method to study the tissue distribution and synthesis sites of BDs. This is the first time this method has been applied to studying BD tissue distribution and biosynthesis. The whole-body mapping of BDs in toads by DESI-MSI was the initial step to gaining insight into BD biosynthesis. In Fig. 6A, BDs are characteristically distributed in the liver but not enriched in the immature parotid gland. BDs were most abundant in the parotid glands of adult toads. By stimulating these glands, they can secrete large quantities of Chansu (the crude extract of BDs), which can be used as the raw material for clinical preparations.

Excessive storage of BDs, which may result in excessive peak height in the gland and hinder the observation of other organs and tissues, will not occur when juvenile toad larvae are used as the samples for tissue distribution research. The results in Fig. 6A suggest that BDs can be synthesized in the liver, and the parotid gland may simply be a storage organ for BDs. Garraffo et al. demonstrated that the glands were unable to produce the synthetic precursor cholesterol of BDs when sodium (1-14C) acetate and (5-3H) mevalonate were co-incubated with toad liver tissue and postauricular gland tissue, respectively [57], which also supported our opinion.

Fig. 6B shows that BDs also had a distinct characteristic distribution in the gallbladder, and the amount of cinobufagin in the gallbladder was much higher than that in the liver. This suggests that the gallbladder may be the first collection organ of BDs, and it is responsible for transporting BDs produced by the liver to the skin and glands, while the BDs of juvenile toads mainly accumulated in the liver and gallbladder. The analysis of horizontal sections also showed that a large number of BDs were enriched in the liver (Fig. 6C), again confirming the importance of the liver in the synthesis of BDs and implicating it as an independent synthesis site of BDs.

Conclusion

In summary, lipidomics and transcriptomic methods were used to study the metamorphic development of bufo gargarizans, demonstrating the diet, energy metabolism, and growth of the organs of toads. The PPAR signaling pathway was mapped as the main lipids regulatory pathway, demonstrated the processes of G31-G46 regulating ketogenesis, adipogenesis, cholesterol metabolism, and fatty acid oxidation. The accumulation trends of seven BDs were identified during the metamorphosis with standards and suggested the key pathways of BDs biosynthesis combined with the transcriptomic data. Furthermore, the study of BDs tissue distribution by DESI-MSI indicated that the liver plays a crucial role in BDs biosynthesis and may be the synthetic site of BDs. Overall, this study provides a valuable data resource for the resource development of B. gargarizans and BDs.

Availability of data and materials

The RNA-seq reads of the toads at the four stages can be obtained from BioProject PRJNA853591. The assembled Bufo gargarizans genome ASM1485885v1 was deposited in NCBI under BioProject PRJNA628553.

Abbreviations

DESI:

Desorption electrospray ionization

LC–MS:

Liquid chromatography-tandem mass spectrometry

DEGs:

Differentially expressed genes

BDs:

Bufadienolides

MSI:

Mass spectrometry imaging

PPAR:

Peroxisome proliferators-activated receptor

G:

Gosner stage

KEGG:

Kyoto Encyclopedia of Genes and Genomes

GO:

Gene ontology

FA:

Fatty acids

ST:

Sterol

PI:

Phosphatidylinositol

LPC:

Lyso-posphatidyl choline

PE:

Phosphatidylethanolamine

CAR:

Carnitines

PG:

Phosphatidylglycerol

PCA:

Orthogonal partial least-squares discrimination analysis

OPLS-DA:

Orthogonal partial least-squares discrimination analysis

RT-qPCR:

Quantitative real-time PCR

GC:

Guanine and cytosine

PPP:

Pentose phosphate pathway

References

  1. Lu B, Jiang J, Wu H, Chen X, Song X, Liao W, Fu J. A large genome with chromosome-scale assembly sheds light on the evolutionary success of a true toad (Bufo gargarizans). Mol Ecol Resour. 2021;21(4):1256–73.

    Article  CAS  PubMed  Google Scholar 

  2. Kaplan W, Khatri HR, Nagorny P. Concise enantioselective total synthesis of cardiotonic steroids 19-hydroxysarmentogenin and Trewianin aglycone. J Am Chem Soc. 2016;138(22):7194–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Li Y, Liu X, Lin X, Zhao M, Xiao Y, Liu C, Liang Z, Lin Z, Yi R, Tang Z, et al. Chemical compound cinobufotalin potently induces FOXO1-stimulated cisplatin sensitivity by antagonizing its binding partner MYH9. Signal Transduct Target Therapy. 2019;4:48.

    Article  CAS  Google Scholar 

  4. Bokony V, Uveges B, Verebelyi V, Ujhegyi N, Moricz AM. Toads phenotypically adjust their chemical defences to anthropogenic habitat change. Sci Rep. 2019;9(1):3163.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Kaczmarski M, Kolenda K, Rozenblut-Koscisty B, Sosnicka W. Phalangeal bone anomalies in the European common toad Bufo bufo from polluted environments. Environ Sci Pollut Res Int. 2016;23(21):21940–6.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Sondheimer F, McCrae W, Salmond WG. Synthesis of bufadienolides. Synthesis of bufalin and resibufogenin. J Am Chem Soc. 1969;3(3):1228–30.

    Article  Google Scholar 

  7. Yoshii E, Koizumi T, Ikeshima H, Ozaki K, Hayashi I. Studies on the synthesis of cardiotonic steroids. I. Efficient synthesis of cardenolides. Chem Pharm Bull (Tokyo). 1975;23(11):2496–506.

    Article  CAS  Google Scholar 

  8. Yin X, Jiang S, Yu J, Zhu G, Wu H, Mao C. Effects of spirotetramat on the acute toxicity, oxidative stress, and lipid peroxidation in Chinese toad (Bufo bufo gargarizans) tadpoles. Environ Toxicol Pharmacol. 2014;37(3):1229–35.

    Article  CAS  PubMed  Google Scholar 

  9. Zhu W, Chang L, Zhao T, Wang B, Jiang J. Remarkable metabolic reorganization and altered metabolic requirements in frog metamorphic climax. Front Zool. 2020;17:30.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Ishizuya-Oka A, Hasebe T, Shi YB. Apoptosis in amphibian organs during metamorphosis. Apoptosis. 2010;15(3):350–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Zhou S, Zheng Q, Huang X, Wang Y, Luo S, Jiang R, Wang L, Ye W, Tian H. Isolation and identification of l/d-lactate-conjugated bufadienolides from toad eggs revealing lactate racemization in amphibians. Org Biomol Chem. 2017;15(26):5609–15.

    Article  CAS  PubMed  Google Scholar 

  12. Zhou SW, Quan JY, Li ZW, Ye G, Shang Z, Chen ZP, Wang L, Li XY, Zhang XQ, Li J, et al. Bufadienolides from the eggs of the toad Bufo bufo gargarizans and their antimelanoma activities. J Nat Prod. 2021;84(5):1425–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Wright ML, Richardson SE, Bigos JM. The fat body of bullfrog (Lithobates catesbeianus) tadpoles during metamorphosis: changes in mass, histology, and melatonin content and effect of food deprivation. Comp Biochem Physiol A Mol Integr Physiol. 2011;160(4):498–503.

    Article  CAS  PubMed  Google Scholar 

  14. Zullig T, Kofeler HC. High resolution mass spectrometry in lipidomics. Mass Spectrom Rev. 2021;40(3):162–76.

    Article  PubMed  Google Scholar 

  15. Zhang L, Zhu B, Zeng Y, Shen H, Zhang J, Wang X. Clinical lipidomics in understanding of lung cancer: opportunity and challenge. Cancer Lett. 2020;470:75–83.

    Article  CAS  PubMed  Google Scholar 

  16. Peng B, Geue S, Coman C, Munzer P, Kopczynski D, Has C, Hoffmann N, Manke MC, Lang F, Sickmann A, et al. Identification of key lipids critical for platelet activation by comprehensive analysis of the platelet lipidome. Blood. 2018;132(5):e1–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Rosenberg S, Verreault M, Schmitt C, Guegan J, Guehennec J, Levasseur C, Marie Y, Bielle F, Mokhtari K, Hoang-Xuan K, et al. Multi-omics analysis of primary glioblastoma cell lines shows recapitulation of pivotal molecular features of parental tumors. Neuro Oncol. 2017;19(2):219–28.

    CAS  PubMed  Google Scholar 

  18. Lu J, Chen B, Chen T, Guo S, Xue X, Chen Q, Zhao M, Xia L, Zhu Z, Zheng L, et al. Comprehensive metabolomics identified lipid peroxidation as a prominent feature in human plasma of patients with coronary heart diseases. Redox Biol. 2017;12:899–907.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Biagiotti M, Dominguez S, Yamout N, Zufferey R. Lipidomics and anti-trypanosomatid chemotherapy. Clin Transl Med. 2017;6(1):27.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Xu Z, You W, Zhou Y, Chen W, Wang Y, Shan T. Cold-induced lipid dynamics and transcriptional programs in white adipose tissue. BMC Biol. 2019;17(1):74.

    Article  PubMed  PubMed Central  Google Scholar 

  21. He J, Luo Z, Huang L, He J, Chen Y, Rong X, Jia S, Tang F, Wang X, Zhang R, et al. Ambient mass spectrometry imaging metabolomics method provides novel insights into the action mechanism of drug candidates. Anal Chem. 2015;87(10):5372–9.

    Article  CAS  PubMed  Google Scholar 

  22. Santoro AL, Drummond RD, Silva IT, Ferreira SS, Juliano L, Vendramini PH, Lemos M, Eberlin MN, Andrade VP. In situ DESI-MSI lipidomic profiles of breast cancer molecular subtypes and precursor lesions. Cancer Res. 2020;80(6):1246–57.

    Article  CAS  PubMed  Google Scholar 

  23. Hou J, Zhang Z, Zhang L, Wu W, Huang Y, Jia Z, Zhou L, Gao L, Long H, Lei M, et al. Spatial lipidomics of eight edible nuts by desorption electrospray ionization with ion mobility mass spectrometry imaging. Food Chem. 2022;371: 130893.

    Article  CAS  PubMed  Google Scholar 

  24. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25(4):402–8.

    Article  CAS  PubMed  Google Scholar 

  25. Jeckel AM, Matsumura K, Nishikawa K, Morimoto Y, Saporito RA, Grant T, Ifa DR. Use of whole-body cryosectioning and desorption electrospray ionization mass spectrometry imaging to visualize alkaloid distribution in poison frogs. J Mass Spectrom. 2020;55(6): e4520.

    Article  CAS  PubMed  Google Scholar 

  26. Nelson KA, Daniels GJ, Fournie JW, Hemmer MJ. Optimization of whole-body zebrafish sectioning methods for mass spectrometry imaging. J Biomol Tech. 2013;24(3):119–27.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Wigger L, Barovic M, Brunner AD, Marzetta F, Schoniger E, Mehl F, Kipke N, Friedland D, Burdet F, Kessler C, et al. Multi-omics profiling of living human pancreatic islet donors reveals heterogeneous beta cell trajectories towards type 2 diabetes. Nat Metab. 2021;3(7):1017–31.

    Article  CAS  PubMed  Google Scholar 

  28. Wang J, Chen X, Chu S, You Y, Chi Y, Wang R, Yang X, Hayat K, Zhang D, Zhou P. Comparative cytology combined with transcriptomic and metabolomic analyses of Solanum nigrum L. in response to Cd toxicity. J Hazard Mater. 2022;423(Pt B): 127168.

    Article  CAS  PubMed  Google Scholar 

  29. Kammerer L, Mohammad G, Wolna M, Robbins PA, Lakhal-Littleton S. Fetal liver hepcidin secures iron stores in utero. Blood. 2020;136(13):1549–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Qin Y, Sun H, Hao P, Wang H, Wang C, Ma L, Wei H, Yu S. Transcriptome analysis reveals differences in the mechanisms of fiber initiation and elongation between long- and short-fiber cotton (Gossypium hirsutum L.) lines. BMC Genomics. 2019;20(1):633.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Kunz C, Hower J, Knoll A, Ritzenthaler KL, Lamberti T. No improvement in vitamin D status in German infants and adolescents between 2009 and 2014 despite public recommendations to increase vitamin D intake in 2012. Eur J Nutr. 2019;58(4):1711–22.

    Article  CAS  PubMed  Google Scholar 

  32. Brehm JM, Celedon JC, Soto-Quiros ME, Avila L, Hunninghake GM, Forno E, Laskey D, Sylvia JS, Hollis BW, Weiss ST, et al. Serum vitamin D levels and markers of severity of childhood asthma in Costa Rica. Am J Respir Crit Care Med. 2009;179(9):765–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Newgard CB. Metabolomics and metabolic diseases: where do we stand? Cell Metab. 2017;25(1):43–56.

    Article  CAS  PubMed  Google Scholar 

  34. Xu Y, Li F, Zalzala M, Xu J, Gonzalez FJ, Adorini L, Lee YK, Yin L, Zhang Y. Farnesoid X receptor activation increases reverse cholesterol transport by modulating bile acid composition and cholesterol absorption in mice. Hepatology. 2016;64(4):1072–85.

    Article  CAS  PubMed  Google Scholar 

  35. Zurkinden L, Solca C, Vogeli IA, Vogt B, Ackermann D, Erickson SK, Frey FJ, Sviridov D, Escher G. Effect of Cyp27A1 gene dosage on atherosclerosis development in ApoE-knockout mice. FASEB J. 2014;28(3):1198–209.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Antica M, Kusic B, Hranilovic D, Dietz AB, Vuk-Pavlovic S. Cloning the cDNA for murine U2 snRNP-A’ gene and its differential expression in lymphocyte development. Immunol Lett. 2002;82(3):217–23.

    Article  CAS  PubMed  Google Scholar 

  37. Dickey AS, Sanchez DN, Arreola M, Sampat KR, Fan W, Arbez N, Akimov S, Van Kanegan MJ, Ohnishi K, Gilmore-Hall SK, et al. PPARdelta activation by bexarotene promotes neuroprotection by restoring bioenergetic and quality control homeostasis. Sci Transl Med. 2017. https://doi.org/10.1126/scitranslmed.aal2332.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Sharma R, Ramanathan A. The aging metabolome-biomarkers to hub metabolites. Proteomics. 2020;20(5–6): e1800407.

    Article  PubMed  Google Scholar 

  39. Boric M, Danevcic T, Stopar D. Viscosity dictates metabolic activity of Vibrio ruber. Front Microbiol. 2012;3:255.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Inoue T, Nakata R, Savitzky AH, Yoshinaga N, Mori A, Mori N. Variation in bufadienolide composition of parotoid gland secretion from three taxa of Japanese toads. J Chem Ecol. 2020;46(10):997–1009.

    Article  CAS  PubMed  Google Scholar 

  41. Pearl EJ, Bilogan CK, Mukhi S, Brown DD, Horb ME. Xenopus pancreas development. Dev Dyn. 2009;238(6):1271–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Chai L, Dong Z, Chen A, Wang H. Changes in intestinal microbiota of Bufo gargarizans and its association with body weight during metamorphosis. Arch Microbiol. 2018;200(7):1087–99.

    Article  CAS  PubMed  Google Scholar 

  43. Zhu W, Zhang M, Chang L, Zhu W, Li C, Xie F, Zhang H, Zhao T, Jiang J. Characterizing the composition, metabolism and physiological functions of the fatty liver in Rana omeimontis tadpoles. Front Zool. 2019;16:42.

    Article  PubMed  PubMed Central  Google Scholar 

  44. de Groot P, Nikolic T, Pellegrini S, Sordi V, Imangaliyev S, Rampanelli E, Hanssen N, Attaye I, Bakker G, Duinkerken G, et al. Faecal microbiota transplantation halts progression of human new-onset type 1 diabetes in a randomised controlled trial. Gut. 2021;70(1):92–105.

    Article  PubMed  Google Scholar 

  45. Nishikimi M, Yagi T, Shoaib M, Takegawa R, Rasul R, Hayashida K, Okuma Y, Yin T, Choudhary RC, Becker LB, et al. Phospholipid screening postcardiac arrest detects decreased plasma lysophosphatidylcholine: supplementation as a new therapeutic approach. Crit Care Med. 2022;50(2):e199–208.

    Article  CAS  PubMed  Google Scholar 

  46. Pinelli C, Sansone A, De Maio A, Morgillo A, Scandurra A, D’Aniello B. Proliferative events and apoptotic remodelling in retinal development of common toad (Bufo bufo). J Exp Zool B Mol Dev Evol. 2016;326(1):19–30.

    Article  CAS  PubMed  Google Scholar 

  47. Xue L, Li M, Chen T, Sun H, Zhu J, Li X, Wu F, Wang B, Li J, Chen Y. PEinduced apoptosis in SMMC7721 cells: involvement of Erk and Stat signalling pathways. Int J Mol Med. 2014;34(1):119–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Sun G, Li X, Wei J, Zhang T, Li B, Chen M, Wang Y, Chen K, Li Y. Pharmacodynamic substances in Salvia miltiorrhiza for prevention and treatment of hyperlipidemia and coronary heart disease based on lipidomics technology and network pharmacology analysis. Biomed Pharmacother. 2021;141: 111846.

    Article  CAS  PubMed  Google Scholar 

  49. Claria J. Regulation of cell proliferation and apoptosis by bioactive lipid mediators. Recent Pat Anticancer Drug Discov. 2006;1(3):369–82.

    Article  CAS  PubMed  Google Scholar 

  50. Xie Z, Zhang D, Chung D, Tang Z, Huang H, Dai L, Qi S, Li J, Colak G, Chen Y, et al. Metabolic regulation of gene expression by histone lysine beta-hydroxybutyrylation. Mol Cell. 2016;62(2):194–206.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Kartha JS, Skordos KW, Sun H, Hall C, Easterwood LM, Reilly CA, Johnson EF, Yost GS. Single mutations change CYP2F3 from a dehydrogenase of 3-methylindole to an oxygenase. Biochemistry. 2008;47(37):9756–70.

    Article  CAS  PubMed  Google Scholar 

  52. Yang YH, Miranda CL, Henderson MC, Wang-Buhler JL, Buhler DR. Heterologous expression of CYP2K1 and identification of the expressed protein (BV-CYP2K1) as lauric acid (omega-1)-hydroxylase and aflatoxin B1 exo-epoxidase. Drug Metab Dispos. 2000;28(11):1279–83.

    CAS  PubMed  Google Scholar 

  53. Shou M, Norcross R, Sandig G, Lu P, Li Y, Lin Y, Mei Q, Rodrigues AD, Rushmore TH. Substrate specificity and kinetic properties of seven heterologously expressed dog cytochromes p450. Drug Metab Dispos. 2003;31(9):1161–9.

    Article  CAS  PubMed  Google Scholar 

  54. Tasaki T, Iwata H, Kazusaka A, Fujita S. Regio- and stereoselectivity in propranolol metabolism by dog liver microsomes and the expressed dog CYP2D15. J Biochem. 1998;123(4):747–51.

    Article  CAS  PubMed  Google Scholar 

  55. Ding X, Peng HM, Coon MJ. Structure-function analysis of CYP2A10 and CYP2A11, P450 cytochromes that differ in only eight amino acids but have strikingly different activities toward testosterone and coumarin. Biochem Biophys Res Commun. 1994;203(1):373–8.

    Article  CAS  PubMed  Google Scholar 

  56. El-Awaad I, Bocola M, Beuerle T, Liu B, Beerhues L. Bifunctional CYP81AA proteins catalyse identical hydroxylations but alternative regioselective phenol couplings in plant xanthone biosynthesis. Nat Commun. 2016;7:11472.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Garraffo HM, Gros EG. Biosynthesis of bufadienolides in toads. VI. Experiments with [1,2–3H]cholesterol, [21–14C]coprostanol, and 5 beta-[21–14 C]pregnanolone in the toad Bufo arenarum. Steroids. 1986;48(3–4):251–7.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank the authors and submitting laboratories of the sequences from the Bufo gargarizans; we thank Yi Shen for her operational assistance during the preparation of whole-body cryosections; we also thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript.

Funding

This work was supported by the CACMS Innovation Fund (CI2021A04012), the Fundamental Research Funds for the Central public welfare research institutes (ZXKT22004), National Administration of Traditional Chinese Medicine Young Qihuang Scholar Project, China; Qianjiang Distinguished Expert Project of Hangzhou, China (2020), and the Toxicology Innovation Team (Z2021180).

Author information

Authors and Affiliations

  1. Institute of Chinese Materia Medica, China Academy of Chinese Medical Sciences, Beijing, 100700, China

    Bo Sun, Shan Jiang, Mingli Li, Yan Zhang, Yanyan Zhou, Xiaolu Wei, Hongjie Wang, Nan Si, Baolin Bian & Haiyu Zhao

Authors
  1. Bo Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Shan Jiang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Mingli Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yan Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yanyan Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Xiaolu Wei
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Hongjie Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Nan Si
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Baolin Bian
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Haiyu Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

N.S., B.-L. B., and H.-Y.Z. devised the project.B.S., S.J. and M.-L.L. prepared materials for lipidome and transcriptome sequencing with assistance from Y.Z., Y.-Y.Z., X.-L.W., and H.-J.W. performed the experiments. B.S., and J.S. analysed and interpreted the results of structural variations with assistance from S.N. B.S., H.-Y.Z., and J.S. wrote and substantively revised the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Nan Si, Baolin Bian or Haiyu Zhao.

Ethics declarations

Ethics approval and consent to participate

All animal experiments were performed following the Guide for the Care and Use of Laboratory Animals. The research protocol was approved by the Ethics Committee at Institute of Chinese Materia Medica, China Academy of Chinese Medical Sciences.

Consent for publication

Not applicable.

Competing of interests

The authors declare that they have no conflict of interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Specific primer pairs for selected genes.

Additional file 2: Table S2.

Differential lipids in four stages.

Additional file 3: Table S3.

BDs semi-quantitative analysis reports.

Additional file 4: Table S4.

After raw data filtering and examination of the sequencing error rate and GC content distribution.

Additional file 5: Table S5.

GO enrichment analysis (Top 20).

Additional file 6: Table S6.

KEGG enrichment analysis (Top 20).

Additional file 7: Table S7.

Differential lipids related DEGs (TOP 500).

Additional file 8: Fig. S1.

Cumulative lipids composition.

Additional file 9: Fig. S2.

Multivariate analysis of the transcriptome in B. gargarizans at different growth stages. (A) The degree of dispersion of the FPKM distribution. The horizontal axis is the sample name, and the vertical axis is for log10(FPKM+1). Box plot elements denote maximum, third quartile, median, first quartile, and minimum. (B) Heatmap for Pearson correlation coefficients (PCC) between samples. The horizontal axis represents the sample name, the vertical axis represents the corresponding sample name, and the color denotes the size of the correlation coefficient. (C) Principal component analysis (PCA) plot of the transcriptomic data from the 12 samples. (D and E) Up- and downregulated differential Venn diagram. Each circle represents a comparison group, the numbers outside the overlaps represent the number of specific DEGs of the comparison group, and the numbers in the circle overlaps represent the number of DEGs common to the comparison group.

Additional file 10: Fig. S3.

GO enrichment bubble diagram (G31 vs. G38, G38 vs. G42, G42 vs. G46). (A) and (B) showed the top 30 pathways of significant enrichment of the upregulated and downregulated DEGs on GO for G31 vs. G38, (C) and (D) showed for G38 vs. G42, and (E) and (F) showed for G42 vs. G46.

Additional file 11: Fig. S4.

KEGG enrichment bubble diagram (G31 vs. G38, G38 vs. G32, G42 vs. G46). (A) and (B) showed the top 30 pathways of significant enrichment of the up-and down- regulated DEGs on KEGG for G31 vs. G38, (C) and (D) showed for G38 vs. G42, and (E) and (F) showed for G42 vs. G46.

Additional file 12: Fig. S5.

Enrichment analysis of DEGs for B. gargarizans metamorphosis. The top 30 GO terms and pathways of significant enrichment for the (A and B) upregulated DEGs and (C and D) downregulated DEGs. The y-axis corresponds to the GO term or KEGG pathway, the x-axis corresponds to the EnrichmentScore, and the size of the point corresponds to the number of differential genes in the term or pathway.

Additional file 13: Fig. S6.

Enrichment analysis of differential lipid-related DEGs. (A) showed the top 30 GO terms of significant enrichment of the differential lipid-related DEGs on GO;(B) showed the top 30 pathways of significant enrichment of the differential lipid-related DEGs on KEGG.

Additional file 14: Fig. S7.

Results of RT-qPCR.

Additional file 15: Fig. S8.

Enrichment analysis of DEGs for BD accumulation. (A) By comparing G46 with G31, G38 and G42 respectively, 3260 common DEGs were identified. (B) and (C) showed the top 30 GO terms and pathways of significant enrichment of the BDs-related DEGs.

Additional file 16: Fig. S9.

General overview of the major steps in DESI-MSI in the entire toad at the G46 stage. Step 1. After removing the limbs of the toads (G46 stage), whole-body cryosections were created along the sagittal and horizontal planes. Step 2. Serial tissue sections were collected while performing the optical microscopy and DESI-MSI techniques, and the appropriate sections were then selected for testing. DESI-MSI was performed, including generation of mass spectra and conversion to image, and identification of organs and tissues with optical microscopy.

Additional file 17: Fig. S10.

Identification of organs in DESI-MSI sections. (A–C) Frozen sections selected for DESI-MSI; (D–F) DESI-MSI at a specific mass-to-charge ratio that presents a clear org an position. (a–c) Microscope images of the corresponding parts of A to C. A, D and a representative of the left sagittal sections (sagittal section A); B, E, and b represent the mid-sagittal sections (sagittal section B); C, F, and c represent the horizontal sections (horizontal section C).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, B., Jiang, S., Li, M. et al. Lipidomics combined with transcriptomic and mass spectrometry imaging analysis of the Asiatic toad (Bufo gargarizans) during metamorphosis and bufadienolide accumulation. Chin Med 17, 123 (2022). https://doi.org/10.1186/s13020-022-00676-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13020-022-00676-7

Keywords