Cell culture and transfection
The experimental BMECs were previously constructed in our laboratory[24]. The basal medium consisted of DMEM/F12 (HyClone, Grand Island, NY, USA) supplemented with 10% (v/v) fetal bovine serum (Pasching, Austria) and 1% (v/v) penicillin and streptomycin (PAA, Pasching, Austria). The cell culture plates (Thermo, Suzhou, China) were then incubated in an incubator (Thermo, Marietta, OH, USA) containing 5% CO2 at 37°C. The medium for transfection consisted of DMEM/F12 supplemented with 10% (v/v) fetal bovine serum.
Total RNA extraction
Total RNA from BMECs with high and low milk fat percentages was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. Total RNA was further treated with RNase-free DNase I (NEB, Beijing, China), and RNA longer than 18 nucleotides was recovered. RNA concentration was measured by NanoDrop 2000 (Thermo Scientific, USA). The quality of RNA was assessed by the Agilent2200 (Agilent, USA).
Conjoint analysis of RNA-seq and Solexa sequencing
(1) RNA-seq analysis: The sequencing library of each RNA sample was prepared using Ion Total RNA-Seq Kit v2 according to the protocol provided by manufacturer (Life technologies, USA). The transcript with poly (A)+ RNA of Cattle were analyzed using high-throughput Life technologies Ion Proton Sequencer. Reads sequenced were filtered and mapped to Cattle genome (download from NCBI) using Mapsplice software. The mapped reads was counted to achieve the expression of each gene based on the gene annotation information from NCBI database. In order to achieve the best quality of the RNA Sequencing, Novelbio applied the reads filtration to filter the reads with lower quality and short sequence under following criteria: read length > 50; over 30% base quality >13. Using the internationally recognized algorithm EBSeq based on the fold change (FC) and false discovery rate (FDR) thresholds to screen the DEGs in the BMECs of Chinese Holstein cows with high and low milk fat percentages[25, 26]. Genes were considered differentially expressed according to log2FC >1 or log2FC <-1and FDR <0.05.
(2) Solexa sequencing analysis: We performed a differential expression analysis according to the data of Solexa sequencing obtained by a previous study in our laboratory and screened the DERs in the BMECs of Chinese Holstein cows with high and low milk fat percentages using the internationally recognized EBSeq algorithm based on the FC and FDR thresholds[27]. MiRNAs were considered differentially expressed according to log2FC >1 or log2FC <-1and FDR <0.05.
(3) Conjoint analysis of RNA-seq and Solexa: All target genes regulated by DERs were predicted by the combination of miRbase (http://www.mirbase.org) and TargetScan (http://www.targetscan.org). All the predicted target genes were intersected with the DEGs obtained by RNA-seq, and then the DEGs and DERs were obtained by a combined analysis.
GO and KEGG pathway enrichment analyses
GO analysis was applied to analyze the main functions of the DEGs according to GO, which is the functional classification of NCBI (http://www.ncbi.nlm.nih.gov). The GO category was classified by Fisher's exact test and the χ2 test, and the FDR was calculated to correct the P value: the smaller the FDR was, the smaller the error in judging the P value[28]. We computed the P values for the GO terms of all the DEGs. Enrichment provides a measure of the significance of the function: as the enrichment increases, the corresponding function is more specific, which helps us identify those GO terms with concrete functional descriptions in the experiment.
Pathway analysis was used to determine the significant pathways of the DEGs. Pathway annotations of microarray genes were downloaded from KEGG (http://www.genome.jp/Wegg/). Fisher’s exact test was used to find the significant enrichment pathways. The resulting P values were adjusted using the BH FDR algorithm (Benjamini and Hochberg 1995). Pathway categories with an FDR < 0.05 were reported.
Vector construction of candidate miRNAs and luciferase reporter
(1)Construction vector of six candidate DERs: The mimics, inhibitor and negative control vectors of 6 miRNAs (bta-miR-21-3p, bta-miR-124a, bta-miR-148a, bta-miR-877, bta-miR-2382-5p and bta-miR-2425-5p) were constructed by GenePharma Company (Suzhou, China).
(2)Construction of luciferase reporter vectors: We predicted the DNA-binding sites of miRNAs to the target gene’s 3'-UTR with TargetScan. For the luciferase reporter assay, the WT and mut DNA sequences were obtained by annealing performed in a 10 μL reaction with 0.5 μL forward primer, 0.5 μL reverse primer, 1 μL 10 × Ex Taq buffer and 8 μL ddH2O under the following procedure: 95°C for 5 minutes, then dropped at a rate of 1.5°C per minute to 22°C. The DNA fragment was then ligated with the linearized pmiR-BR-REPORT plasmid to construct the WT and mut vectors. All DNA oligo nucleotides used for annealing were synthesized by Sangon (Shanghai, China), and the sequence details are shown in Table 1.
Determination of triglyceride content
The Triglyceride Extraction Kit (APPLYGEN, China) was used to separate the cells following the manufacturer's instructions. Then, the cell lysates were collected in 1.5 mL centrifuge tubes, heated at 70°C for 10 minutes, centrifuged at 2000 g for 5 minutes at room temperature, and the supernatant was extracted to detect the contents of triglyceride and protein on a SpectraMax M5 Microplate Reader. The protein content was used as a reference to correct the variation in cell number, each sample was analyzed in triplicate.
Quantitative RT-PCR of DEGs and DERs
Total RNA was extracted and reverse transcribed to cDNA using the PrimeScriptTM RT Reagent Kit with gDNA Eraser (Takara Biological Company, Japan) for the quantitative RT-PCR analysis of 8 DEGs (GOT2, SH3BP5, PRDX6, ILDR2, PDE4D, PEG10, SOD3 and ADAMTS1) and 6 DERs (bta-miR-2382-5p, bta-miR-124a, bta-miR-21-3p, bta-miR-148a, bta-miR-877 and bta-miR-2425-5p). Quantitative RT-PCR was performed in a 10 µL reaction with 0.5 µL forward primer, 0.5 µL reverse primer, 5 µL SYBR Green Real-Time PCR Master Mix (Roche, Germany), 1 µL cDNA and 3 µL ddH2O according to the manufacturer's instructions under the following procedure: 95°C for 30 s, 40 cycles of 95°C for 5 s and 60°C for 30 s on a PCRmax (Eco, United Kingdom). Each sample was analyzed in triplicate, and the average threshold cycle (Ct) was calculated with the 2-△△Ct method. The reverse transcription primers of the miRNAs and primers for the quantitative RT-PCR analysis of DERs and target DEGs were designed using Primer 6.0 and synthesized by Sangon (Table 2).
Luciferase reporter assay
The luciferase reporter experiments were performed with BMECs at a density of 1 × 105 cells in 24-well plates (BIOFIL, Guangzhou, China). To validate the target relationship of the candidate miRNAs and genes, 1.5 µg miRNA mimic plasmid was cotransfected with 1.5 µg WT, 1.5 µg mut or 1.5 µg pmiR-BR-REPORT plasmid into BMECs. Forty-eight hours after transfection, the firefly (hluc+) and Renilla (hRluc) luciferase activities were measured with the Dual-Glo luciferase assay system (Promega) according to the manufacturer's instructions. The hluc+ gene was used as a reference gene to correct the variation in transfection efficiency, and the relative luciferase activity was calculated as hRluc/hluc+. All experiments for the luciferase reporter assay were performed in triplicate.
Statistical analysis
The data were analyzed by GraphPad 6.0 statistical software with a completely random design of Student's t test for the comparison of two independent groups. The mean and standard deviation (SD) or standard error of mean (SEM) among triplicate samples were calculated, and p<0.05 was considered to indicate statistical significance.