Before the spawning season, wild female S. o’connori individuals were captured in the YLTR in Tibet, and cultivated S. o’connori individuals were collected at the Institute of Fisheries Science, Tibet Academy of Agricultural and Animal Husbandry Sciences (Tibet, China). Cultivated S. o’connori individuals were fed with a formula feed containing 35% crude protein and 8% lipids twice a day at 3% body weight each time for more than 1 year. A total of four wild-type and four cultured female S. o’connori individuals were randomly selected for the experiment. All experiments with fish in this study were conducted in accordance with the guidelines on the care and use of animals for scientific purposes, established by the College of Animal Science and Technology of Henan Agricultural University, Zhengzhou, China (approval number: No. 20190521). After euthanizing, the liver and oocyte were collected. Part of the tissues was rapidly snap-frozen in liquid nitrogen and stored at −80°C, while the other part was fixed in Bouin’s fluid overnight at 4°C.

The fixed samples were dehydrated in an alcohol gradient and xylene and then embedded in paraffin prior to being cut into 5-μm sections for HE staining. HE staining was performed as previously described (Feldman and Wolfe, 2014). Three discontinuous sections for each sample were selected and observed with a light microscope at 400X and 20X (Olympus, Japan).

For each tissue, approximately 50 mg of powdered sample was vortexed in 400 μl of cold mix (acetonitrile:methanol:water, 4:4:2) and then incubated for 60 min at −20°C. After centrifugation at 14,000 g for 20 min at 4°C, the supernatant was collected and dried in a vacuum desiccator. To end derivatization, the sample was vortexed in 100 μl of an acetonitrile:water mixture (1:1) for 1 min. After centrifugation at 14,000 g for 15 min at 4°C, 2 μl of supernatant was immediately analyzed using liquid chromatography electrospray ionization tandem mass spectrometry (LC-ESI-MS/MS).

Chromatographic separation was implemented on an ultra-high performance Liquid Chromatography (Agilent 1290 Infinity, Waldbronn, Germany). The column oven was maintained at 25°C. The flow rate was 0.3 ml/min. The injection volume for each sample was 2 μl. The mobile phase consisted of solvent A (water + 25 mM aqueous ammonia + 25 mM ammonia) and solvent B (acetonitrile). Gradient elution conditions were set as follows: 0–1 min, 95% phase B; 1–14 min, 95–65% B; 14–16 min, 65–40% B; 16–18 min, 40% B; 18–18.1 min, 40–95% B; 18.1–23 min, 95% B. The samples were continuously analyzed in a 4°C autosampler. The quality control samples were inserted into the sample queue to evaluate system stability throughout the experiment.

The metabolites were then analyzed using a high-resolution tandem mass spectrometry TripleTOF 6600 (AB SCIEX, United States) with an electronic spray ionization (ESI) source, operated in both positive and negative ion modes. The ESI source operation parameters were as follows: source temperature, 600°C; IonSpray Voltage Floating (ISVF), ±5500 V; ion source gas I (GSI), gas II (GSII), and curtain gas (CUR), 60, 60, and 30 psi, respectively. The time-of-flight (TOF) mass range was from 60 to 1000 Da, and the scan time was 0.2 s/spectra. The product ion scan range was from 25 to 1000 Da, and the product ion scan time was 0.06 s/spectra. The operation parameters of the high-sensitivity mode were as follows: declustering potential, ±60 V; collision energy, 35 ± 15 eV; excluding isotopes, within 4 Da; and candidate ions to monitor per cycle, 6.

The raw data were converted into mzXML format by ProteoWizard. Peak identification, filtration, and alignment were performed using XCMS program. The ion peaks with >50% missing values. The identification of each metabolite was performed according to the chemical formula generated from the accurate mass, with an accuracy of <25 ppm and MS/MS spectrogram matching. Metabolite identification was performed using a self-built database (Personalbio Technology Co., Ltd., Shanghai, China).

Principal component analysis (PCA), partial least square discriminant analysis (PLS-DA), and orthogonal partial least squares discriminant analysis (OPLS-DA) were performed after pattern recognition and pareto-scaling using SIMCA-P 14.1 (Umetrics, Umea, Sweden). The variable influence on projection (VIP) values were obtained from the OPLS-DA. Differentially accumulated metabolites (DAMs) were defined according to the following criteria: a VIP value > 1 and a q-value < 0.05. Clustering analyses of DAMs were applied based on Euclidean distance. Metabolic pathway analyses were performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.kegg.jp/). KEGG pathways with FDR < 0.05 were considered to be significantly enriched.

