Animals and experimental design
Experimental procedures and animal care were conducted in compliance with the European Communities Council Directive 2010/63/EU. The experiment was conducted at the Ducks and Goose Experimental Facility – INRAE, UEPFG, (Benquet, France) that received the accreditation number B40–037-1. The protocol and procedures were approved by the French Minister of Higher Education, Research and Innovation (authorization APAFIS#1847-2015092213418825v2).
The experimental design has already been described [18]. Briefly, 60 female common ducks received adequate level of Met until the age of 10 weeks and were then divided into two groups. They were fed experimental diets from 10 to 51 weeks of age. Met is the first main limiting amino acid in a typical corn−soybean diet for poultry [36]. In laying ducks, the recommendations vary from 0.40 to 0.45% [37,38,39,40]. Thus, two levels of total Met were used for the experimental diets: 0.25% for Met-restricted diets (R Group) and 0.40% for control diets (C Group) that meets Met requirements for female laying ducks (Additional Table 2). The mule duckling production was performed by 2 artificial inseminations per week between 34 and 36 weeks of age, with the semen of 15 Muscovy drakes not subjected to any dietary treatment and fed commercial diets. The eggs were incubated for 28 days at 37.6 °C and 60% average relative humidity throughout the entire incubation period (incubator Sologne, La Nationale, Briaire, France). They were then put in a hatcher (hatcher Bretagne, La Nationale, Briaire, France) for 4 days at 37.3 °C and 80% average relative humidity. Ducklings traits were recorded at hatching. The ducklings that were the offspring of the females from the R and C groups are subsequently assigned to R and C groups, respectively.
Phenotypic traits of ducklings were recorded at hatching on 180 and 190 ducklings from the R and C groups respectively, as reported in Bodin et al., 2019 [18]. Moreover, a total of 58 ducklings were sacrificed by cervical dislocation at hatching (12 females and 16 males from the C group and 15 females and 15 males from the R group). These ducklings did not receive any feed before being sacrificed. Their liver weight was recorded. In some cases, the gallbladder has been damaged and for this reason, only the livers from 8 females and 13 males from the C group and 15 females and 15 males from the R group, that were satisfactorily retrieved were immediately immersed in liquid nitrogen before being transferred to a − 80 °C freezer.
RNA extraction and reverse transcription
Frozen liver samples from newly hatched ducklings of both sexes and of both control and restricted maternal diet groups (10 males and 8 females in the C group and 10 males and 10 females in the R group) were ground using a Retsch grinder at 30 Hz for 45 seconds, in liquid nitrogen. Next, 80 to 100 mg of tissue powder were processed as previously described [41] for RNA extraction and purification using the TRIzol® method (Invitrogen, California, USA) followed by a column from Nucleospin RNA kit (Macherey Nagel, France) and following the manufacturer’s instructions. The on-column DNAse treatment was done with 20 μl of rDNAse (Macherey Nagel) and 80 μl of reaction buffer for 20 min to avoid DNA contamination as recommended [42, 43]. The total RNA was quantified using a NanoDrop 8000 spectrophotometer (Thermo Fisher, Illkirch, France) and stored at − 80 °C. Its integrity was controlled by electrophoresis and using an Agilent 2100 Bioanalyzer, with the RNA 6000 Nano Lab Chip Kit (Agilent Technologies, Massy, France). Reverse transcription was carried out immediately after the quality control evaluation, and the same amount of total RNA was used for all experimental samples, in accordance with recommendations [44]. The reaction used the RNase H-MMLV reverse transcriptase SuperScriptTM II (Invitrogen, California, USA), RNasin® Ribonuclease inhibitor (Promega Corporation, USA) and oligo (dT)15 (Sigma Aldrich, France). The cDNAs were then diluted in RNase free water and stored at − 80 °C.
Primer design and qPCR validation
The study targeted 100 genes selected from the literature for being related to energy metabolism as well as to amino acid transport, oxidative stress, apoptotic activity and susceptibility to liver damage. Sequences were obtained either from Anas platyrhynchos when available, or from Gallus gallus on NCBI [45] and/or Ensembl [46] databases. The two primers used for each gene (Additional Table 3) were designed in exons, but each on either side of an intron, and for an annealing temperature of 60 °C, using either Primer3Plus [47] or LightCycler® Probe Design Software 2.0 (Roche Applied Science). The primer sequences were blasted to the databases to confirm that they were specific to the gene in question. PCR products were first subjected to 2% agarose gel electrophoresis to confirm amplicon size. Next, the primer pairs showing a specific band and no primer dimers were selected to be tested by qPCR, using SYBR green fluorescence detection (Applied Biosystems) and a QuantStudio6 (Thermo Fisher Scientific). Each primer pair was tested on four serial dilutions of a pool of cDNA (cDNA of all the animals used in the study) to obtain a standard curve and check the PCR efficiency, each point being done in duplicate. The conditions were: 50 °C for 2 min, denaturation at 95 °C for 10 min, 40 cycles of 15 s at 95 °C and 1 min at 60 °C. A gradual temperature increasing from 60 °C to 95 °C was added to analyze the melting curves and detect primer dimers. The Cq values and the PCR efficiency were obtained directly from the QuantStudio Real-Time PCR software v1.3.
Identification of potential reference genes
The expressions of 9 potential reference genes were tested for their stability in livers of newly hatched mule ducklings. Primer sequences were either from Chapman et al., 2016 [48] (designed in Anas platyrhynchos: ALB, GAPDH, NDUFA10 and RPS13), from Staines et al., 2016 [49] (designed in Gallus gallus: HMBS and TBP) or newly designed in Anas platyrhynchos (HPRT1, POLA1 and TUBA1C). These 9 primer pairs were tested by qPCR on liver cDNA from 8 ducklings of both sexes and of both maternal diet groups (C group and R group) with SYBR green fluorescence detection, on a QuantStudio6. The selection of the 5 most stable genes (GAPDH, HMBS, NDUFA10, RPS13, TBP), highlighted in dark grey in Additional Table 3) was done with the package SLqPCR on RStudio [50].
Quantitative PCR and gene expression analysis
The quantification of gene expression was performed using the 96.96 Dynamic Array integrated fluidic circuits (IFCs) and the BioMark HD system from Fluidigm as previously described [51].
In this paper we present the study of transcripts of genes involved in energy metabolism in the liver of newly hatched ducklings but the full experiment was conducted on 168 liver samples (38 samples of newly hatched duckling livers and 130 samples of older duck livers) and a total of 170 genes, either targeting energy metabolism for this study (100 genes) or playing a role in one-carbon metabolism (70 genes), as well as five reference genes identified in the previous step. As the technology used did not allow all samples and genes to be analyzed at the same time, care was taken to randomize the samples onto two arrays and the genes into each specific target amplification (STA). Thus, a total of four chips were made. In this study, we focus our interest on the 38 liver samples from newly hatched ducklings and 100 genes involved in energy metabolism.
For each chip, the first step consisted of a STA with 14 cycles performed on the cDNA samples, a calibrator sample (a pool of the 38 cDNA samples), a cDNA pool of the 168 cDNA samples in five-fold serial dilutions (to determine the PCR amplification efficiency), a duck genomic DNA control, an internal control (human genomic DNA) and a negative control (TE). The first STA was performed with a pre-amplification primer mix that contained primers for the five potential reference genes and 85 target genes (among those 53 are some of this study). A second STA contained primers for the five potential reference genes and other 85 target genes (among those 47 are some of this study). The resulting STA cDNA samples were treated according to Bonnet et al., 2013 [51].
Data were analyzed with the Fluidigm Digital PCR Analysis software (version 4) using the linear (derivative) baseline correction method and the auto (global) cycle threshold (Cq) method. The data were pre-processed before expression analysis. Indeed, the cycle threshold values (Ct) recorded from amplifications whose melting curves showed either abnormal Tm (melting temperature) or double peaks (corresponding to a mixture of expected and aberrant PCR products) or a high baseline, were removed. To measure the efficiency (E) of PCR for each gene, the slope of the standard curve obtained with serial dilutions of the 168 cDNA pool was used and genes with less than three dilution points were eliminated. Finally, all genes with an efficiency greater than 2.2 or less than 1.7 were removed from the analysis.
The relative expression (RE(i,j)) for each gene (i) and each sample (j) was calculated as proposed by Pfaffl [52]: RE(i,j) = Eff(i)^(Cq(i,cal) – Cq(i,j)) where Cq(i,cal) is the Cq for the gene (i) determined for the calibrator sample (a pool of the 38 cDNA samples). To determine the reference genes, the stability of the five potential genes was tested using the GeNorm (version 3.4) algorithm [50]. The three most stable genes were chosen as references genes (GAPDH, RPS13 and TBP). Finally, the relative expression (RE(i,j)) for each gene (i) and each sample (j) was normalized with the geometric average expression of the three most stable reference genes as proposed by Vandesompele et al. [50].
In the end, 13 of the 100 genes studied and 2 of the 38 liver cDNA samples showed more that 25% of missing data and were removed from the study. Moreover, another cDNA sample was also removed from the data set because it showed data points that differed significantly from other observations. The 13 removed genes are highlighted in light grey in Additional Table 3. The 35 remaining liver cDNA samples are from 9 male and 8 female ducklings from the C group and 10 male and 8 female ducklings from the R group.
Statistical analyses
For the 87 remaining genes, the few missing values were imputed inside each group with similar sex and maternal diet using the function imputPCA with 3 principal components of the missMDA package of the R software [53]. These normalized and imputed relative expressions were then transformed using the function qqnorm(Y)$x to make the data follow a centered reduced normal distribution. These transformed data were then used to describe the whole dataset.
First a heatmap was obtained by using the R gplots package thus defining groups of animals and genes with similar expression patterns and a Principal Component Analysis (PCA) was performed with the package MixOmics of R software [54] on the 87 genes. The individuals were plotted on the two first principal components of the PCA score plot.
Then, ANOVAs were conducted on the transformed normalized relative expressions of the 87 genes using a linear mixed model fitted with the ASReml software [55] that included the maternal diet, the sex of the duckling, and the interaction between them as fixed effects, as well as the duckling associated to its relationship matrix as a random effect. Genes with a significant difference - diet P-value < 0.05 assessed after a Benjamini-Hochberg (1995) [56] correction, in order to account for multiple tests and named (diet P-Value (BH)) - were selected and considered as differentially expressed genes (DEGs). The effect of sex on gene expression was also evaluated (sex P-value (BH)), as well as the interaction between sex and diet effects (sex*diet P-value (BH)). For each gene, least square means and standard deviations were calculated for the 2 groups of maternal diet (R group and C group), for the 2 sexes, and finally for the four subgroups of interest, i.e., males in the R group and C group and females in the R group and C group (Additional Table 1).
Further analyses were performed on the 16 DEGs for the diet effect. Another PCA was performed on the 16 DEGs with the package MixOmics of R software [54], using the qqnorm transformed normalized relative expressions. The biplot of variables and samples from the two first principal components was obtained with the package Factoextra of R. Concentration ellipses were plotted around each group mean points with a confidence level of 0.75 using the package ellipse.
Then, correlation matrices were performed using phenotypic traits of the ducklings and the 16 DEGs for the diet effect, using the imputed normalized relative expressions to be consistent with the phenotypic data which were also not transformed.
The correlation matrices were plotted with the package Hmisc of R using the functions rcorr and corrplots [57, 58] and only the correlations with a P-value < 0.05 were reported.
Lastly, this study was conducted and is reported in accordance with the ARRIVE guidelines [59].
