extfieldsofstudy: Medicine Chemistry

s2fieldsofstudy: Biology Chemistry

text:
Ethanolic Extracts from Azadirachta indica Leaves Modulate Transcriptional Levels of Hormone Receptor Variant in Breast Cancer Cell Lines

Breast Cancer (BC) encompasses numerous entities with different biological and behavioral characteristics, favored by tumor molecular complexity. Azadirachta indica (neem) presents phenolic compounds, indicating its potential as an antineoplastic compound. The present study aimed to evaluate the cellular response of MCF10, MCF7, and MDA-MB-231 breast cell lines to ethanolic extracts of neem leaves (EENL) obtained by dichloromethane (DCM) and ethyl acetate (EA) solvent. Extracts’ antiproliferative activities were evaluated against MCF 10A, MCF7, and MDA-MB-231 for 24 and 48 h using MTT assay. ESR1, ESR2, AR, AR-V1, AR-V4, and AR-V7 transcripts were quantified through qPCR for 0.03125 μg/mL of DCM and 1.0 μg/mL for EA for 48 h. The EENL was tested on Drosophila melanogaster as a sole treatment and then also together with doxorubicin. Antiproliferative effect on tumor cell lines without affecting MCF 10A were 1.0 µg/mL (P < 0.001) for EA, and 0.03125 µg/mL (P < 0.0001) for DCM, both after 48 h. Transcriptional levels of AR-V7 increased after treatment. In vivo assays demonstrated that EENL induced fewer tumors at a higher concentration with doxorubicin (DXR). The behavior of AR-V7 in the MDA-MB-231 tumor lineage indicates new pathways involved in tumor biology and this may have therapeutic value for cancer.

Introduction
According to the World Health Organization (WHO) [1], Breast Cancer (BC) is the most common malignant tumor among women. An increase from 14 million in 2012 to 22 million new cases is estimated by 2022. By 2035, the number of deaths is expected to grow by 70% [2]. BC is a disease that encompasses numerous entities, with peculiar biological and behavioral characteristics favored by a complex molecular microenvironment [3]. The treatment with EENL-DCM extract for 24 h reduced the proliferation of MCF7 compared to MCF 10A (P < 0.05) at 0.015625 μg/mL ( Figure 1C). After 48 h ( Figure 1D), DCM extract inhibited the proliferation of the triple-negative tumor cell line (MDA-MB-231), with extracts' concentrations ranging from 0.0625 to 0.25 μg/mL, when compared to MCF7 lineage. Compared to MCF 10A, the ESR + tumor cell line (MCF7) showed a reduction in viability at a concentration of 1.0 μg/mL. Only at the concentration 0.03125 μg/mL and 0.0625 were there significant differences (P < 0.001) between the non-tumorigenic lineage (MCF 10A) and the triple-negative cells' (MDA-MB-231) viability.
Based on these results, the concentrations that showed antiproliferative effect on tumor cell lines, but did not decrease MCF 10A viability, were chosen for molecular assays. For EENL-EA extract we used 1.0 μg/mL (P < 0.001), and for EENL-DCM extract 0.03125 μg/mL (P < 0.0001) after 48 h of treatment, as demonstrated in the time-course Figure 2. The treatment with EENL-DCM extract for 24 h reduced the proliferation of MCF7 compared to MCF 10A (P < 0.05) at 0.015625 µg/mL ( Figure 1C). After 48 h ( Figure 1D), DCM extract inhibited the proliferation of the triple-negative tumor cell line (MDA-MB-231), with extracts' concentrations ranging from 0.0625 to 0.25 µg/mL, when compared to MCF7 lineage. Compared to MCF 10A, the ESR + tumor cell line (MCF7) showed a reduction in viability at a concentration of 1.0 µg/mL. Only at the concentration 0.03125 µg/mL and 0.0625 were there significant differences (P < 0.001) between the non-tumorigenic lineage (MCF 10A) and the triple-negative cells' (MDA-MB-231) viability.
Based on these results, the concentrations that showed antiproliferative effect on tumor cell lines, but did not decrease MCF 10A viability, were chosen for molecular assays. For EENL-EA extract we used 1.0 µg/mL (P < 0.001), and for EENL-DCM extract 0.03125 µg/mL (P < 0.0001) after 48 h of treatment, as demonstrated in the time-course Figure 2.

Transcriptional Profile of Hormone Receptors after Treatment with EENL
Gene expression of ESR1, ESR2, AR, AR-V1, AR-V4, and AR-V7 in breast cells was evaluated before and after 48 h of treatment with 1.0 μg/mL of EENL-EA extract ( Figure 3). AR-V7 expression increased 2.85-fold (P < 0.01) in treated MDA-MB-231 cells after 48 h ( Figure 3H). The remaining genes in MDA-MB-231, and all genes in MCF7, did not display statistically significant differences in gene expression levels. Although not significant, the mean AR, AR-V1, ESR1, and ESR2 relative gene expression levels were higher in 48 h-treated MDA-MB-231 cells (31.74, 2.03, 7.07, and 4.37-fold, respectively).