Total RNA was extracted using TRIzol reagent (Takara, Dalian, China) following the manufacturer’s instructions and treated with DNase I. RNA quantity and purity were assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States). Sequencing libraries were generated using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, United States) following the manufacturer’s recommendations. The library was purified (AMPure XP system) and quantified using the Agilent high sensitivity DNA assay on a Bioanalyzer 2100 system (Agilent Technologies, CA, United States). The libraries were then sequenced on a HiSeq platform (Illumina, San Diego, CA, United States) by Shanghai Personal Biotechnology Cp. Ltd., and 150 bp paired reads were generated. The raw sequence files were deposited to the National Center for Biotechnology Information Sequence Read Archive with accession number SRP300881.

Raw reads were filtered to remove adaptor-containing reads, low quality reads (mean Q-value <20), and reads containing more than 5% ambiguous nucleotides using SOAPnuke (V2.0.7) (Chen et al., 2018). Clean reads were obtained and aligned to the Schizothorax o’connori genome (Xiao et al., 2020) using Bowtie2 v2.2.5 (Langmead and Salzberg, 2012). The gene expression level was estimated based on the unique mapped reads by the reads per kilobase million mapped reads (RPKM) method using RSEM v1.2.8 (Li and Dewey, 2011). Hierarchical clustering analyses were performed based on the Pearson correlation and average linkage method. Differentially expressed genes were identified by DESeq2 (v1.16) with a fold change (FC) > 2 and a false discovery rate (FDR) < 0.05 (Love et al., 2014). GO enrichment analysis was performed using Blast2GO software (Conesa et al., 2005). Pathway analysis of DEGs was performed using the KEGG PATHWAY database (http://www.kegg.jp) and Blast2GO (Conesa et al., 2005). GO and KEGG terms with FDR < 0.05 were considered to be significantly enriched. The Pathview package was used to visualize the results of association analysis of differential metabolites and related transcripts.

qRT-PCR analysis was performed by using the ChamQ SYBR qPCR Master Mix (Vazyme Biotec, Nanjing, China) on a Light Cycler 480 system (Roche, Basel, Switzerland). A total of 10 DEGs related to fatty acid biosynthesis, biosynthesis of unsaturated fatty acids, and starch and sucrose metabolism, which are three significantly enriched pathways, are selected. Specific primers are listed in Supplementary Table S5. The expression levels of genes were normalized to the expression of the ACTB1 and calculated by the 2^−ΔΔCt method. All samples were carried out in three technical replicates on the same plate. Data were expressed as mean ± standard deviation. Pearson correlations were computed to evaluate the correlation between the RNA-seq and qRT-PCR data.

Compared with wild S. o’connori, captive-reared S. o’connori has a yellow body color and nonuniformly sized whitish ova (Figure 1). Histological examination further showed that compared with wild S. o’connori, there were many atresia oocytes in the ovaries of captive-reared S. o’connori. The chorion and follicular cell layer also showed irregular shrinkage and folding (Figures 2A–D). The liver cells were evenly distributed in wild S. o’connori, while the nucleus was absent or pushed to the edge of the cell (Figures 2E,F, red arrows), and the cytoplasm was liquefied in the liver of captive-reared S. o’connori (Figures 2E,F, black arrows).

Metabolites were further analyzed using a liquid chromatography electrospray ionization tandem mass spectrometry (LC-ESI-MS/MS) system. The retention time and peak intensities of the quality control (QC) samples were reproducible, indicating that LC-ESI-MS/MS was reliable (Supplementary Figure S1). After raw data were preprocessed, 5,129 metabolite peaks (2,671 peaks in ESI positive mode and 2,458 peaks in ESI negative mode) were obtained for further analysis. The results of PCA, PLS-DA, and OPLSDA showed that the wild and captive-reared groups were clearly different in both the liver and the oocytes, reflecting significant metabolic differences among the different groups (Table 1, Supplementary Figures S2, S3). All samples fell within the 95% confidence interval.

After comparison with the self-built database, a total of 233 metabolites were annotated. There were 12 and 25 differentially accumulated metabolites (DAMs) between oocyte from domestic S. o’connori (OD) and oocyte from wild S. o’connori (OW) in the ESI+ and ESI− modes, respectively (Supplementary Table S1). However, there were 25 and 45 DAMs between liver from domestic S. o’connori (LD) and liver from domestic S. o’connori (LW) in the ESI+ and ESI− modes, respectively (Supplementary Table S1). Of the 37 DAMs identified between OD and OW, 18 (48.65%) and 19 (51.35%) metabolites were downregulated and upregulated, respectively. These DAMs included many metabolites related to ovarian development, such as pregnenolone, prostaglandin F3a (PGF3a), arachidonic acid (ARA), adenosine 5′- diphosphate (ADP), and adenosine monophosphate (AMP). Of the 70 DAMs in LD compared to LW, 45 (64.29%) and 25 (35.71%) metabolites were downregulated and upregulated, respectively (Figures 3A,B). A hierarchical clustering analysis was further performed to assess the DAMs accumulation patterns (Figures 3C–F).

According to the previous KEGG analysis, only two pathways were significantly enriched between the OD and OW group. These included arginine and proline metabolism and linoleic acid metabolism. There were five significantly enriched pathways between LD and LW, including the biosynthesis of unsaturated fatty acids, pyrimidine metabolism, ABC transporters, starch and sucrose metabolism, and fatty acid biosynthesis (Supplementary Table S2).

To provide an overview of the system-wide changes in gene expression during artificial breeding, transcriptomic profiles were investigated using RNA-Seq. After data filtering, a total of 658.68 M clean reads were acquired, ranging from 39.4 to 45.1 M clean reads per sample. Among these, more than 88.56 and 72.16% were mapped to the S. o’connori genome and the exonic regions for each sample, respectively. However, the average percentage of uniquely mapped reads of exonic regions was only 21.10%, indicating large number of duplicated genes (Table 2). Saturation curve analysis showed that the results of gene quantification were reliable (Figure 4A). Correlation analysis and cluster analysis were performed to gain insight into the relationships among different samples. Almost all the biological duplicates were also clustered together, confirming the relevance of the samples (Figures 4B,C).

A total of 2,841 DEGs (1,285 upregulated and 1,556 downregulated) were found in LD compared with LW, and 159 DEGs (80 upregulated and 79 downregulated) were found in OD compared with OW (Figures 5A–C, Supplementary Table S3). All DEGs between the wild and captive-reared S. o’connori were classified using hierarchical clustering (Figures 5D,E).

To further understand the functions of these DEGs, GO and KEGG analyses were performed. Between LD and LW, a total of 449 GO terms and 102 pathways were significantly enriched, respectively (Supplementary Table S4). The top 10 significantly enriched GO terms and pathways are shown in Figure 6, including carbon fixation pathways in prokaryotes, PPAR signaling pathway neuroactive ligand-receptor interaction, propanoate metabolism, and carbon metabolism. However, only one GO term and one pathway (cell growth and death) were significantly enriched between OD and OW.

Based on the DEG pathway enrichment results, qRT-PCR of 10 DEGs (ACACA, ELOV6, FADS2, FAS, PGM1, SCD5, LYAG, G6PC3, HXK2, and ACBG1) related to fatty acid biosynthesis, biosynthesis of unsaturated fatty acids, and starch and sucrose metabolism was performed to validate the differential gene expression obtained by RNA-seq. The results showed a strong correlation between the RNA-seq and qRT-PCR data (R² = 0.90), indicating the reliability of the RNA-seq data (Figure 7).

Comparing the metabolomic and transcriptomic data showed that there were three common significant enrichment pathways between LD and LW, including biosynthesis of unsaturated fatty acids, starch and sucrose metabolism, and fatty acid biosynthesis. The DAMs and DEGs related to the fatty acid synthesis pathway are shown in Figure 8A. However, no common pathway was detected between OD and OW. Notably, two DEGs (CDC25B and FZR1) in the progesterone-mediated oocyte maturation pathway were enriched, although the entire pathway was not significantly enriched. We further analyzed the metabolite-gene pairs using pathway-based integration method. Detailed information about the metabolite-gene pairs is shown in Supplementary Table S6. Top 15 metabolite-gene pairs between LD and LW are shown in Figure 8B. Of these, two metabolite-gene pairs including pregnenolone-STS and pregnenolone-SULT2B were negatively correlated (Figure 8B). These two metabolite-gene pairs are related to steroid hormone biosynthesis pathway (Figure 8C). In addition, the levels of pregnenolone and LDLR were all decreased in captive-reared S. o’connori.