Animals and sampling
All animals used in the experiments were female Lushi blue-shell-egg chickens obtained from the Animal Center of Henan Agricultural University. The chickens were housed in cages separately under the same environmental conditions with ad libitum access to food and water. The corn soybean basal diets containing 15% crude protein, 5% crude fat, 6% Crude Fiber, and 2750 kcal/kg energy were obtained from Shandong Newhope Liuhe Group Co., Ltd. (Liuhe, Shandong, China). A total of 30 chickens with similar body weight (680 ± 24.6 g) at 10 weeks of age were selected and subdivided randomly into four groups. Three groups (n = 8) were injected intramuscularly with 17β-estradiol (dissolved in olive oil) (Sigma, St. Louis, MO, USA) at 0.5, 2.0 or 8.0 mg / kg of body weight. One group (n = 6) serving as a control was given vehicle (olive oil) only. All chickens were euthanized by cervical dislocation at 12 h after treatment. A portion of the fresh liver tissue was embedded with optimum cutting temperature compound for oil red O staining. The rest of the liver was snap-frozen in liquid nitrogen and stored at − 80 °C in a freezer until use. Three livers from the control group and three livers from the 8.0 mg/kg 17β-estradiol group were subjected to RNA-Seq. Livers from the 8.0 mg/kg 17β-estradiol group were also used for ChIP-Seq.
Primary hepatocytes culture and estrogen administration
The isolation, purification, and culture of chicken primary hepatocytes were performed as described previously [19, 42]. The chicken hepatocytes were isolated from 18-day chicken embryonic livers. Cells were divided into four groups with triplicates. The cells were starved for 6 h after they grew to 80% confluence. Then, we added 17β-estradiol dissolved in 0.1% ethanol to final concentrations of 25 nM, 50 nM, or 100 nM. Cells treated with ethanol only were used as controls. After 12 h incubation, cells were collected and stored at − 80 °C until use.
Oil red O staining
The optimum cutting temperature compound-embedded liver tissue sample was sectioned with a freezing microtome. The sections were washed with PBS, fixed with 4% paraformaldehyde for 30 min at room temperature and then stained with 0.2% oil red O solution (Sigma) for 15 min. After staining, sections were washed with PBS and counterstained with hematoxylin for 5 min. The washed sections were observed under a microscope at 200× and 400× magnifications.
Total RNA was extracted from the chicken liver tissues using the TRIzol® reagent following manufacturer instructions (Invitrogen, Carlsbad, CA, USA). Libraries were constructed using the TruSeq Stranded Total RNA LT (with Ribo-Zero™ Gold) Set B (cat. # RS-122-2302; Illumina, San Diego, CA, USA). The libraries were sequenced on the Illumina HiSeq 4000 platform via the 150 bp pair-end sequencing strategy following manufacturer instructions. The obtained raw reads were cleaned using the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html). FASTQ files of clean paired-end reads were aligned to the reference genome using TopHat. The reference genome for the chicken (Gallus 5.0) was downloaded from the UCSC Genome Browser. The transcript abundance and putative novel mRNA isoforms were analyzed in the Cufflinks software. Next, the FPKM values were used to quantify the gene expression levels. In this study, the FDR was used to determine the threshold of the p value in multiple tests and analyses.
After annotation, the unannotated transcripts were employed to identify lncRNAs. We applied the following filtering criteria: (1) exon number ≥ 2; (2) transcript length ≥ 200 bp; (3) predicted open reading frame < 300 bp; (4) the transcript aligned to Pfam  without significant hits; (5) coding–noncoding index score < 0 ; and (6) coding potential calculator score < 0 .
Small RNA was extracted from chicken liver tissue samples using the mirVana™ miRNA Isolation Kit (cat. # AM1561, Ambion, Austin, TX, USA). The extracted RNA samples were ligated sequentially with 3′ and 5′ RNA adapters by means of T4 RNA ligase (cat. # M0242l; BioLabs, Beverly, MA, USA). The ligated RNA samples were reverse-transcribed and amplified by PCR to generate cDNAs. The cDNAs of appropriate lengths were purified from an agarose gel to construct sequence libraries. After purification, the small-RNA libraries were quantified on a Qubit Fluorometer (Invitrogen) and used for cluster generation and 50 bp single-end sequencing analysis with an Illumina HiSeq 2000 system.
An initial filtering step was performed using the Fastx-toolkit to remove adaptor sequences, low-quality reads (base quality less than 10), and short reads (shorter than 18 nt) before the clean reads were collected and subjected to bioinformatics analysis. Clean reads were first compared with the miRBase database 21.0 using the CLC Genomics Workbench 5.5 commercial software (CLC Bio, Aarhus, Denmark) to identify known miRNAs. The unmatched sequences were screened against the noncoding RNA database, Rfam, and piRNA database to filter out rRNAs, tRNAs, snRNAs, and snoRNAs. After the elimination of repeat-associated small RNAs, degradation fragments of mRNAs, and known miRNAs, the remaining reads were mapped to the chicken genome using Burrows–Wheeler Alignment to obtain pre-miRNA sequences . The mapped sequences were then utilized to predict novel miRNAs by means of MIREAP. The expression of the identified miRNAs was normalized by calculation of TPM. The differential significance was identified using EdgeR software with the following thresholds: FDR < 0.05 and fold change ≥1.5 or fold change ≤0.667.
The ChIP assay was performed as described previously . In brief, fresh liver tissue samples collected from chickens treated with 17β-estradiol were cross-linked with 1% formaldehyde, after which the cross-linking was stopped by the addition of 2.0 M glycine. Isolated DNA was fragmented and incubated with a monoclonal specific ERα antibody (0.2 mg/mL, cat. # MA5–13065; Invitrogen) or IgG antibody (Cell Signaling Technology, Danvers, MA, USA). The immunoprecipitated DNA fragments were quantified and used to construct the DNA libraries using the ChIP-Seq Sample Prep Kit (Illumina). The libraries were then sequenced on the Illumina HiSeq 2500 platform using 50 bp single-end sequencing analysis. Raw reads were subjected to a quality check and trimming, after which the obtained clean reads were aligned to the chicken genome (Gallus 5.0) via Burrows–Wheeler Alignment. Peak calling procedures were performed in MACS2  with threshold p values ≤0.005. To perform a motif search, all positive peaks were used, and the peak sequences were extended on both sides to obtain a 500 bp sequence. The extended 500 bp sequences were then used to discover motifs using MEME and FIMO. MAST was employed for motif alignment to exclude duplicate motifs, whereas Tomtom was used for annotation of the discovered motifs.
GO and Kyoto Encyclopedia of KEGG pathway enrichment analyses of genes were performed by DAVID. p value < 0.05 was used as the cut-off criterion for GO and KEGG pathway enrichment analyses.
miRNA target gene prediction
The miRNA target prediction software programs miRDB, TargetScan 7.1 and PicTar were used to predict miRNA target genes. Only when miRNA-mRNA pairs were positively predicted by > 2 of these software programs were the pairs taken as a positive result.
A target fragment was amplified by PCR with special primers containing XhoI and NotI restriction sites. Seven bases of binding sites in 3’UTR were deleted to create a mutant version by means of special primers designed for overlap extension PCR. Then, the PCR product was cloned into the XhoI-and-NotI double-digested psi-CHECK™-2 plasmid (Promega, Madison, WI, USA) by means of T4 ligase (Biolabs, Beverly, MA, USA). All of the constructed vectors were confirmed by PCR and sequencing (BGI, Shenzhen, China). A chicken embryonic fibroblast cell line (DF1) was cultured and cotransfected with the constructed luciferase vector and a miRNA mimic or mimic negative control (mimic NC) with Lipofectamine™ 2000 (Thermo, Waltham, MA, USA). The cell lysates were harvested 48 h after transfection. The Renilla luciferase and firefly luciferase activities were measured using a Dual Luciferase Reporter Assay System (Promega). For each transfected group, the procedure was performed in triplicate in at least three independent experiments.
Quantitative real-time PCR (qRT-PCR)
The expression levels of some selected mRNAs and lncRNAs were validated by qRT-PCR. The PrimeScript™ RT Reagent kit with gDNA Eraser (TaKaRa, Dalian, China) was used to synthesize the cDNA according to manufacturer instructions. The qRT-PCR was performed using the SYBR Green method in a LightCycler® 96 instrument. Each reaction contained 5 μL of SYBR Green PCR Master Mix (TaKaRa), 3.5 μL of RNase-free water, 0.5 μL each of forward and reverse primers, and 0.5 μL of extracted cDNA. The reactions were amplified using the following conditions: denaturation at 95 °C for 5 min; followed by 40 PCR cycles at 95 °C for 30 s, 60 °C for 30 s, and 72 °C for 20 s; and then a further 10-min extension at 72 °C. All reactions were performed in triplicate. The expression levels were measured in terms of the cycle threshold (Ct) and then normalized to the expression of β-actin using the 2−△△Ct method. The primers used were designed using the NCBI Primer-BLAST tool and synthesized by Sangon Biotech (Shanghai, China). All primer sequences are listed in Table S6.
Expression of the miRNA was detected by stem-loop qRT-PCR. Reverse transcription of miRNAs was performed using miRNA-specific stem-loop primers and the PrimeScript RT Reagent Kit (TaKaRa). The primers used for reverse transcription and qRT-PCR were designed and purchased from GenePharma Co., Ltd. (Shanghai, China). The miRNA expression levels were normalized to the expression of U6, and the other protocols were the same as described above.
ChIP-quantitative PCR (ChIP-qPCR) analysis
The immunoprecipitated DNA fragments were also used to verify the interaction of ERα and selected binding sites by ChIP-qPCR. Meanwhile, 20% of starting chromatin without chromatin immunoprecipitation served as input to represent the unselected DNA content. The fold enrichment method was chosen to normalize the ChIP-qPCR data: Fold enrichment = log2−ΔΔCt, ΔCt = Ct (IP) − Ct (Input) − log25, ΔΔCt = ΔCt − ΔCt (IgG). Gene-specific primers for the putative binding-site regions were designed using the NCBI Primer-BLAST tool and synthesized by Sangon Biotech (Shanghai, China). All ChIP-qPCR primer sequences used are listed in Table S7.
Statistical analyses were carried out using SPSS version 20.0 (IBM, Chicago, IL, USA). One-way ANOVA were used for statistical analysis, followed by Dunnett’s test. The results were presented as Mean ± SEM of more than 6 replicates, p < 0.05 was considered statistically significant.