Ontogenetic Changes in the Digestive Capacities of the Naozhou Stock of Large Yellow Croaker (Larimichthys crocea)

TL;DR

This study tracks how young large yellow croaker develop their digestive systems as they transition from relying on yolk to eating live prey and then artificial diets.

Contribution

The study provides new insights into the ontogenetic changes in digestive and metabolic capacities during early development of large yellow croaker.

Findings

Protease activity increases sharply with the introduction of live prey, while amylase and lipase activities rise during transition to formulated feeds.

Transcriptomic and metabolomic analyses reveal coordinated changes in digestion, metabolism, and energy use across developmental stages.

Digestive enzyme activities, gene expression, and metabolite profiles enable a smooth transition from yolk dependency to artificial diets.

Abstract

Young fish are fragile in their early days, and many do not survive because their bodies are still figuring out how to handle different kinds of food. In this study, the young large yellow croaker grows from the stage where they rely on their yolk for nourishment to the point where they can eat tiny live animals and eventually a prepared diet. The study tracked how their digestive system developed and how they used nutrients at each step. The study found out that they start by depending on simple internal food reserves, then quickly strengthen their ability to break down proteins and fats once they begin eating live prey, and finally, they learn to handle more complex diets as their organs mature. The study also finds that changes in digestion were closely tied to shifts in growth, energy use, and organ development. These findings provide new insights into the early nutritional development of NZ large yellow croaker and provide a scientific basis for the improvement of artificial aquaculture seed production.

This study examined the digestive and metabolic responses of Naozhou (NZ) stock large yellow croaker (Larimichthys crocea) larvae and juveniles under five developmental feeding stages (DAH3, DAH7, DAH12, DAH19, DAH49) to clarify mechanisms of early nutritional adaptation. Digestive enzyme assays, transcriptome sequencing, and metabolomics were integrated to compare physiological changes across diets. Protease activity increased sharply from DAH7–19 with the introduction of rotifers, Artemia, and copepods, while amylase and lipase activities rose at DAH19–49, reflecting enhanced carbohydrate and lipid utilization during transition to formulated feeds. Transcriptomic analysis showed that differentially expressed genes were enriched in pathways involving protein digestion, lipid and energy metabolism, and cell cycle regulation. The metabolomic analysis further highlighted dynamic changes in amino acid, lipid, carbohydrate, and vitamin metabolism, consistent with transcriptomic findings. The integrated analysis suggests that the coordinated modulation of digestive enzyme activities, gene expression, and metabolite profiles enabled a smooth transition from yolk dependency to live prey feeding and a subsequent use of artificial diets. These findings provide new insights into the early nutritional development of NZ large yellow croaker and provide a scientific basis for the improvement of artificial aquaculture seed production.

Figures

Full text

1. Introduction

The large yellow croaker (Larimichthys crocea) is an economically important marine species in East Asian countries, particularly China. The human population in China were divided into three geographical groups based on morphology: the Naozhou stock (South China Sea), the Min-yuedeong stock (Eastern South China Sea and Taiwan Strait), and the Diqu stock (East China Sea) [1,2,3]. Although disparate gene pools among these stocks have been described, recent population-based genomic analyses suggest that they were formed through climate changes whereby stock boundaries shifted, and genetic divergence at different latitudes might have also contributed to stock differentiation [4,5,6]. Within this framework, the Naozhou stock (hereafter NZ stock) was shown to be a relatively independent population at both genetic and phenetic levels. Establishing a comprehensive genomic resource is required to advance breeding programs and facilitate conservation [5,7,8,9]. As the NZ stock has a high level of genetic diversity with distinctive phenotypes, it is considered a valuable candidate for selection [8,9]. However, the digestive physiological response to varying diets in this stock is not well understood, presenting a barrier to the development of scientifically formulated diets.

Larimichthys crocea undergoes an ontogenetic shift from endogenous to exogenous feeding at 4–6 DAH [10], during which the buccopharyngeal and intestinal tracts and hydrolytic enzymes rapidly develop [11,12,13,14]. Rotifers are introduced at 3–8 DAH [15], followed by Artemia, copepods, and microdiets, reflecting progressive gut maturation [16]. As larvae develop, gastric glands begin secreting acid pepsin and intestinal villi expand [13]. Trypsin and lipase activities rise at 10–20 DAH [17], and juveniles achieve full digestive capacity by ~50 DAH [12,13]. Nevertheless, this stage is sensitive due to immature digestion and fluctuating enzyme activity, contributing to high mortalities [7,18,19]. Balanced larval diets are therefore essential to improve survival and seed quality [16,17].

Live feeds remain standard but present challenges of cost, supply, and pathogen risks [16,17]. Plant protein sources, especially SBM and FSBM, provide sustainable FM replacements, though high SBM levels may induce enteritis and immune activation [20,21]. Fermentation reduces anti-nutrients and improves safety [22,23], while multi-omics studies show plant proteins reshape metabolic and immune pathways [24,25]. Functional additives such as phospholipids and oryzanol enhance survival, enzyme activity, and gut health [25,26,27]. Given phenotypic and genetic distinctions of the NZ stock [5,8,9], diet responses may be stock-specific, yet few studies target this stock [28], underscoring the need to assess its digestive physiology under alternative diets.

Despite advances in genomic characterization of L. crocea, the digestive physiology and early nutritional adaptation of the NZ stock remain poorly understood. This knowledge gap presents a critical limitation for stock-specific feed formulation, larval rearing, and selective breeding programs. Early developmental stages of L. crocea are particularly sensitive, with digestive enzyme activity, nutrient absorption, and metabolic regulatory mechanisms undergoing rapid ontogenetic transitions. These shifts influence larval survival, feed efficiency, and the success of artificial seed production. While previous studies have examined general digestive development or dietary responses in L. crocea, they have rarely accounted for stock identity, and no study has integrated multi-omics approaches to understand the NZ stock’s unique digestive physiology under changing diets. This study represents the first comprehensive multi-omics investigation combining digestive enzyme assays, transcriptomics, and metabolomics to characterize the developmental and diet-dependent digestive physiology of the NZ stock of large yellow croaker. Unlike earlier research that focused on single omics layers or non-stock-specific observations, our integrated approach enables a holistic understanding of how genes, enzymes, and metabolites jointly shape nutritional adaptation from endogenous yolk reliance to live prey feeding and finally to formulated diets. Such a multi-layered perspective is essential for addressing key bottlenecks in NZ stock aquaculture. Stock-specific physiological insights can guide the design of more appropriate formulated diets, improve larval survival and growth, and support selective breeding strategies aligned with the stock’s unique genomic and phenotypic characteristics. This study aims to investigate the stock-specific nutrient requirements needed to maintain intestinal mucosal health in the Naozhou (NZ) stock by integrating enzyme assays with transcriptomic, metabolomic, and differential gene expression analyses. The findings will support the development of sustainable aquaculture feeds that align with the NZ stock’s genomic and phenotypic traits [24,29,30] and guide the formulation of efficient, biosecure diets that promote healthy fish populations and environmentally responsible aquaculture.

