Animals and samples preparation
Animal work was conducted according to the guidelines for the care and use of experimental animals established by the Ministry of Science and Technology of the People’s Republic of China (approval number: 2006-398), and was approved by the Laboratory Animal Management Committee of Northeast Agricultural University (Harbin, China). The experimental birds were obtained from the Avian Farm of Northeast Agricultural University (Harbin, Heilongjiang, China). These broilers under divergent selection over 19 generations were employed from Northeast Agricultural University broiler lines divergently selected for high and low abdominal fat content (NEAUHLF), exhibiting a large difference in abdominal fat content as previously described [3]. In total, ten male birds (lean line, n = 5, and fat line, n = 5) at 7 weeks of age from the 19th generation of NEAUHLF were used for RNA-seq and iTRAQ analysis, and these birds were kept under the same environmental conditions and had free access to feed and water. Abdominal fat tissues were collected right after these birds were euthanized by intramuscular injection of pentobarbital (Sigma, St. Louis, MO, USA) (40mg/kg) under deep anesthesia, and then immediately frozen in liquid nitrogen and stored at -80℃. The detailed information of selected chickens’ body weights, abdominal fat weights (AFW) and abdominal fat percentages (AFP) were listed in Table S1.
Transcriptomic data collection and analysis
Total RNA from abdominal fat tissues was extracted using the TRIzol reagent (Invitrogen, New Jersey, NJ, USA), and genomic DNA was removed by DnaseI treatment. RNA purity, concentration and integrity were checked by NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA), and RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA), respectively. After removal of ribosomal RNA and cleaning-up of rRNA free residue by a Ribo-Zero™ rRNA Removal Kit (Epicentre, USA), the sequencing libraries were generated using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations. cDNA fragments of 150-200 bp in length were selected and purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then, library quality was assessed by the Agilent Bioanalyzer 2100 system. Finally, after cluster generation (cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS, Illumina), the libraries were sequenced on an Illumina Hiseq 4000 platform, and 150 bp paired-end reads were produced.
After demultiplex and quality filtering of raw data, clean reads were obtained and aligned to the G. gallus 5.0 reference genome assembly using HISAT2 (v.2.0.4) [11]. The mapped reads of each sample were assembled by StringTie (v1.3.1) [12] in a reference-based approach. The software Cuffdiff (v2.1.1) [13], which provides statistical routines for determining gene expression data using a model based on the negative binomial distribution, was used to calculate FPKMs (fragments per kilobase of exon per million fragments mapped). Transcripts with a p-value ≤ 0.01 and fold changes ≥ 1.5 or ≤ 0.67 were assigned as significantly differentially expressed.
Proteomics
Protein was extracted according to Damerval et al [14], checked by SDS-PAGE (Fig.S1) and determined concentration by the Bradford method [15]. Following reduction, cysteine alkylation, and trypsin digestion, total proteins were treated to obtain peptides, and labeled with iTRAQ 8-plex or iTRAQ 4-plex reagents (AB SCIEX, USA), as 113 (LL1), 114 (LL2), 115 (LL3), 116 (LL4), 117 (LL5), 118 (FL1), 119 (FL2), and as 117 (FL3), 118 (FL4), 119 (FL5), respectively. We pooled all samples and labeled as 121 in 8-plex and 4-plex iTRAQ, to calibrate the two iTRAQ experimental data sets. Then, the iTRAQ-labeled peptide mixture was reconstituted and loaded on SCX (strong cation exchange) column, which were subjected to nanoelectrospray ionization, followed by tandem mass spectrometry (MS/MS) in a TripleTOF 5600 system (AB SCIEX, USA).
The MS/MS data were processed with ProteinPilot Software v. 5.0 (AB SCIEX, USA) against Gallus gallus database using the Paragon algorithm [16]. The experimental data from tandem mass spectrometry (MS) were utilized to match the theoretical data to identify proteins, which was performed by the search option (with an emphasis on biological modifications). An automatic decoy database search strategy was used to estimate the false discovery rate (FDR) calculated as the false positive matches divided by total matches, using the PSPEP (Proteomics System Performance Evaluation Pipeline Software, integrated into the ProteinPilot Software). The significantly differentially expressed proteins were identified using the following criteria: 1) peptide groups considered for quantification required at least 2 peptides, and a global FDR less than 1% was used; and 2) a fold change ≥ 1.5 or ≤ 0.67 and with p-value ≤ 0.05 (t-test).
RT-qPCR analysis
To validate RNA-Seq results, 20 DEGs with higher expression levels and larger fold changes were validated by RT-qPCR. Ten male birds (n=5 for each line) from the 19th generation of NEAUHLF were used. Total RNA from abdominal fat tissue was reversely transcribed into cDNA using a PrimeScript™ RT Reagent Kit (Takara, Dalian, China). FastStart Universal SYBR Green Master kit (Roche) and the ABI 7900 PCR detection system were used to perform RT-qPCR. The program began at 95 ℃ for 30 s for activation, followed by 40 cycles of amplification at 95 ℃ for 5 s and 58 ℃ for 30 s. An additional 15 s at 95 ℃, 1 min at 60 ℃, and 15 s at 95 ℃ were performed for the melt curve stage. The housekeeping gene TATA-Box binding protein (TBP) was used as the control. RT-qPCR primer pairs were designed by Primer Premier 6.0 and the detailed information were listed in Table S2. The comparative 2-ΔΔCt method was used to determine the statistical significance.
PRM-MS analysis
To verify the protein expression level obtained by iTRAQ analysis, 10 DEPs with higher expression levels and larger fold changes were selected for validation. Signature peptides for the target proteins were defined according to the iTRAQ data, and only were unique peptide sequences selected for the PRM-MS analysis (Table S3). We randomly selected 6 male birds from the 19th generation (n=3 for each line) of NEAUHLF, and extracted the proteins from abdominal fat tissues, which were prepared following the iTRAQ protocol. Tryptic peptides were loaded on C18 stage tips for desalting prior to reversed-phase chromatography on an Easy nLC-1200 system (Thermo Scientific).
One-hour liquid chromatography gradients with acetonitrile ranging from 5 to 35% were used, and PRM was performed on a Q-Exactive Plus mass spectrometer (Thermo Scientific). Methods were optimized for collision energy, charge state, and retention times for the most significantly regulated peptides experimentally, using unique peptides of high intensity and confidence for each target protein. The mass spectrometer was operated in positive ion mode and with the following parameters: the full MS1 scan was acquired with the resolution of 70000 (at 200 m/z), automatic gain control (AGC) target values 3.0×106, and a 250 ms maximum ion injection times. Full MS scans were followed by 20 PRM scans at 35000 resolution (at m/z 200) with AGC 3.0×106 and maximum injection time 200 ms. The targeted peptides were isolated with a 2.0 Th window and fragmented at a normalized collision energy of 27 in a higher energy dissociation (HCD) collision cell. The raw data were analyzed using Skyline (MacCoss Lab, University of Washington) [17] to get the signal intensities of individual peptide sequences.
For PRM MS data, each sample’s average base peak intensity was extracted from the full scan acquisition using RawMeat (version 2.1, VAST Scientific). The normalization factor for sample N was calculated as fN = the average base peak intensity of sample N divided by the median of average base peak intensities of all samples. The area under curve (AUC) of each transition from sample N was multiplied by this factor. After normalization, the AUC of each transition was summed to obtain AUCs at the peptide level. Relative protein abundance was defined as the intensity of a certain peptide.
Gene enrichment analysis
DEGs and DEPs were submitted for the gene ontology (GO) analysis by DAVID (Database for Annotation, Visualization and Integrated Discovery, https://david.ncifcrf.gov/) version 6.8, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (http://kobas.cbi.pku.edu.cn/kobas3). The thresholds for significant enrichment was set at p-value < 0.05.
Statistical analysis
Chicken body measurement data were shown as mean ± SD. Student’s t-test was used to compare the differences between two groups, and the threshold of significance was set at p < 0.05.