Transcriptional Profile of Hormone Receptors after Treatment with EENL
Gene expression of ESR1, ESR2, AR, AR-V1, AR-V4, and AR-V7 in breast cells was evaluated before and after 48 h of treatment with 1.0 µg/mL of EENL-EA extract ( Figure 3). AR-V7 expression increased 2.85-fold (P < 0.01) in treated MDA-MB-231 cells after 48 h ( Figure 3H). The remaining genes in MDA-MB-231, and all genes in MCF7, did not display statistically significant differences in gene expression levels. Although not significant, the mean AR, AR-V1, ESR1, and ESR2 relative gene expression levels were higher in 48 h-treated MDA-MB-231 cells (

Transcriptional Profile of Hormone Receptors after Treatment with EENL
Gene expression of ESR1, ESR2, AR, AR-V1, AR-V4, and AR-V7 in breast cells was evaluated before and after 48 h of treatment with 1.0 μg/mL of EENL-EA extract ( Figure 3). AR-V7 expression increased 2.85-fold (P < 0.01) in treated MDA-MB-231 cells after 48 h ( Figure 3H). The remaining genes in MDA-MB-231, and all genes in MCF7, did not display statistically significant differences in gene expression levels. Although not significant, the mean AR, AR-V1, ESR1, and ESR2 relative gene expression levels were higher in 48 h-treated MDA-MB-231 cells (31.74, 2.03, 7.07, and 4.37-fold, respectively).  No significant gene expression changes were found in MCF7 cells treated with EENL-DCM extract at 0.03125 μg/mL for 48 h (Figure 4). However, the TNBC cell line MDA-MB-231 had a 28.41fold increase on the expression of AR-V7 upon treatment (P < 0.05) ( Figure 4F).

In Vivo Experiments with Drosophila Melanogaster
The results showed the carcinogenic action of EENL isolate at the concentrations 0.03125%, 0.0625%, and 0.125%, and the modulating effect of the EENL on the carcinogenic action of doxorubicin (DXR at 0.4 mM). DXR may intercalate on DNA and induce formation of DNA adducts at active promoter sites, increasing torsional stress and enhancing nucleosome turnover. Furthermore, it may trap topoisomerase II at breakage sites, causing double strand breaks. Enhanced nucleosome turnover might increase the exposure of DNA to reactive oxygen species (ROS) resulting in DNA damage and cell death [44]. The frequency of the tumor clone per segment of Drosophila melanogaster is demonstrated in Table 1 and Figure 5. The analyses were performed using the comparative Cq calibrated with data from MCF10 * P < 0.05.

In Vivo Experiments with Drosophila Melanogaster
The results showed the carcinogenic action of EENL isolate at the concentrations 0.03125%, 0.0625%, and 0.125%, and the modulating effect of the EENL on the carcinogenic action of doxorubicin (DXR at 0.4 mM). DXR may intercalate on DNA and induce formation of DNA adducts at active promoter sites, increasing torsional stress and enhancing nucleosome turnover. Furthermore, it may trap topoisomerase II at breakage sites, causing double strand breaks. Enhanced nucleosome turnover might increase the exposure of DNA to reactive oxygen species (ROS) resulting in DNA damage and cell death [44]. The frequency of the tumor clone per segment of Drosophila melanogaster is demonstrated in Table 1 and Figure 5.
The study design included a negative control, which was flies with mutated gene; and DXR as a positive control. For the negative control, the frequency of 0.02 of tumors per fly were observed, and this discrete tumor induction occurs due to the genetic predisposition of the test organism. On the other hand, the positive control induced a frequency of 0.46 of tumors per fly, proving that the organism lineage responded to tumor induction.   The study design included a negative control, which was flies with mutated gene; and DXR as a positive control. For the negative control, the frequency of 0.02 of tumors per fly were observed, and this discrete tumor induction occurs due to the genetic predisposition of the test organism. On the other hand, the positive control induced a frequency of 0.46 of tumors per fly, proving that the organism lineage responded to tumor induction. Figure 5. Frequency of tumors in Drosophila melanogaster. Drosophila melanogaster heterozygous for the warts tumor suppressor gene were treated with different concentrations of isolated Ethanolic Extract from Neem Leaves (EENL) (0.03125%, 0.0625%, and 0.125%) and associated with doxorubicin (DXR) (0.03125%, 0.625%, and 0.125%). The frequency of tumors was analyzed in different segments (eyes, head, wings, body, legs, and halteres). Statistical analysis was according to Mann-Whitney U-test with a significance level of P < 0.05. + value considered statistically different from the negative control (P < 0.05) and ** value considered statistically different from the positive control (DXR 0.4 mM) (P < 0.05).
Larvae that were exposed to the EENL isolate, at concentrations of 0.03125%, 0.0625%, and 0.125% displayed frequencies of 0.21, 0.16, and 0.20 tumors per fly, respectively. Compared to the negative control, a statistically significant increase in tumor induction confirmed the carcinogenic effect of the extract on D. melanogaster. When 0.03125%, 0.0625%, and 0.125% of EENL were applied together with doxorubicin at 0.4 mM, the frequencies of tumors per fly were 5.43, 10.24, and 0.69, respectively. These results demonstrate that EENL with DXR increased tumor frequency at the lowest concentrations, when compared to the positive control. However, in the highest concentration used (0.125%), even when associated with DXR 0.4 mM, the tumor frequency decreased.

Discussion
The association between phenolic compounds in plants and antioxidant activity is due to phenolic hydroxyl groups that have a strong free radical scavenger activity [45][46][47][48][49][50]. Abdelhady et al. (2015) [40], during a phytochemical search for active substances, demonstrated that several phenolic compounds are present in Azadirachta indica.
In fact, natural phytochemicals contain phenolic compounds with the ability to prevent cancer metastasis [36]. Several studies have shown that different flavonoids and polyphenols exert an anti-metastatic effect [37][38][39]. Here we evaluated the cellular effects of ethanolic extracts of neem in breast cell lines. The presence of phenolic compounds in the EENL-DCM and EENL-EA demonstrate their potential antiproliferative activity, which may reduce tumor aggressiveness. The different cellular responses, evaluated by the MTT reduction, after treatments with EENL extracts obtained using DCM and EA, indicated this behavior. For the assays with Drosophila melanogaster, the EENL extracts induced higher tumors at a higher concentration when applied with DXR, probably due to the toxic effect of EENL together with the chemotherapeutic DXR, inducing cell death, which makes the appearance of tumors unfeasible.
According to Paul et al. (2011) [51], the major secondary metabolites present in the Azadirachta indica leaves are nimbolide, vilasinin, nimbinene, 6-deacetyl nimbinene, nimbandiol, nimocinol, β-sitosterol, β-sitosterol-β-D-glicoside, neem leaf glycoprotein, quercetin, glycoside of quercetin, glycoside of kaemferol, quercetin-3-galactoside (hyperin), and rutin. Nimbolide is the most abundant tetranortriterpenoid isolated from leaves of A. indica, showing apoptotic and antiproliferative activity acting in the following pathways: (1) oxidative stress and caspase activation; (2) reduction of the expression of anti-apoptotic proteins (Bcl-xl; Bcl-2) and increasing the expression of pro-apoptotic proteins (Bax, Bad, Bid, cytochrome c); (3) activation of the tumor suppressor p53; (4) activation of the extrinsic apoptosis pathway; (5) inhibition of IGF-1; (6) reducing levels of cyclin-dependent kinase (CDKs) and cyclins, promoting cell cycle arrest; and (7) inhibition of NFκB and its pro-tumorigenesis pathway [52][53][54]. The cytotoxicity profile verified in the MTT assays may also be attributed to the probable induction of apoptosis by the EENL. MCF7 cells were more sensitive to high concentrations of EENL. This behavior was previously observed in prostate cancer, in which hormone-responsive cells (LNCaP) were more sensitive to treatment with this extract [55]. Therefore, analysis of the gene expression, especially of hormone receptor genes, may indicate molecular patterns involved in EENL modulation in breast tumor cells.
The results of Aleskandarany et al. (2016) [56] have shown that there is an association between the expression of AR and good prognosis in BC. Comparing AR expression in HER+, TNBC, and luminal tumors, they observed that luminal BC presented a higher receptor expression. However, AR is an oncogene in triple-negative tumors by "replacing" the estrogen receptor and stimulating tumor growth [57]. After treatment with EENL we observed a decrease of the AR transcripts in the MDA-MB-231 line. This effect, therefore, strongly suggests a possible action of EENL in triple-negative tumors and patients treated with this phenolic compound would have a better prognosis considering the behavior of AR signaling in TNBC [58].
After treatment with EENL-EA and EENL-DCM we detected an increase in AR-V7 transcripts in MDA-M-231. The AR-V7 gene is expressed in primary BC and in breast tumor cell lines, which can promote growth and mediate resistance in androgen deprivation therapies (ADT) in BC subsets [59][60][61]. In this context, the expression of AR and its variants emerge as a new strategy for BC treatment.
Taken together, our results demonstrate a differential effect of EENL in breast tumor cell lines (MCF7 and MDA-MB-231), specially modulating nuclear receptor expression. The behavior of AR-V7 in the MDA-MB-231 tumor cell line indicates new pathways involved in tumor biology, especially as a therapeutic target. Further studies are needed to better understand the role of these compounds in the modulation of such receptors.

Plant Material
The Azadirachta indica (neem) leaves were collected at latitude 18 • 34 13.99 S and longitude 46 • 29 52.57 . The present study was registered in the Genetic Patrimony Management Board (CGEN) (number A2FFD84) and a voucher specimen was recorded in the Herbarium of Institute of Biology of the Federal University of Uberlandia with the number 71869. After identification, the neem leaves were sanitized and oven-dried (Lucaderma), for 72 h at 40 • C. The dry leaves were further milled in a knife mill (Solab).

Ethanolic Extraction and Determination of Phenolic Compounds
Flavonoids and phenolic compounds were extracted according to Cechinel-Filho and Yunes method (1998) [68]. Dried and milled leaves (100 g) were stirred for homogenization in a solution of methanol 50% and protected from light for 20 days.
Extractions were further performed with triple addition of 30 mL of ethyl acetate (EA) P.A. (Proquímicos), and dichloromethane (DCM) P.A. (Alphatec, Carlsbad, CA, USA), separately. After filtration, the filtrates were evaporated in rotary evaporator at 55 • C with reduced pressure and resuspended with DMSO (Sigma Aldrich, St. Louis, MO, USA). All experiments were conducted in six replicates. For in vivo assays, EENL were obtained according to Silva (2010) [69].
Total phenolic compounds (mg gallic acid equivalents (GAE)/g DW) were determined by the Folin-Ciocalteu method [70]. Absorbance analysis was performed in spectrophotometer UV-VIS at the wavelength of 760 nm and dosage was adopted from Swain and Hillis (1959) [71]. Medium (Sigma Aldrich), respectively. MCF 10A was supplemented with 10% of fetal bovine serum (FBS), 20 ng/mL of epidermal growth factor (EGF), 500 ng/mL of hydrocortisone, 10 µg/mL of insulin, and 50 µg/L of gentamycin. The malignant lineages were supplemented with 10% of FBS and 50 µg/L of gentamycin. All lineages were incubated at 37 • C in a humidified atmosphere containing 5% CO 2 atmosphere until they reached 80% of confluence.
After incubation, 10 µL of MTT (5 mg/mL in PBS) (Sigma Chemical Co., St. Louis, MO, USA) was added in each well and incubated for 4 h. Formazan crystals were dissolved in a Sodium Dodecyl Sulfate (SDS) solution containing N-dimethylformamide at 50%. The absorbance was measured at the wavelength of 570 nm in a Thermo Plate Reader (TP-Reader), and the percentage of viable cells was calculated by using the formula in Equation (1). Total RNA was extracted before and after treatments with EENL-EA (1.0 µg/mL) and EENL-DCM (0.03125 µg/mL) for 48 h. Trizol Reagent (Thermo Fisher Scientific, Waltham, MA, USA) was used for extraction according to the manufacturer's recommendations. After isolation, RNA quality was electrophoretic confirmed (1.5% agarose gel) and the concentration was determined using the NanoDrop System (Thermo Fisher Scientific). Reverse transcription was performed as previously described [73].
Quantitative qPCR was performed to evaluate gene expression profiles of ESR1, ESR2, AR, AR-V1, AR-V4, and AR-V7. For normalization, Beta-2-microglobulin (B2M) gene was used and Cq method was validated through relative standard curve construction.

In Vivo Assays
Drosophila melanogaster is a well-established insect model for human diseases and toxicological research due to its well-documented genetics, developmental biology, and capacity to activate enzymatically promutagens and procarcinogens in vivo [75].
The Drosophila tumor suppressor gene warts (wts) encodes a homolog of human myotonic dystrophy kinase and is required for the control of cell shape and proliferation. Homozygous loss of the wts gene of Drosophila, caused by mitotic recombination in somatic cells, leads to the formation of tumors, with aspects of warts appearing on the body of adult flies [76].

Test for the Detection of Epithelial Tumor Clones (Warts)
The test warts (wts) was performed using two strains of Drosophila melanogaster (Diptera: Drosophilidae): (1) Warts lineage (wts/TM3, Sb 1 ) presents a lethal allele wts, on chromosome three, balanced by a chromosome TM3, which is characterized by multiple inversions and marked by dominant mutated stubble (Sb). The stubble mutation is phenotypically characterized by short and thicker hair identified all over the body of the fly. This lineage was kindly provided by the Bloomington Drosophila Stock Center, from Indiana University-USA, registered: Bloomington/7052; (2) Lineage flies multiple wing hairs (mwh) has the marker gene on the chromosome 3 (3-0.3) in a distal position, which is characterized by expressing three or more hairs in each cell. This lineage was kindly provided by Dr. Ulrich Graf (Physiology and Animal Husbandry, Institute of Animal Science, ETH Zurich, Schwerzenbach, Switzerland) [75].
Stocks are grown at the Laboratory of Cytogenetics and Mutagenesis (University Center of Patos de Minas-UNIPAM, Brazil), kept in jars with culture medium for D. melanogaster. The medium is composed of 820 mL of water, 25 g of yeast powder (Sacchoromyces cerevisae), 11 g of agar, 156 g of banana, and 1 g of nipagin. The lineages are maintained inside an incubator B.O.D. 411 D, under light/dark cycles (12 h/12 h), at 25 ± 1 • C and approximately 60% humidity.

Experimental Procedure
To obtain wts +/+ mwh heterozygotic larvae, virgin females wts/TM3, Sb 1 [82] are crossed with mwh/mwh males [83]. The eggs of the descendants were collected for 8 h in flasks containing a culture medium suitable for posture (4% w/v agar-agar based, topped with a thick layer of live baker's yeast supplemented with sucrose). After 72 ± 4 h, third-instar larvae were washed in reverse osmosis water and collected using a fine mesh sieve. The larvae were transferred to a 25 mL vial containing 1.5 g of instant mashed potatoes (HIARI ® brand, São Paulo, Brazil) culture medium. Five milliliters of EENL (in the concentrations 0.03125%, 0.0625%, and 0.125%), diluted in 2.38% ethanol, were added to each tube on its own as well as together with DXR. A co-treatment system was used to apply DXR and EENL. Ethanol 2.38% was used as negative control, and DXR 0.4 mM as the positive control ( Table 3). The treatments were performed in quadruplicate.
These larvae were submitted to a chronic treatment for 48 h. After 120 h of the larval feeding, following metamorphosis, the flies were collected and preserved in 70% ethanol.

Analysis of the Flies
Adult flies of the wts +/+ mhw genotype with long hairs were analyzed for the presence of tumor (wart). The adults with short and thick hairs were discarded, since they did not have the gene being studied. The flies were observed using a stereoscopic magnifying glass and entomological tweezers. Only tumors large enough to be unequivocally classified were recorded. The tumor frequency was calculated as the number of tumors divided by number of wts +/+ mwh flies [75]. A total of 1600 flies were analyzed, 200 flies for each control (positive or negative), 200 flies for each isolated concentration, and 200 flies for each associated group.

Statistical Analysis
Statistical analysis was carried out using GraphPad Prism 7.0 (GraphPad Software, La Jolla, CA, USA). P < 0.05 was considered significant. Paired Student t test was used to compare gene expression and the unpaired for phenols dosage. MTT assays were tested using ANOVA. Finally, the carcinogenic potential of EENL in the D. melanogaster model was validated by the Mann-Whitney and Wilcoxon nonparametric U test, using α = 0.05 level of significance.  Acknowledgments: This work was supported by the following grants: Fundação de Amparo a Pesquisa de Minas Gerais (FAPEMIG) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).

Conflicts of Interest:
The authors declare no conflict of interest.


extfieldsofstudy: Biology Medicine

text:
A demonstration of the H3 trimethylation ChIP-seq analysis of galline follicular mesenchymal cells and male germ cells

Objective Trimethylation of histone 3 (H3) at 4th lysine N-termini (H3K4me3) in gene promoter region was the universal marker of active genes specific to cell lineage. On the contrary, coexistence of trimethylation at 27th lysine (H3K27me3) in the same loci—the bivalent H3K4m3/H3K27me3 was known to suspend the gene transcription in germ cells, and could also be inherited to the developed stem cell. In galline species, throughout example of H3K4m3 and H3K27me3 ChIP-seq analysis was still not provided. We therefore designed and demonstrated such procedures using ChIP-seq and mRNA-seq data of chicken follicular mesenchymal cells and male germ cells. Methods Analytical workflow was designed and provided in this study. ChIP-seq and RNA-seq datasets of follicular mesenchymal cells and male germ cells were acquired and properly preprocessed. Peak calling by Model-based analysis of ChIP-seq 2 was performed to identify H3K4m3 or H3K27me3 enriched regions (Fold-change≥2, FDR≤0.01) in gene promoter regions. Integrative genomics viewer was utilized for cellular retinoic acid binding protein 1 (CRABP1), growth differentiation factor 10 (GDF10), and gremlin 1 (GREM1) gene explorations. Results The acquired results indicated that follicular mesenchymal cells and germ cells shared several unique gene promoter regions enriched with H3K4me3 (5,704 peaks) and also unique regions of bivalent H3K4m3/H3K27me3 shared between all cell types and germ cells (1,909 peaks). Subsequent observation of follicular mesenchyme-specific genes—CRABP1, GDF10, and GREM1 correctly revealed vigorous transcriptions of these genes in follicular mesenchymal cells. As expected, bivalent H3K4m3/H3K27me3 pattern was manifested in gene promoter regions of germ cells, and thus suspended their transcriptions. Conclusion According the results, an example of chicken H3K4m3/H3K27me3 ChIP-seq data analysis was successfully demonstrated in this study. Hopefully, the provided methodology should hereby be useful for galline ChIP-seq data analysis in the future.

INTRODUCTION
Functional cells of multi-cellular organisms maintain their mutual state by adopting specific gene expression patterns. They hereby fix their transcriptional activities by imprinting a particular chromatin structure into genome to allocate and control target gene regions-also recognized as the chromatin packaging [1]. In eukaryotic cell, a chromatin consists of several units of 147 base pairs DNA twining around nucleosome formed by histone octamer-which is the molecule responsible for the chromatin packaging process. A histone octamer comprised of four basic core histones-H2A, H2B, H3, and H4 which are targets for histone modification-methylation, phosphorylation and acetylation. When a precursor cell is committed a lineage-specific development, it modifies particular N-termini amino acids of specific core histones. This event thus results in chromatin packaging to control gene a transcription pattern specific to its lineage-committed cell [1,2]. Different histone modification can relax or condense chromatin structures differently. Since each specific chromatin structure can affect the accessibility of gene transcription factors and RNA polymerase II differently at a gene promoter region, they will result in a gene transcription pattern specific to the lineage-committed cell [1][2][3]. Among well-recognized histone modifications, trimethyl modification of histone 3 (H3) at 4th lysine N-termini (H3K4me3) is very common in active gene transcriptions of several lineage-committed cells [4,5]. This phenomenon, however cannot conclusively comply with several genes in stem cells and germ cells, of which H3K4me3 regions are presented without transcription. Interestingly, such genes always demonstrate not only trimethylation of H3 on 4th lysine (H3K4me3), but also the 27th lysine (H3K27me3) in the same loci of N-termini-the bivalent H3K4m3/ H3K27me3 chromatin or poised chromatin [2,6,7].
In the embryonic stem cells, the presence of H3K4me3 in the poised gene helps marking the gene as the active one. However, the concurrent H3K27me3 repressed the transcription and thus guarantees the cells' pluripotent state [6,7]. Actually, such characteristics were proved to be inherited from germ cells even prior to fertilization. According to this reason, the study of H3K4me3/H3K27me3 bivalent (poised) epigenetic state in male germ cells became a very popular topic in epigenetic evolution studies [2,8].
Chromatin immunoprecipitation (ChIP) assays with sequencing or ChIP-Seq is a sequencing technology generally utilized to study genome-wide H3 trimethylation-H3K4me3 and H3K27me3. In ChIP-Seq, the chromatin regions presented with modified histones of interest are immunoprecipitated and sequenced to study epigenetic control of gene expressions in several cell types [2,5,9,10]. Not only mammals, ChIP-Seq of H3 trimethylation were also conducted in chicken (Gallus gallus) as the representative of galline species [2,[11][12][13]. These studies provided biologists and veterinarians the valuable insights into evolution and cell development of galline species. Of note, an increase of H3K4me3 and H3K27me3 ChIP-seq databases in chicken also allowed both biologists and veterinarians a new opportunity to observe and meta-analyze ChIP-seq data of various cell types acquired from different experiments.
Analytical guidelines are generally concerned in human and mouse ChIP-Seq studies [10]. On the contrary, the issue was limitedly described in H3 trimethylation ChIP-Seq data analysis in galline species [2,12,13]. Since the procedures could affect the analytical results differently, a practical demonstration of such approach should prove beneficial, especially for those who are not acquainted with ChIP-Seq data analysis. Due to such necessity, we therefore demonstrated a chicken ChIP-seq meta-analysis in this study. The H3K4me3 ChIP-Seq data of chicken follicular mesenchymal cells and male germ cells were acquired, and we also included H3K27me3 ChIP-Seq data of chicken germ cells to demonstrate the effect of bivalent H3K4me3/H3K27me3 on gene transcription.

Sample datasets
As far as we knew, the follicular mesenchymal cell was the only lineage-committed cell type with H3 trimethylation ChIP-seq data available in public database. A list of follicular mesenchymal cell and germ cell datasets used in this study is provided (Table 1). Follicular mesenchymal cells used in this study consisted of two populations isolated from median (Med) and lateral mesenchyme (Lat) of feather follicles, accordingly. Male germ cells used in this study derived from two time points during development: pachytene spermatocyte (PS) acquired during prophase of meiosis I, and round spermatid (RS) acquired after completing meiosis. Both germ cell types were representatives of male germ cells to illustrate the poised genes-where their poised states were already known to be retained throughout spermatogenic development [2]. For ChIP-seq data, the sample file and its corresponding control file were provided together. Of note, the control file was not required for RNA-seq data analyses, and thus not provided.

Analytical workflow
In this study, we categorized our analytical processes into 3 steps-data preprocessing, peak calling with visualization, and integrative genomics viewer (IGV) exploration ( Figure  1). In brief, the sample datasets would be aligned to chicken genome in order to acquire ChIP files and RNA-seq files in the data preprocessing step ( Figure 1A). The ChIP files were then applied in peak calling step to identify genome regions enriched with H3 trimethylation-the significant peaks (Figure 1B). Only peaks overlapped with gene promoter region would be considered in IGV exploration of which, the enriched peaks along with the mapped sequences ChIP and RNA-seq files were integrated and visualized together in IGV platform ( Figure 1C). Data preprocessing: The sequence read archive (SRA) files of germ cells and follicular mensenchymal ChIP-Seq and RNA-Seq datasets were retrieved from SRA database (https://www. ncbi.nlm.nih.gov/sra) (Table 1), subsequently extracted, and preprocessed as previously described [14]. The ChIP files were aligned to chicken genome (Gallus gallus 5.0) by Bowtie2 aligner [15], while STAR aligner was applied with the RNA-seq files [14,16]. The aligned files were converted into "bam" format, indexed by gene annotation, removed for duplicated sequences [17], and bias corrected for guanine-cytosine base balance [18] to acquire ChIP files and RNA-seq files. The qualities of both raw and preprocessed datasets were then determined by FastQC as previously described [14]. For ChIP quality determination, normalized strand cross-correlation coefficient (NSC), relative strand cross-correlation coefficient (RSC) and DNA fragment length were additionally included. In more details, NSC is genome-wide correlation between positive and negative strand read counts when shifted by half of the fragment length relative to background. NSC value thus represents enrichment of clustered ChIP fragments around target sites. RSC is the relative enrichment of fragment-length cross-correla- tion to read-length cross-correlation. The RSC thus represented the clustering of relatively fixed sized fragment around target sites.
Peak calling with visualization: Peak calling was, in case of this study, the computational method to identify genomewide chromatin regions with targeted H3 trimethylation significant to the other regions. The peaks therefore represented the genome area with significant presences of H3K4me3 or H3K27me3 when compared to the same region presented in the corresponding control datasets (Table 1). Model-based analysis of ChIP-seq 2 (MACS2) algorithm [19] was used in this study. Since the broad peak regions were suggested in the previous quality analysis, we applied broad peak calling func-tion in MACS2 with all ChIP files. Further details for such function calling was already described [19]. Significant peaks were those with false discovery rate (FDR) ≤0.01 and Foldchange ≥2 presented with chicken gene promoter regions (2,000 base pairs upstream to 200 base pairs downstream of gene promoter) (Peak files in Figure 1B). Peak files were also utilized for counting aggregated DNA within peak region, and counting the DNA fragments distributed from the transcription start site of chicken genes ( Figure 2B).
Integrative genomics viewer (IGV) exploration: The peak files along with their corresponding ChIP files and RNA-seq files were integratively observed together in IGV platform ( Figure 1C). The Gal_gal5.0 genome (UCSC version galGal5 Figure 2. H3 trimethylation loci. The total number of peaks significantly identified were categorized according to their regional chromosomes (A). Each color bar presented in the legend was described as following:  released in December 2015) was administrated as the reference genome. The gal_Gal5.0 consisted of 34 chromosomes, 1 linkage group and 15,411 unplaced scaffolds. Further details of IGV platform can be acquired online (http://software. broadinstitute.org/software/igv/). In brief, DNA aggregations along the genome of sample ChIP files were expressed by blue heatmap. Similarly, the density of transcript fragments from RNA-seq files acquired from the corresponding cell sample were expressed by red heatmap. Peak regions of each cell type were also included in the figure the aggregated DNA fragments from ChIP files. To demonstrate bivalent H3K4me3/ H3K27me3 chromatin of germ cell samples, cellular retinoic acid binding protein 1 (CRABP1), growth differentiation factor 10 (GDF10), and gremlin 1 (GREM1) genes were chosen for illustration ( Figure 4). We select such genes due to the fact that transcription of these genes was crucial in topological control of chicken feather's generation of follicular mesenchyme [11].

All sample datasets were successfully preprocessed and aligned to chicken genome
The raw datasets of all cell types were retrieved and preprocessed according to the methodology ( Figure 1A). All preprocessed datasets (both ChIP-seq and RNA-seq data in Table 1) acquired quality scores across all bases ≥30 without presence of contaminated adapters. Alignment percentages of ChIP-seq and RNA-seq datasets were ≥85% and ≥70%, accordingly. NSC, RSC, and DNA fragment length of ChIP-seq samples are provided in Table 2.

Follicular mesenchymal cells and germ cells shared several gene promoter regions with H3K4me3 and H3K27me3
As previously described in materials and methods, the peaks were the chromatin loci in gene promoter regions enriched with target H3 trimethylations (H3K4me3 or H3K27me3). When categorized the acquired peaks according to their localized chromosomes, chromosome 1 contained the largest number of peaks comparing to the others (chr1 in Figure 2A).
Most of aggregated DNA fragments of ChIP files were found distributed within 0-1 kilobases from the transcription start sites in all cell types ( Figure 2B). The results indicated that that most H3K4me3 peaks (5,704 peaks) were uniquely shared among follicular mesenchymal cells and germ cells (lat_K4, med_K4, PS_K4, and RS_K4), and the second largest intersected peaks were uniquely shared among all H3K4me3 and H3K27me3 (1,909 peaks) (Figure 2A, 3). CRABP1, GDF10, and GREM1 genes of chicken germ cells adopted the bivalent H3 trimethylation chromatin pattern with presence of transcription's inhibition Follicular mesenchymal cells (Lat RNA-seq and Med RNAseq in Figure 4) expressed CRABP1, GDF10, and GREM1 genes higher than germ cells (PS RNA-seq and RS RNA-seq in Figure 4). Of note, the bivalent H3K4me3/H3K27me3 pattern near transcription start site of these gene in germ cells was demonstrated by the significant peak regions (PS_K4 Peak, PS_K27 Peak, RS_K4 Peak, and RS_K27 Peak in Figure 4) with presence of remarkable ChIP-seq DNA aggregations (PS_K4 ChIP-seq, PS_K27 ChIP-seq, RS_K4 ChIP-seq, and RS_K27 ChIP-seq in Figure 4).

DISCUSSION
Data quality was among the most critical issues in ChIP-seq analysis due to its later effect on measurement and comparison. According to this reason, we preprocessed the ChIP-seq data to acquire adequate quality and alignment percentages. Another crucial issue in ChIP-seq quality was the determination of signal-to-noise ratio represented by NSC and RSC values. ChIP-seq data with few genuine binding site regions usually rendered high NSC (≥1.5) and RSC (≥0.8), accordingly. However, marks that tend to be enriched at repeat-like regions with diffused genome-wide patterns like H3 trimethylation usually show lower NSC and RSC values-which also presented in our study [10,11] (Table 2).
Since H3 trimethylation (H3K4me3 or H3K27me3) at gene promoter regions were strongly associated with gene transcription control [4,5,12], only peaks within gene promoter region were hereby considered in this study. Inheritance mechanism of H3K4me3/H3K27me3 bivalent state from germ cells to embryonic stem cells, and activation mechanism of lineage-specific genes in lineage-committed cells were comprehensively described elsewhere [20], and are not included in this discussion. According to the result, most H3 trimethylation peaks in both follicular mesenchymal cells and germ cells were located in chromosome 1-agreeing with the largest gene numbers in the chromosome (2,162 coding genes and 524 non-coding genes) (chr1 in Figure 2). Interestingly, the largest number of intersected H3K4me3 peaks among ChIPseq samples (orange bar in Figure 3) implied the conserved  role of H3K4me3 among chicken germ cells and mesenchymal cells. Of note, shared presented among H3K4me3 and H3K27me3 (red bar in Figure 3) could only partially imply several unique genes controlled by H3K4me3/H3K27me3 , growth differentiation factor 10 (GDF10), and gremlin 1 (GREM1) genes. The heatmaps of CRABP1, GDF10, and GREM1 genes were acquired from the integrative genomics viewer (IGV) interface. The chicken full genome sequences was used to identify the chromatin loci of these genes with transcription direction indicated by the arrow heads presented in the genome along with the associated promoters in the area (galGal5 and Promoter regions at the bottom of the figure). The figure show 3 types of data-ChIP-seq, RNA-seq, and Peak files acquired from each sample cell type. In brief, ChIP-seq data show the density of DNA fragments aligned to a chromatin region (aggregated DNA fragments) in blue heatmap. The RNA-seq data show the density of transcript fragments aligned to the gene regions representing level of gene expressions by red heatmap. Finally, the peak file of each cell sample indicated the chromatin region with significant peaks identifed by peak calling analysis.