2. Materials and Methods

2.1. Nutritional Enrichment and Artificial Spawning of the NZ Stock

One thousand one hundred and ninety five (1195) wild Naozhou stock broodstock individuals of the NZ stock of large yellow croaker were collected from the coastal waters of Naozhou Island, Zhanjiang, Guangdong Province. After acclimation, the fish were temporarily maintained in land-based concrete tanks (30 m^3^) with controlled rearing water conditions (23–25 °C, 28–30 salinity (ppt), and >6.0 mg·L^−1^ dissolved oxygen) with continuous aeration and daily water exchange at 50% of the tank volume. Shading was provided during acclimation, maintaining a photoperiod of 12 h light and 12 h dark. Only broodstock individuals in good health, with no external injuries and confirmed gonadal development status by biopsy, were selected for the nutritional conditioning phase. Seven hundred and ten (710) parent fish were able to spawn artificially after domestication (including 181 males and 529 females).

To promote gonadal development and enhance reproductive performance, the broodstock was subjected to a 4-week nutritional conditioning period prior to induced spawning. For this, the balanced diet was mainly composed of fresh ragworms (Nereis spp.; 60% of the total daily ration) and supplemented with frozen fish (20%) and formulated feed (20%). The fish were fed twice daily (at 08:30 and 17:30) at 3–4% of their body weight, and a vitamin premix (100 mg·kg^−1^ feed) and fish oil (2%) were regularly added to the diet.

Robust and healthy broodstock individuals with gonadal development at stage IV were selected for induced spawning, with a female-to-male ratio of 1:2. Hormonal induction was performed via intramuscular injection at the base of the pectoral fin using a combination of human chorionic gonadotropin (HCG, 800–1000 IU kg^−1^ body weight) and luteinizing hormone-releasing hormone analog (LHRHa, 20–30 µg.kg^−1^ body weight). Females received a single injection, while males were administered half the dosage. After injection, the fish were kept under flow-through and aerated conditions (as above but with a slightly higher temperature of 24–25 °C). Typically, 12–16 h post-injection, the females exhibited abdominal swelling and softness with clear signs of ovulation, while the males reached synchronous spermiation. At this stage, the fish were examined every 2 h to determine the optimal timing for gamete collection. Mature eggs and milt were obtained from females and males, respectively, by gentle abdominal stripping. Fertilization was conducted using the dry method, whereby the eggs and milt were mixed at a ratio of 1:100–150, followed by activation with a small volume of clean seawater. The fertilized eggs were immediately transferred to indoor concrete tanks for incubation. At each developmental stage, three pooled biological replicates (each derived from one rearing tank) were used for both digestive enzyme assays and omics analyses. The number of individuals per replicate varied with fish size smaller larvae required pooled samples of tens of individuals, while larger juveniles required fewer (5–10) fish per replicate.

2.2. Egg Hatching and Seedling Cultivation

Hatching buckets were maintained with continuous aeration and a flow-through (as above after injection) and the rearing water was filtered through an 80 μm mesh and exchanged daily at 50% of the tank volume to ensure stable water quality. Hatching was completed within 24 h after fertilization, and by 3 DAH, the larvae had fully absorbed their yolk sacs and initiated exogenous feeding. From 4–7 DAH, larvae were fed 2–3 individuals.mL^−1^ of small rotifers (Brachionus plicatilis), which were enriched for 12 h prior to feeding, four times per day. At 8 DAH, 1–2 individuals.mL^−1^ of Artemia nauplii were gradually introduced as the main diet, while the proportion of rotifers was reduced. By 12 DAH, Artemia completely replaced rotifers. Before the transition from the larval to juvenile stage (around 13 DAH), 2–3 individuals.mL^−1^ of copepod nauplii were gradually introduced. Initially, larvae were co-fed with copepods and Artemia to facilitate the diet transition, and by 19 DAH, copepods became the dominant live feed. From 20 DAH onward, the larvae were gradually weaned onto a formulated diet (Supplementary Materials) comprising a combination of live feed and a microdiet. The microdiet contained 45% crude protein and 10% crude lipids (particle size 150–250 μm), of which small amounts were mixed with copepods, progressively increasing until a complete copepod replacement was achieved (Supplementary Materials). Both larvae and juveniles were fed 3–4 times daily at approximately 5% of their body weight. Throughout the rearing period, water quality was maintained at 24–26 °C with the same salinity and dissolved oxygen content as above. Total ammonia nitrogen was kept below 0.2 mg·L^−1^ and nitrite below 0.02 mg·L^−1^. Feces and uneaten feed were removed daily, and probiotics were regularly added to the tanks to stabilize water quality. Growth performance (body length and weight) and survival rate were monitored periodically to evaluate feeding efficiency and success at different dietary stages.

2.3. Larval and Juvenile Sampling

