Chemicals and reagents
Water, methanol, acetonitrile, and formic acid were purchased from CNW Technologies GmbH (Düsseldorf, Germany). L-2-chlorophenylalanine was purchased from Shanghai Hengchuang Biotechnology Co., Ltd. (Shanghai, China). DNeasy PowerSoil Kit (Cat. No. 12888) and QIAamp 96 PowerFecal QIAcube HT Kit (Cat. No. 51531) were purchased from QIAGEN (Germany). Qubit dsDNA Assay Kit (Cat. No. Q32854) was purchased from Life Technologies (America). Tks Gflex DNA Polymerase (Cat. No. R060B) was purchased from Takara (Japan).
Fish culture conditions, treatments and sample collection
Two-month-old fish were raised in a greenhouse tank recirculating aquaculture system (RAS) and underwent sixth-generation selection in indoor concrete tanks at Hongze Fishseeds Biotech. Inc., Jiangsu, China. RAS can ensure the excellence of the water quality due to the physical filtration, oxygenation and bio-filtration (Liu et al., 2021). The system was a specially designed dual-drain RAS equipped with an online infrared dome camera to monitor the growth, maturation and natural spawning behavior of shad broodstocks. Other two-month-old fish were also raised in a greenhouse pond aquaculture system at Hongze Fishseeds Biotech. Inc., Jiangsu, China. Temperature and dissolved oxygen were measured by a YSI 550A (Y.S.I. Environmental Inc.) every morning in both locations. Water quality was monitored by testing ammonia and nitrite levels every four days in both locations. Both recirculating aquaculture system and the pond were maintained under a natural photoperiod in the greenhouse, with the sunlight intensity attenuated by a roof cover with 40% light penetration. The daily average light intensity ranged from 800 to 1600 lux. Fish were fed commercial pellets three times a day for 20 min each time during the entire culture period.
Sixty fish (two months postfertilization) were randomly selected and divided into three tanks or three ponds (20 fish per tank or pond). After thirteen month cultured, the ten-individual fish from each tank and pond (n = 3) sampled to measure the body weight. Then these fish were washed with purified water and euthanized. The intestinal parts were dissected and gathered for metabolites and RNA. Four biological replicates per group were used for LC-MS analysis and three biological replicates per group were used for RNA-Seq analysis. The sediment at the bottom of the tank aquaculture and pond aquaculture was also collected. Four tubes (>2mg/tube) of sediment were gathered for each group to ensure a sufficient number of DNA.
Sample preparation and analysis
Metabolite extraction
The process of metabolite extraction referenced Lu et al., (2016). The ACQUITY UHPLC system (Waters Corporation, Milford, USA) coupled to an AB SCIEX Triple TOF 5600 System (AB SCIEX, Framingham, MA) was used to reveal the metabolic profiles with ESI. The modes of positive and negative ion executed by using the C18 column. Water/formic acid and acetonitrile/formic acid were used to elute, and the separation parameter followed the instructions.
RNA extraction, library construction and sequencing
RNA was extracted from 3 mg of intestinal tissue using a mirVana miRNA Isolation Kit (Ambion). The RNA-Seq process was performed as described in Zhang et al. (2021). RNA integrity was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Samples with an RNA integrity number (RIN) ≥ 7 were used for subsequent analysis. Libraries were constructed using a TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. Then, these libraries were sequenced on an Illumina sequencing platform (Illumina HiSeq X Ten), and 125 bp/150 bp paired-end reads were generated.
DNA extraction and library construction
DNA was extracted from gut samples using TruSeq Nano DNA LT Sample Prepararion Kit (Illumina). The V4 region of bacterial 16S rRNA was amplified using DNA samples. Agencourt AMPure XP (BECKMAN COULTER, USA) was used for purifying PCR products. Sequencing libraries were constructed using KAPA Library Quantification Kits (Illumina, USA). The Illumina Miseq platform was used for paired-end sequencing (Kapa Biosystems, USA).
Data preprocessing and statistical analysis
LC-MS metabolomics data analysis
LC-MS metabolomic analysis was done by OE Biotech Co. Ltd. (Shanghai, China). Progenesis QI software was used for acquired LC-MS raw data analysis (Waters Corporation, Milford, USA). Progenesis QI Data Processing Software was used to analyze the metabolites (Waters Corporation, Milford, USA). The public databases were referenced to explore, including http://www.hmdb.ca/; http://www.lipidmaps.org/ and in-house-developed databases. The changes of metabolite were analyzed by Principal component analysis (PCA) and orthogonal partial least-squares-discriminant analysis (OPLS-DA). The Hotelling’s T2 region shows the 95% confidence interval in the model. Variable importance in the projection (VIP) shows the contribution rate of each variable in the OPLS-DA model. The differential metabolites were obtained based on the combination of a statistically significant threshold of variable influence on projection (VIP) values selected from the OPLS-DA model and p values from a two-tailed Student’s t-test on the normalized peak areas. The differential metabolites were with VIP>1 and p values<0.05.
RNA sequencing data analysis
Transcriptome sequencing and analysis were also done by OE Biotech Co., Ltd. Trimmomatic was used to process raw data (raw reads) (Bolger et al., 2014). Clean reads were obtained by removing poly-N sequences and low-quality reads. Hisat2 was used to map the clean reads and genome (Kim et al., 2015). For transcript-level quantification, bowtie2 was used to analyze the values of the transcript. DESeq (2012) functions were used to identify the DEGs based on SizeFactors and nbinom Test. The significantly differential expression was showed with the thresholds of P < 0.05 and fold change >2 or fold change < 0.5. The transcript expression patterns were executed by hierarchical cluster analysis of DEGs. DEGs were analyzed by GO term enrichment and KEGG pathway enrichment following the hypergeometric distribution (Kanehisa et al., 2008).
16S rRNA metagenetics and bioinformatics analysis
FASTQ format was used for the raw sequence data. Trimmomatic software (Bolger et al., 2014) was used to preprocess the paired-end reads and eliminated ambiguous bases (N). FLASH software was used to assemble the left paired-end reads (Reyon et al., 2012) as follows: overlapping was from 10bp to 200bp, and the mismatch rate was from 0% to 20%. QIIME software was used to denoise and clean the reads (Caporaso et al., 2010) (version 1.8.0). The primer sequences were removed from the clean reads, and Vsearch software was used to execute the OTUs (Rognes et al., 2016). QIIME package was used to select the representative read in every OTU. RDP classifier was used to annotate and blast the usual reads (Wang et al., 2007). All representative reads were annotated and blasted against Unite database (ITSs rDNA) using blast (Lobo et al., 2008).
Functional analysis of integrated metabolomic and transcriptomic data
KEGG database (https://www.kegg.jp/) was used to define the metabolic pathways and integrated metabolomic and transcriptomic datasets. The markers list was introduced into the KEGG database to find out metabolic pathways induced by ENR treatment. The specific metabolite is related to a particular gene if they shared one common KEGG pathway.