Based on the developmental stages and feeding characteristics of the fish larvae and juveniles, five representative sampling points based on number of DAH were established: DAH3 (unfed, yolk sac not fully absorbed), DAH7 (primarily feeding on rotifers, DAH12 (feeding mainly on Artemia nauplii), DAH19 (feeding on copepods), and DAH49 (feeding on a formulated diet). Prior to sampling, fish were starved for 24 h to empty the digestive tract and minimize interference from feed residues. At each time point, individuals were anesthetized with MS-222 (100 mg·L^−1^) manufactured by Syndel (a subsidiary of Aquatic Life Sciences, Inc.) Ferndale, MI, United States and was sourced in China., immediately frozen in liquid nitrogen, and stored at –80 °C for subsequent digestive enzyme activity assays and omics analyses.

2.4. Digestive Enzyme Activity Assays

The digestive enzyme activities of both larvae and juveniles were measured using commercial assay kits (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). The activities of alkaline phosphatase (AKP), trypsin, pepsin, lipase (LPS), amylase (AMS), and leucine aminopeptidase (LAP) were determined.

For sample preparation, the frozen specimens stored in liquid nitrogen were slowly thawed on ice and accurately weighed. Each sample was homogenized for 3–5 min in ice-cold 0.9% NaCl solution at a ratio of 1:9 (tissue weight (g): buffer volume (mL)) using a glass homogenizer in an ice-water bath until fully disrupted. The homogenates were centrifuged at 2500 r/min for 10 min at 4 °C and the supernatant was collected for enzymatic assays.

The enzyme activities were measured following the manufacturer’s instructions: AKP using p-nitrophenyl phosphate (pNPP) as the substrate, with the absorbance measured at 405 nm; trypsin and pepsin using casein as the substrate, with the absorbance measured at 280 nm; LPS colorimetrically using p-nitrophenyl ester substrates, with the absorbance measured at 420 nm; AMS using soluble starch as the substrate and iodine-starch colorimetry at 660 nm; LAP using L-leucyl-p-nitroanilide as the substrate, with the absorbance measured at 405 nm.

All enzyme activities were expressed as specific activity (U/mg prot), normalized to protein content. Protein concentrations were determined using the Bradford method with Coomassie Brilliant Blue G-250, using bovine serum albumin (BSA) as the standard. Three replicates were sampled at each sampling point, and all assays were conducted at a constant 25 °C to ensure reproducibility and accuracy.

2.5. Transcriptome Sequencing and Analysis

2.5.1. RNA Extraction, Library Preparation, and High-Throughput Sequencing

A total of 6 biological replicates were used for the metabolomic analysis. Total RNA was extracted from the fish samples using the TRIzol reagent (Life Technologies, Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocol, and genomic DNA was extracted with DNase I (TRANS). RNA integrity and potential DNA contamination were first examined by agarose gel electrophoresis. RNA quality was further assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), its purity was determined by a NanoDrop spectrophotometer Guangzhou, China. (OD260/280 and OD260/230 ratios), and the concentration and yield were quantified with a Qubit 2.0 Fluorometer (Toshima City, Japan) [31].

RNA qualifying for analysis was diluted to 1.5 ng.μL^−1^ for cDNA library construction. Poly(A)^+^ mRNA was enriched using oligo(dT)-attached magnetic beads and then fragmented into shorter pieces through ultrasonic disruption. Using the fragmented mRNA as a template, first-strand cDNA was synthesized with reverse transcriptase, followed by second-strand cDNA synthesis with DNA polymerase I and dNTPs after RNase H treatment. The resulting double-stranded cDNA was then purified, subjected to end repair, A-tailing, and adaptor ligation. Target fragments of approximately 200 bp were screened using AMPure XP beads, Santa Clara, CA, USA and sourced in China, followed by PCR amplification and purification to obtain high-quality cDNA libraries.

Three replicate libraries were constructed for each of the five sampling time points. After quality control, all libraries were sequenced on the Illumina NovaSeq 6000 platform, San Diego, CA, USA to generate high-throughput transcriptome data (Reference Genome Information: Larimichthys crocea isolate HW-2024, whole-genome shotgun sequencing project. ACCESSION JBIBPI000000000, VERSION JBIBPI000000000.1 and BioProject: PRJNA1169539 with BioSample: SAMN44078122. Web link: https://www.ncbi.nlm.nih.gov/nuccore/JBIBPI000000000.1/, accessed on 23 December 2025). Library preparation and sequencing were conducted by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China).

2.5.2. Raw Data Quality Control and Filtering

To ensure reliability, image data generated by Illumina NovaSeq 6000 were first processed by the system’s built-in software for base calling and converted into FASTQ format raw sequencing reads. Prior to downstream analyses, strict quality control and filtering of the raw data by removing (1) reads containing adaptor sequences, (2) reads with more than 10% unknown bases (N), (3) reads in which more than 50% of the bases had a Phred quality score (Q) below 20, and (4) low-complexity reads composed of a single nucleotide (e.g., continuous poly-A). The resulting high-quality sequences were defined as clean reads and used as the basis for subsequent analyses. In addition, the quality of the clean reads was evaluated using FastQC software Version 0.12.0 and in-house scripts, assessing indicators such as Q20, Q30, and GC content to ensure the accuracy and reliability of the sequencing data.

2.5.3. Differentially Expressed Gene (DEG) Analysis

After transcriptome sequencing, reads were mapped to the reference genome, gene expression levels were quantified, and DEGs were identified using DESeq2 (version 1.30.0) [32]. Raw read count data were first normalized to eliminate any biases caused by sequencing depth and gene length. A negative binomial statistical model was applied to assess significant expression differences between groups. Correction for multiple comparisons was performed using the Benjamini–Hochberg method to obtain the false discovery rate (FDR), thereby reducing false positives. Genes with FDR < 0.05 and |log_2_(FC)| > 1 (i.e., fold change, FC > 2 and statistically significant after FDR correction) were defined as DEGs. Subsequent multidimensional functional analyses were then performed including volcano plots to visualize the magnitude and significance of gene expression changes, a Gene Ontology (GO) enrichment analysis to identify significantly enriched categories in biological processes (BP), molecular functions (MF), and cellular components (CC), and a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis to determine significantly enriched pathways related to metabolism, signal transduction, and physiological functions.

2.5.4. Quantitative Real-Time PCR (RT-qPCR)

To validate the reliability of the RNA sequencing (RNA-seq) data, ten DEGs were randomly selected for the RT-qPCR analysis, with β-actin acting as the internal reference gene. Gene-specific primers were designed based on the cDNA sequences using Primer Premier version 5.0 (Table 1) and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). RT-qPCR assays were performed on an ABI QuantStudio 3 Real-Time PCR system, Waltham, MA, USA and sourced in China. Each 10 μL of reaction mixture contained 5 μL of 2× PerfectStart Green qPCR SuperMix, 3.2 μL of nuclease-free water, 1 μL of cDNA template, and 0.4 μL each of both forward and reverse primers (10 µmol·L^−1^). Amplification was carried out through initial denaturation at 95 °C for 30 s, 40 cycles at 95 °C for 5 s and at 60 °C for 15 s, then at 95 °C for 15 s, 60 °C for 1 min, and 95 °C for 1 s. All reactions were performed in triplicate, and the Ct value was analyzed using the 2^−ΔΔCt^ method [33].

2.6. Metabolomic Sequencing and Analysis

2.6.1. Metabolite Extraction

Samples were thawed slowly at 4 °C, and an appropriate amount of tissue (0.3 g) was transferred into pre-cooled methanol/acetonitrile/water solution (2:2:1, v/v). After vortexing, samples were subjected to low-temperature ultrasonication for 30 min, followed by incubation at −20 °C for 10 min. The mixture was then centrifuged at 14,000× g for 20 min at 4 °C, the supernatant was collected and vacuum-dried, and these dried extracts were reconstituted in 100 μL of acetonitrile/water solution (1:1, v/v), vortexed, and centrifuged at 14,000× g for 15 min at 4 °C. The resulting supernatant was used for the liquid chromatography–mass spectrometry (LC-MS) analysis.

Chromatographic separation was performed using an Agilent 1290 Infinity LC UHPLC system equipped with a HILIC column. A 2 μL aliquot of each sample was injected with an autosampler, and separation was carried out under gradient elution with a flow rate of 0.5 mL.min^−1^ and a column temperature of 25 °C. The mobile phases consisted of (A) water +25 mM ammonium acetate +25 mM ammonium hydroxide, and (B) acetonitrile and the gradient program followed 95% B for minute 0–0.5, a linear decrease in B from 95% to 65% for minute 0.5–7, a linear decrease in B from 65% to 40% for minute 7–8, 40% B for minute 8–9, a linear increase in B from 40% to 95% for minute 9–9.1, and 95% B for minute 9.1–12. During this time, the samples were kept in the autosampler at 4 °C. To minimize signal drift, samples were randomly analyzed, with quality control (QC) samples inserted into the sequence to monitor instrument stability and data reliability.

2.6.2. Identification of Differential Metabolites

The method was slightly modified based on the method described by [34]. Raw LC-MS data were converted into MzML format using ProteoWizard software (v3.0.8789), and peak detection, retention time alignment, and peak area extraction were performed using XCMS software (Version 3.7.0) [35], with parameter settings as follows: for peak picking centWave m/z tolerance = 10 ppm, peakwidth = c(10, 60), prefilter = c(10, 100); for peak grouping bw = 5, mzwid = 0.025, minfrac = 0.5. The extracted data were checked for completeness, and total ion current (TIC) normalization was applied to ensure comparability across samples and metabolites. Metabolites with more than 50% missing values were excluded from further analysis and remaining missing values were imputed using the K-Nearest Neighbors (KNN) method [36,37]. Outliers were removed, and normalized peak areas were used for subsequent statistical analyses.

A data matrix was generated containing retention time (RT), mass-to-charge ratio (m/z), and peak intensity information. Metabolite annotation was achieved through matching them against public databases including MassBank, Metlin, and MoNA, in combination with an in-house secondary mass spectrometry database.

To evaluate biological reproducibility among samples, a correlation analysis was performed, followed by a principal component analysis (PCA) and orthogonal partial least squares-discriminant analysis (OPLS-DA) comparing experimental and QC samples [38]. A 200-time permutation test was used to assess potential model overfitting. Differential metabolites between groups were identified based on the intersection of three criteria: FC > 1, p < 0.05 (Student’s t-test), and variable importance in projection (VIP) > 1 derived from the OPLS-DA. Pathway annotation of differential metabolites was conducted using the KEGG database and identified metabolites were classified using the KEGG Compound and Human Metabolome Database (HMDB) for compound category statistics [39].

2.7. Data Analysis

The enzyme activity data were statistically analyzed using GraphPad Prism 10 (GraphPad Software, Boston, MA, USA, www.graphpad.com). Differences among groups were evaluated by one-way analysis of variance (ANOVA), followed by post hoc multiple comparisons when the assumption of homogeneity of variance was met. A significance level of p < 0.05 was used and all results are presented as mean ± standard deviation (SD).

3. Results

3.1. Diet-Based Effects on Digestive Enzyme Activities

Changes in the digestive enzyme activities of NZ large yellow croaker at different developmental stages are shown in Table 2. As indicated, AKP, LAP, AMS, pepsin, and trypsin in the DAH49 group (formulated diet stage) were significantly higher than those in the DAH3, DAH7, DAH12, and DAH19 groups (p < 0.05). During the live feed stages (DAH7–DAH19), trypsin activity was significantly higher than in the DAH3 group (unfed stage) (p < 0.05), suggesting that ingestion of exogenous feed markedly enhanced protein digestion capacity. In contrast, LPS activity did not differ significantly among developmental stages (p > 0.05).

3.2. Transcriptome Analysis

3.2.1. RNA-Seq Quality Assessment

Transcriptome sequencing was performed for all groups, with three biological replicates per group. A total of 676,738,462 bp of raw sequencing data were generated and after QC and removal of low-quality reads, 676,519,612 bp of clean data were retained (Table 3). A base quality analysis showed that Q20 values exceeded 99.36% and Q30 values exceeded 97.57% for all samples, indicating high sequencing accuracy. The GC content ranged from 47.63% to 48.64%, falling within a reasonable range. Overall, the data quality was sufficient for subsequent transcriptome assembly and differential expression analysis.

3.2.2. Identification of DEGs

Significant DEGs were detected in all pairwise comparisons across the different feeding stages (Figure 1) comprising 8071 between DAH3 and DAH7 (5399 upregulated and 2672 downregulated), 3426 between DAH7 and DAH12 (2079 upregulated and 1347 downregulated), 2773 between DAH12 and DAH19 (1138 upregulated and 1635 downregulated), and the most pronounced difference between DAH19 and DAH49 with 11,278 DEGs (2638 upregulated and 8640 downregulated).

Overall, throughout development and across varying diets, the number of significant DEGs decreased, followed by an increase, and then decreased again. During early developmental stages (DAH3–DAH12), the number of upregulated genes was consistently higher than downregulated genes, suggesting activation of genes related to digestive system development and adaptation to exogenous nutrition. In contrast, during later developmental stages (DAH12–DAH49), the number of downregulated genes was higher, particularly during the transition from DAH19 to DAH49, indicating that certain genes with high expression in early developmental stages were progressively suppressed.

3.2.3. Gene Ontology Analysis of Differentially Expressed Genes

To investigate the potential biological functions of DEGs and their roles in larval–juvenile development of the NZ stock, a GO enrichment analysis was performed (Figure 2). The results showed that DEGs from all five feeding stage comparisons were classified into the three major GO categories (BP, MF, and CC), encompassing a total of 48 subcategories.

For BP, DEGs were mainly enriched in fundamental biological activities such as cellular and metabolic processes, biological regulation, and responses to stimuli. For MF, enrichment was dominated by binding and catalytic activity, indicating that many DEGs are involved in substrate binding and catalytic reactions. For CC, DEGs were primarily associated with cellular anatomical structures and protein-containing complexes. Overall, MF DEGs contained the most enriched genes, suggesting that enzyme activity and related molecular binding functions play pivotal roles for the larval-juvenile stages in regulating digestive metabolism and development across different feeding stages.

3.2.4. Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis

To investigate the functional roles of DEGs during the larval–juvenile development stages of the NZ stock, DEGs were mapped to the KEGG database for an enrichment analysis (Figure 3).

When comparing DAH3 and DAH7 (Figure 3A), 366 enriched pathways were identified, of which 47 were significantly enriched. These were mainly associated with muscle cytoskeleton, motor proteins, glutamatergic synapse, and protein digestion and absorption, which are closely related to cell proliferation, tissue differentiation, and energy supply, providing the material and energy required for early larval development.

When comparing DAH7 and DAH12 (Figure 3B), 364 pathways were enriched, 37 of which were significantly enriched. These were primarily related to metabolic pathways, the cell cycle, and pancreatic secretion, and the associated genes were involved in the metabolism of steroids, carbohydrates, and lipids, and lipid digestion and bile secretion (digestive processes), as well as immune-related pathways including antigen processing and presentation. These changes ensured efficient nutrient use and supported the establishment of the immune defense system.

Comparing DAH12 and DAH19 (Figure 3C), 360 pathways were enriched, with 47 showing significant enrichment. These were mainly associated with muscle cytoskeleton, extracellular matrix–receptor interactions, and pyrimidine metabolism, suggesting that larvae-enhanced protein digestion and absorption support rapid increases in body length and weight, cell cycle regulation, and DNA replication ensure normal cell proliferation and division, and genome stability is maintained by genetic material repair. These processes are all critical for somatic growth and cellular development.

Comparing DAH19 and DAH49 (Figure 3D), 371 enriched pathways were identified, of which 99 were significantly enriched. These included metabolic pathways, cell adhesion molecules, muscle cytoskeletons, oxytocin signaling pathways, and adrenergic signaling in cardiomyocytes, driving the development of major organs (nervous and circulatory systems) and nutrient digestion and utilization.

Overall, the identified DEGs in both larvae and juveniles were mainly enriched in pathways related to digestion and absorption, energy metabolism, and tissue development. Early developmental stages were characterized by the establishment of nutrient metabolism and absorption, whereas later developmental stages also included cell cycle regulation, signal transduction, and the functional maturation of multiple organ systems. This indicates that both larvae and juveniles adapt to sudden environmental changes through coordinated regulation of diverse metabolic and signaling pathways throughout development.

3.2.5. Quantitative Reverse Transcription Polymerase Chain Reaction Validation Results

To validate the accuracy of the RNA-seq data, 10 genes were randomly selected from each developmental stage for the RT-qPCR analysis. The results showed that the expression trends of these genes were highly consistent with those obtained from RNA-seq (Figure 4), supporting the reliability and reproducibility of the transcriptome sequencing data used in this study.

3.3. Metabolomic Results

3.3.1. Sample Quality and Multivariate Statistical Analysis

A PCA was performed for all sample groups under both positive and negative ion (NEG) modes. The results showed that the QC samples clustered closely together (Figure 5), indicating data stability and reproducibility, validating the reliability of subsequent metabolomic analyses.

3.3.2. Orthogonal Partial Least Squares Discriminant Analysis

To investigate the metabolic profile differences among feeding groups, OPLS-DA models were constructed. Under both positive and negative ion modes, all groups exhibited clear distinction (Figure 6, Figure 7, Figure 8 and Figure 9A,B), indicating significant metabolic differences associated with larval-juvenile development and dietary changes.

To verify the robustness and reliability of the models, 200 permutation tests were conducted by randomly shuffling the class variable Y. The R^2^ and Q^2^ of the randomized models were markedly lower than those of the original models (Figure 6, Figure 7, Figure 8 and Figure 9C,D), and the Q^2^ regression line and y-axis intercept were less than zero. This demonstrates the reliability of the OPLS-DA models through good fitting performance and predictive ability, validating its use for the identification and functional interpretation of differential metabolites.

3.3.3. Screening and Identification of Differential Metabolites

In the positive ion mode (POS; Figure 10A), 360 differential metabolites were identified between DAH3 and DAH7 (218 upregulated and 142 downregulated), 342 between DAH7 and DAH12 (78 upregulated and 264 downregulated), and 288 between DAH12 and DAH19 (75 upregulated and 213 downregulated). The most pronounced differences were observed between DAH19 and DAH49 with 435 identified differential metabolites, of which 346 were upregulated and 89 were downregulated.

In the NEG mode (Figure 10B), 362 differential metabolites were identified between DAH3 and DAH7 (262 upregulated and 100 downregulated), 369 between DAH7 and DAH12 (83 upregulated and 286 downregulated), 308 between DAH12 and DAH19 (87 upregulated and 221 downregulated), and 399 between DAH19 and DAH49 (327 upregulated and 72 downregulated).

Overall, the number of differential metabolites, as well as those up- and downregulated, represented stage-specific changes. From DAH3 to DAH7, the number of upregulated metabolites was greater than downregulated ones, from DAH7 to DAH19, downregulated metabolites predominated, and from DAH19 to DAH49, the number of upregulated metabolites increased again. This pattern reflects the dynamic metabolic shifts in both the larvae and juveniles during different developmental stages and feeding conditions, transitioning from early active regulation to mid-stage adaptation to enhanced utilization in later stages.

3.3.4. Differential Metabolite and Pathway Analysis

To investigate the metabolic characteristics of larvae at different feeding stages, a KEGG pathway enrichment analysis was performed on the identified differential metabolites (Figure 11). These differential metabolites were mainly grouped into six major biological categories, including metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, and human diseases. These involved multiple key pathways, including the metabolism of amino acids, lipids, carbohydrates, nucleotides, cofactors, and vitamins, as well as membrane transport and the digestive system.

From DAH3 to DAH7, differential metabolites were significantly enriched in amino acid and lipid metabolism pathways. In particular, elevated levels were observed for taurochenodeoxycholate, taurolithocholate, isoleucine, taurine, N-acetyl-L-phenylalanine, oleic acid, 1-palmitoyl-sn-glycero-3-phosphocholine, and creatine, while lower levels were observed for palmitic acid and 1-stearoyl-2-hydroxy-sn-glycero-3-phosphocholine (Table 4). This suggests that larvae in the early developmental stages relied heavily on amino acid and lipid metabolism to support rapid development. From DAH7 to DAH12, differential metabolites were significantly enriched in nucleotide and lipid metabolism pathways. In particular, hypoxanthine levels increased, while choline cation, 1-palmitoyl-sn-glycero-3-phosphocholine, adenine, leucine, 1-hexadecyl-sn-glycero-3-phosphocholine, trimethylamine N-oxide, cadaverine, taurochenodeoxycholate, 8S-hydroxy-5Z,9E,11Z,14Z-eicosatetraenoic acid, taurocyamine, and pantothenate significantly decreased (Table 5), reflecting gradual improvement of nucleotide synthesis and digestive functions accompanied by decreased activity in certain lipid and amino acid metabolic pathways. From DAH12 to DAH19, differential metabolites were primarily associated with energy and vitamin metabolism where acetyl-DL-carnitine, vitamin C, and taurochenodeoxycholate levels increased, while leucine, nicotinamide, lithocholic acid, pantetheine, succinate, and glycocholic acid levels decreased (Table 6). This indicates that vitamin and energy metabolism pathways were activated in the larvae to aid in rapid growth, even though some amino and bile acid metabolic activities were reduced. From DAH19 to DAH49, differential metabolites were enriched in lipid and steroid metabolism pathways where nicotinamide, glycochenodeoxycholate, trimethylamine N-oxide, glycocholic acid, thiamine cation, 11-dehydrocorticosterone, taurolithocholate, taurocholate, isoleucine, parsilin, and taurine levels significantly increased, while deoxycholic acid levels decreased (Table 7). This suggests that after transitioning to a formulated diet, lipid metabolism and bile acid biosynthesis were enhanced by the larvae, thereby improving digestion and supporting organ system development.

Overall, the KEGG enrichment analysis revealed that as feeding regimes shifted during larval development, metabolic regulation transitioned from amino acid and energy metabolism in early stages toward lipid and vitamin metabolism, and multi-system coordinated development in later stages. This reflects the dynamic metabolic adaptations underlying rapid growth and nutritional adjustment.

At different larval-juvenile developmental stages, the differential metabolites exhibited distinct pathway enrichment patterns (Figure 12). From DAH3 to DAH7, 152 pathways were enriched (27 significantly), mainly involving aminoacyl-tRNA, amino acid, and plant secondary metabolite biosynthesis, carbon and 2-oxocarboxylic acid metabolism, and protein digestion and absorption. From DAH7 to DAH12, 137 pathways were enriched (28 significantly), primarily glycerophospholipid, 2-oxocarboxylic acid, and α-linolenic acid metabolism, and amino acid biosynthesis. From DAH12 to DAH19, 138 pathways were enriched (25 significantly), mainly 2-oxocarboxylic acid metabolism, amino acid and aminoacyl-biosynthesis, and protein digestion and absorption. From DAH19 to DAH49, 151 pathways were enriched (17 significantly), mainly related to aminoacyl-tRNA and plant secondary metabolite biosynthesis, protein digestion and absorption, carbon metabolism, and mineral absorption.

Overall, metabolic pathway enrichment exhibited stage-specific shifts with the early stages dominated by amino acid and carbon metabolism, the midway stages showing progressively enhanced lipid metabolism, and the later stages transitioning toward improved nutrient absorption and multi-system metabolism.

The joint analysis revealed significant alterations in metabolic pathways during the early developmental stages of the Dazhu strain of yellow croaker, with particularly pronounced changes in the “metabolic pathway” category. As shown in Figure 13, nearly 870 genes were enriched in this pathway, significantly higher than the enrichment level of metabolites, suggesting the potential existence of a complex metabolic regulatory network. In contrast, pathways such as “neurological ligand-receptor interactions,” “cAMP signaling,” and “PI3K-Akt signaling” exhibited moderate enrichment of metabolites and genes but demonstrated relatively limited regulatory effects. Furthermore, pathways like “Dopaminergic Synapse” and “Cholinergic Synapse” showed lower involvement. These findings suggest that “metabolic pathways” play a key regulatory role in transcriptomics-metabolomics association analysis and represent the most prominent pathways for metabolic changes during early developmental stages.

The analysis of the results further indicate (Figure 14): Genes farther from the origin (such as ABCA12, MEP1B.3, Msln, SLC38A3.1, sec61al1) and metabolites (e.g., 4,4′-Sulfonylbis(2-methylphenol), Fexinidazole, Xanthine) indicate stronger associative effects on another omics dataset (metabolome or transcriptome), potentially serving as key factors in cross-omics regulation. The dispersed distribution of genes and metabolites suggests significant variations in contribution across different variables within the joint analysis. For instance, certain genes exhibit high loadings on both principal components, while some metabolites show lower loading values, potentially reflecting their lesser contribution to the analysis or interference from other factors.

4. Discussion

4.1. Diet-Based Effects on Digestive Enzyme Activities

In this study, the digestive enzymatic activity, transcriptomic expression, and metabolomics of larvae and juveniles of the NZ stock of L. crocea were closely associated with diet. Protease activity was significantly enhanced in the early-stage post-exogenous feeds (DAH7–19) for live high protein prey digestion, particularly during the copepod diet and formulation stage, indicating better carbohydrate and lipid utilization. These enzyme-driven changes align with the transcriptomic enrichment of protein digestion, lipid metabolism, and energy pathways, as well as changes in amino acid, lipid, and carbohydrate content. The integration of the enzymatic, transcriptomic, and metabolomic data highlights the adaptive flexibility of L. crocea regarding diet shifts during early developmental stages.

Previous studies have suggested that ingesting live prey partly influences the digestive enzyme activity of fish larvae [40,41]. Accumulating evidence indicates that exogenous contributions are minimal to the digestive enzyme system of larvae and juveniles, which is predominantly regulated by endogenous genetic mechanisms and the progression of digestive organ development [42,43]. For the red drum (Sciaenops ocellatus), the presence of digestive enzymes precedes the onset of initial feeding behavior, indicating that gene expression regulation is the core mechanism governing the early physiological preparation for digestion [44]. Enzyme secretion depends on developmental stage, nutritional requirements, and life mode [45,46], with key changes occurring at diet (yolk to feed) and developmental (larvae to juvenile) shifts [47,48]. Together, these studies indicate that enzyme activity is influenced by both genetic and environmental effects. At DAH3-7, AKP, LAP, LPS, AMS, trypsin, and pepsin levels remained low, supporting yolk nutrient contribution, as seen in long-spine porcupinefish (Diodon holocanthus) where these maternal reserves were used [49]. These basal levels primarily mobilize the yolk but provide pre-adaptation for subsequent feeding stages. During DAH7-19, the enzyme levels increased, indicating the transition into exogenous feeding. Here, trypsin and pepsin increased with protein ingestion, LAP improved the absorption of amino acids, and AMS and LPS were adjusted for the digestion of carbohydrates and lipids. In other marine species, increased trypsin appears to be an adaptation for incomplete gastric function [50,51,52]. Furthermore, essential amino acids have been shown to support both protein synthesis and energy production [53,54]. During DAH19-49, the juveniles fully established adult-like digestion. Here, protein-digesting enzymes, trypsin and pepsin, reached a peak, high AKP and LAP activity was observed, and AMS and LPS aligned with diet components; these reflect intestine maturation and the efficient application of formulated diets [55,56,57,58,59].

For the olive flounder (Paralichthys olivaceus), low fishmeal substitution with insect meal (black soldier fly) and/or other protein alternatives (FM35–45) did not result in disease after feeding [60] despite potential changes in gut microbes. Growth and health measurements of olive flounder in experimental groups were the same as those from the control group consuming 70% fishmeal, indicating a gut response to dietary reformulation [60]. These results, in contrast to L. crocea where dietary changes resulted in rapid enzymatic and metabolic shifts, especially in larval stages, highlight increased developmental sensitivity. For juvenile southern catfish (Silurus meridionalis), formulated diets resulted in higher peak metabolic rates, faster gastric drainage, and a greater specific growth rate compared to those fed natural prey (loach meat), and the discrepancies were linked to macronutrient equilibrium and microbial community diversity [61]. Likewise, L. crocea demonstrated developmental stage-dependent digestion optimized for natural prey, which was slowly modified under formulated diets. In contrast to S. meridionalis, L. crocea showed progressive dietary transitions, indicative of species-specific digestive plasticity. Larvae of European seabass (Dicentrarchus labrax) exhibited adaptive feeding by transitioning from live food to microparticulate diets in response to increased availability of microparticulate feed in their environment [62]. As a result, decreased growth rates initially led to enzyme (synthetic) upregulation, but this response eventually enabled adaptation. Similarly, an induction with enzymes facilitated the transition to formulated diets in L. crocea, although the mortality risk was high during this time. Enzyme activity exhibited these developmental stage-specific patterns including early endogenous control, mid-stage diet stimulation, and late convergence to digestion as an adult. The increased enzyme activities facilitated rapid growth rates and provided direction for the formulation of feed and larval survival.

4.2. Transcriptome Enrichment Pathway Analyses

In L. crocea, the metabolic pathways, vital for energy production and biosynthesis, were enriched throughout development. In the early developmental stages, larvae relied on yolk proteins/lipids, which comprised low enzyme activity [63,64,65]. At the onset of feeding bile secretion, cAMP signaling, DNA replication, and cell cycle pathways were enriched, highlighting proliferation through increased protein intake. During DAH3–19, live feed (rotifers, Artemia, and copepods) supplied high-quality proteins to the diet [66,67], resulting in an increase in AMS activity, increasing carbohydrate use. During this stage, enriched pathways included steroid biosynthesis, fat digestion, antigen presentation, and cell cycle regulation, representing immune and endocrine maturation. By DAH49, the larvae had matured digestive systems; AMS enabled carbohydrate use, thereby sparing proteins, lipids provided polyunsaturated fatty acids (PUFAs) for neural and immune functions [68], and enriched insulin secretion and olfactory transduction reflected the energy balance and adaptive feeding behavior [69]. Thus, the pathways shifted from yolk mobilization (early stage) to exogenous protein/lipid use with immune development (mid stage) to efficient carbohydrate utilization and energy regulation (late stage). These developmental stages matched enzyme activity changes with basal activity in early stages (DAH3–7), protease-dominated activity in the mid stage (DAH7–19), and carbohydrate/lipid digestion in the late stage (DAH19–49). Together, transcriptome analyses highlighted developmental stage-dependent nutritional adaptation. A cross-species comparison of fishes indicated that carnivores have increased growth rates, gastrointestinal evacuation, and postprandial metabolic responses compared to omnivores and herbivores [70]. Omnivores, by contrast, have a greater diversity of intestinal microbes. For L. crocea, the increased early upregulation of proteases is consistent with that in other carnivorous fish species, while subsequent increases in AMS and LPS indicated omnivorous digestion adaptation to ensure easy consumption of formulated meals.

4.3. Nutritional Adaptation from a Metabolomic Perspective

This study demonstrated that the metabolic profiles of larvae and juveniles of the NZ stock of large yellow croaker exhibited clear developmental stage-specific characteristics across feeding regimes, suggesting that metabolic regulation is driven by digestive system development. Generally, differential metabolites were mainly enriched in pathways related to amino acid, lipid, carbohydrate, nucleotide, and cofactor/vitamin metabolism, consistent with the high energy supply and biosynthetic demands of young Atlantic salmon (Salmo salar) and grouper (Epinephelus spp.) [71,72,73,74,75].

During the early developmental stage, larvae were dependent on yolk nutrients and had limited digestive function. The significantly enriched pathways included DNA replication, the cell cycle, and bile secretion, reflecting the mobilization of energy and substrates supporting rapid cell division and preparation for exogenous feeding. Active amino acid metabolism underscores the central role of proteins in cell proliferation and tissue formation and suggests that essential amino acids may function as energy substrates during this period. This finding is consistent with other species, such as Atlantic cod (Gadus morhua) and European seabass (Dicentrarchus labrax), where essential amino acids served as key energy sources during early development [76,77,78,79,80,81]. Combined with low digestive enzyme activity, it can be inferred that metabolism, at this stage, was primarily genetically regulated, ensuring mobilization and utilization of endogenous reserves.

At the onset of live feeding (DAH3–19), metabolites were significantly enriched in pathways related to fat digestion and absorption, steroid biosynthesis, and antigen processing and presentation, indicating that the metabolic profiles were shaped by exogenous protein and lipid intake. Increases in trypsin and pepsin were observed, consistent with the enrichment of protein digestion and absorption pathways. The activation of lipid-related pathways highlights the importance of unsaturated fatty acids from live feeds for membrane synthesis, hormone production, and immune development [82,83,84]. These results are consistent with those from other populations of large yellow croaker [18] where amino acid and lipid metabolism influenced rapid growth and immune system maturation during live feeding stages.

By the formulated diet feeding stage (DAH49), metabolomic profiles were significantly enriched in pathways including carbon metabolism, insulin secretion, and energy regulation, reflecting digestive system adaptation to the nutritional composition of formulated diets. Increased AMS and LPS activities indicated enhanced carbohydrate and lipid digestive capacity, with carbohydrates becoming increasingly important energy sources, thereby contributing to the protein-sparing effect. Similar findings have been reported in rainbow trout (Oncorhynchus mykiss), where enhanced carbohydrate metabolism during formulated diet feeding reduced protein use for energy, thus improving growth efficiency [85,86,87]. Furthermore, olfactory transduction pathway enrichment suggests that larvae may rely on sensory mechanisms to adjust feeding behavior, improving adaptation to formulated diets.

This study revealed developmental stage-dependent nutritional shifts in larvae and juveniles of the NZ stock of large yellow croaker. Early larval stages were dependent on endogenous proteins and lipids for energy. During the mid-developmental stages, they mainly utilized exogenous protein and phosphorus balanced by the inclusion of some lipids, with the gradual development of immunity. During the later developmental stages, they showed improved carbohydrate digestion and energy dispensation, adapting well to the formulated diet. These shifts also revealed coordinated changes between digestive enzyme activity and metabolic pathways, indicating that metabolomics is a valuable tool for evaluating nutritional needs of fish and for future diet adaptation. Importantly, these findings suggest that gene expression changes concurrently observed with dietary changes may directly alter nutrient utilization and energy supply. Furthermore, when diets are changed, the body responds within minutes at the mRNA level and by producing gene products like digestive enzymes. During the transition from endogenous to exogenous feeding, transcriptomics showed a cross-stage expression of DEGs with enriched protein and fat digestion. This approach also identified multiple DEGs over the different feeding periods and the transition from endogenous to exogenous feeding involved DEGs enriched in protein and fat digestion, the cell cycle, and energy metabolism pathways. Metabolomics results were consistent with these findings, with stage-specific changes in amino acid, lipid, and derivative levels observed, with overlapping KEGG pathways. This agreement between gene expression and metabolite levels demonstrated that dietary changes activated digestion- and metabolism-related genes and directly altered nutrient utilization and energy supply. For example, from DAH7 to DAH19, protein and lipid digestion genes were upregulated with shifts in amino acid and bile acid metabolites, and from DAH19 to DAH49, active carbon metabolism and energy regulation were indicated by transcriptomics and increased carbohydrate utilization, and energy-related metabolites were observed through metabolomics.

Collectively, these comparisons underscore the similarities and differences in digestive physiology between fish species. Like olive flounder, L. crocea exhibits enzymatic adaptation, but unlike yellow catfish (Pelteobagrus fulvidraco) or flounder, the larval stages are exceptionally sensitive to food adaptation. A study on seabass larvae indicated that a functional enzymatic system does not always translate to growth success unless the artificial diet is nutritional and palatable. Furthermore, genera are studied at different scales implying that digestive enzyme maturation, the degree of host microbiota modulation, and specific dynamic actions are determining factors for aquaculture performance. For L. crocea, adjusting the feed to suit the developmental digestive stage may decrease mortality rates, improve seed quality, and enhance aquaculture sustainability.

5. Conclusions

This study systematically investigated the digestive physiology and nutritional adaptation of larvae and juveniles of the NZ stock of large yellow croaker across five representative feeding stages. By integrating digestive enzyme assays, transcriptomic profiling, and metabolomic analyses, we revealed developmental stage-dependent changes in digestive capacity and metabolic regulation. Early larvae primarily relied on yolk reserves, supported by low but detectable levels of digestive enzymes and enrichment of pathways related to cell proliferation and yolk mobilization. At the onset of exogenous feeding, protease activities increased, and protein and lipid metabolism pathways were significantly enriched, highlighting the critical role of dietary proteins and lipids for rapid growth and immune development. In later stages, AMS and LPS activities were enhanced, and carbohydrate and energy regulation pathways were enriched, reflecting adaptation to formulated diets and the establishment of a mature digestive system.

The consistency between the transcriptomic and metabolomic results demonstrated that dietary transitions induced transcriptional regulation and metabolite-level changes, providing a comprehensive perspective of the molecular and physiological mechanisms of nutritional adaptation. These findings deepen our understanding of digestive and metabolic strategies of L. crocea during early development and offer practical guidance for hatchery practices. This study has provided a scientific basis for developmental stage-specific diet design, feeding regime optimization, methods for larval survival improvement, and efficient transition from live prey to formulated feeds. Collectively, this work contributes to advancing seedling production efficiency and supports the sustainable development of large yellow croaker aquaculture farms.