Identification of Genes Affected Blue Eggshell Coloration in Xuefeng Black-Bone Chickens

ABSTRACT

Blue eggshells have become economically valuable due to their higher popularity among consumers when compared with eggs of other colors, particularly in East Asia. However, eggshell colors vary widely and gradually become lighter with the increase of age. Therefore, to determine the association between DNA methylation and gene expression on the molecular mechanism of eggshell color variation in Xuefeng Black-bone chickens, we collected the shell glands from chickens that produced dark blue eggshells and light blue eggshells, respectively, and performed whole genome bisulfite sequencing (WGBS) and RNA sequencing (RNA-Seq) on them. The results showed that in the context of CG, DNA methylation levels were the highest compared to CHG and CHH contexts, and were negatively correlated with gene expression levels in three different gene regions in all samples, while there was no significant correlation in the contexts of CHG and CHH. Furthermore, when disregarding the location of coding genes, a total of 55 genes will show differences not only in expression levels but also in DNA methylation levels, and among the majority of differential expression genes (n=50), the DNA methylation levels in the gene body region were the highest. GO analysis and KEGG pathway analysis were performed in these genes and the results showed that the genes NR1H4 and HEPHL1, ABCA4 and ABCA12, and the pathways, cell adhesion molecular, ABC transporters, ECM-receptor interaction, and neuroactive ligand-receptor interaction, may affect the color variation of the blue eggshell of Xuefeng Black-bone chickens.

Keywords:
Blue eggshell; DNA methylation; eggshell glands; transcriptome; Xuefeng Black-bone chicken

INTRODUCTION

The color of poultry eggshells displays enormous diversity, with colors ranging from brown or white to red or blue-green, and the potential for speckling and other patterns (Li et al., 2019). Eggshell colors have multiple functions, such as avoiding predation (Kilner, 2006), strengthening the structure of the shell (Samiullah et al., 2016), filtering harmful solar radiation (Bakken et al., 1978) and decreasing trans-shell microbial contamination (Ishikawa et al., 2010). Despite there being no association between the nutritional content of eggs and the color of their shell, blue eggshell eggs are more popular, since consumers believe eggs with blue eggshell are more organic and nutritious than white and brown eggs laid by commercial hens (Chen et al., 2022).

The main pigments responsible for eggshell coloration are protoporphyrin IX, biliverdin, and Zn-biliverdin chelate (Lang et al., 1987); the difference in their proportions resulting in different eggshell colors. Pink and brown eggshell colors are related to the deposition of protoporphyrin IX, while blue and blue-green eggshell colors are related to the deposition of biliverdin IX and biliverdin IX zinc chelate. Most studies have shown that protoporphyrin IX and biliverdin are synthesized in the eggshell gland, and then secreted and deposited in the eggshell. SLCO1B3 is the key gene regulating the molecular formation mechanism of the blue eggshell trait in chickens (Wragg et al., 2013). More specifically, EAV-HP, an endogenous retrovirus, inserts in the 5’ flanking regions of SLCO1B3, which promotes the expression of theSLCO1B3gene (GenBank Accession No. NM_001318449.1) in the shell gland. This expression causes bile salts to enter the shell gland, from where it is secreted into the eggshell, forming a blue eggshell (Wang et al., 2013).

There are ranges of chicken breeds known to lay blue eggs, such as the native Chinese Dongxiang, Lushi, and Xuefeng Black-bone chickens. Xuefeng Black-bone chickens are the most famous native breed in Hunan province, P.R. China. It is classified as a meat-egg chicken and can lay blue-shelled eggs. However, the eggshell pigment of the blue-shelled eggs laid by Xuefeng Black-bone chickens varies widely, and it tends to become more and more light as age increases, which is typical of epigenetic characteristics. The phenomenon observed in Lushi chickens showed that the methylation of the promoter region of SLCO1B3 was significantly negatively correlated with bothSLCO1B3 expression in the shell gland tissue and eggshell coloration (Li et al., 2019), which substantiated an epigenetic regulation on the eggshell pigment of blue-shelled eggs. However, the molecular mechanism underlying the eggshell coloration in Xuefeng Black-bone chickens remains unknown.

In the present study, RNA sequencing (RNA-seq) and whole-genome bisulfite sequencing (WGBS) jointly used to reveal the association between DNA methylation and gene expression on the molecular mechanism of the eggshell color variation in Xuefeng Black-bone chickens.

MATERIALS AND METHODS

Ethics Statement

All animal experiments followed the Chinese guidelines for animal welfare and the guidelines of the College of Animal Science and Technology, Hunan Agricultural University (Changsha, China). The protocol number was HUNAU2020084.

Animals and Tissue Collection

All Xuefeng Black-bone chickens were picked from the blue eggshell line of Hunan Yunfeifeng Agriculture Co., Ltd (Huaihua, China). Chickens were fed with the same diet, in the same environment. At the age of 64 weeks, 20 chickens that produced dark blue eggshells (group S) and another 20 chickens that produced light blue eggshells (group Q) were selected. Then their eggs were collected for a week to measure eggshell color. Shell gland samples were collected from chickens with a significant difference in eggshell color and relatively uniform laying times. The samples used for further RNA and DNA extraction were cleaned quickly on ice and immediately put into a liquid nitrogen tank. They were then transferred to a refrigerator at -80ºC until they were used.

RNA-seq

RNA extraction, Library Construction, and Transcriptome Sequencing

According to the protocol, the total RNA from 4 shell gland samples was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA). RNA quality was estimated on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase-free agarose gel electrophoresis. After the total RNA was extracted, eukaryotic mRNA was enriched with Oligo(dT) beads. The enriched mRNA was then fragmented into short fragments using fragmentation buffer, and reversely transcribed into cDNA by using NEBNext Ultra RNA Library Prep Kit for Illumina (NEB #7530, New England Biolabs, Ipswich, MA, USA). The purified double-stranded cDNA fragments were end-repaired, and a base was added and ligated to Illumina sequencing adapters. The ligation reaction was purified with AMPure XP Beads (1.0X). Ligated fragments were subjected to size selection by agarose gel electrophoresis and polymerase chain reaction (PCR) amplified. The resulting cDNA library was sequenced using Illumina Novaseq6000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Sequence Read Quality Control, Align, and Quantification of Gene Abundance

Read were further filtered by fastp to obtain high-quality, clean reads (Chen et al., 2018) (version 0.18.0). Raw reads that contained adapter, more than 10% unknown nucleotides (N), or low-quality reads were removed. The short reads alignment tool Bowtie2 (Langmead et al., 2012) (version 2.2.8) was subsequently used for mapping reads to the ribosome RNA (rRNA) database, and the rRNA mapped reads were removed. The remaining clean reads were aligned to the reference genome (Anas platyrhynchos) by using HISAT2. 2.4 (Kim et al., 2015) with “-RNA-strandedness RF”. The mapped reads of each sample were then assembled by using StringTie v1.3.1 (Pertea et al., 2015) in a reference-based approach. For each transcription region, an FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated to quantify its expression abundance and variations, using the RSEM (Li et al., 2011) software. The FPKM formula was as follows:

$FPKM=\frac{{10}^{6}C}{NL/{10}^3}$

Where C is the number of fragments mapped to gene A, N is the total number of fragments mapped to reference genes, and L is the number of bases on gene A.

Principal Component Analysis (PCA)

Principal component analysis (PCA) was performed with R package gmodels2.18.1 (http://www.r-project.org/). We input the data to the R package, and the software gives the PCA result through the following steps:

Step 1. Standardization of the range of continuous initial variables.

Step 2. Computation of the covariance matrix to identify correlations.

Step 3. Calculation of the eigenvectors and eigenvalues of the covariance matrix to identify the principal components.

Step 4. Creation of a feature vector to determine which principal components should be kept.

Step 5. Recasting of the data along the principal components’ axes.

Differentially Expressed Genes (DEGs)

RNAs differential expression analysis was used through the DESeq2 (Love et al., 2014) software between two different groups (and through edgeR (Robinson et al., 2010) between two samples). Genes/transcripts with the parameter of false discovery rate (FDR) below 0.05 and absolute fold change ≥2 were considered differentially expressed genes/transcripts.

Gene Ontology Term Enrichment and Kyoto Encyclopedia of Genes and Genomes Analysis

All DEGs were also analyzed through GO and KEGG tools to recognize the main biological functions they exercise. The terms and pathways with a P-value not greater than 0.05 (p≤0.05) were defined as significantly enriched GO terms and pathways in DEGs.

cDNA Synthesis and Quantitative Real-Time PCR (qRT-PCR)

RNA was extracted following the same procedure described above, and was subsequently converted into cDNA by using the TAKARA PrimeScriptTM RT reagent kit (TaKaRa). qRT-PCR primers were designed by NCBI (https://www.ncbi.nlm.nih.gov/) (Supplementary Table S1). b-actin was used as an internal control to normalize the genes’ expression level and the 2^-△△Ct method was used to analyze each gene’s relative expression levels.

WGBS

DNA Extraction, Library Construction, and Whole Genome Bisulfite Sequencing (WGBS)

Genomic DNA was extracted using a DNA kit (Tiangen, China) according to the manufacturer’s protocols. DNA concentration was detected by NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and integrity was tested by agarose gel electrophoresis. Then, the genomic DNAs were fragmented into 100-300bp by Sonication (Covaris, Massachusetts, USA) and purified with a MiniElute PCR purification kit (QIAGEN, MD, USA). The fragmented DNAs were then end-repaired and a single “A” nucleotide was added to the 3’ end of the blunt fragments. The genomic fragments were subsequently ligated to methylated sequencing adapters. Fragments with adapters were bisulfite converted using the methylation-gold kit (ZYMO, CA, USA), and unmethylated cytosine was converted to uracil during sodium bisulfite treatment. Finally, the converted DNA fragments were PCR amplified and sequenced using Illumina HiSeqTM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Methylation Level Analysis

High quality clean reads were obtained after filtering raw reads. The reads that contained more than 10% unknown nucleotides (N) and 40% low quality (Q-value ≤ 20) were removed. The clean reads obtained from the raw reads were mapped to the species reference genome using BSMAP software (Xi et al., 2009) (version: 2.90) by default. Methylated cytosines were called using a custom Perl script and tested with the correction algorithm described by Lister R. et al (2009). The methylation level was determined based on the percentage of methylated cytosine in the overall genome, per chromosome, and in the different regions of the genome in different sequence contexts (CG, CHG and CHH). To estimate different methylation patterns in different genomic regions, the methylation profile at flanking 2 kb regions and gene body (or transposable elements) was drawn based on the average methylation levels per window.

Differentially Methylated Cytosines (DMCs) Analysis

Pearson’s chi-square test (χ2) in methylKit (Akalin et al., 2012) (version: 1.7.10) was carried out to determine differential DNA methylation between the four samples at each locus. To identify DMCs, the minimum read coverage to call a methylation status for a base was set to 4. Differentially methylated cytosines for each sequence context (CG, CHG and CHH) between two samples were identified according to different criteria: 1) For CG and CHG, absolute value of the difference in methylation ratio| ≥ 0.25, and q ≤ 0.05; 2) For CHH, absolute value of the difference in methylation ratio| ≥ 0.15, and q ≤ 0.05; 3) For all C, absolute value of the difference in methylation ratio| ≥ 0.2, and q ≤ 0.05.

Differentially Methylated Regions (DMRs) Analysis

To identify DMRs between four samples, the minimum read coverage was set to 4, which can be called a methylation status for a base. DMRs for different sequence contexts (CG, CHG and CHH) were determined according to different criteria: 1) For CG and CHG, numbers of GC and CHG per window ≥ 5, absolute value of the difference in methylation ratio ≥ 0.25, and q ≤ 0.05; 2) For CHH, numbers in a window≥ 15, value ratio ≥ 0.15, and q ≤ 0.05; 3) For all C, numbers in a window ≥ 20, value ratio ≥ 0.2, and q ≤ 0.05.

Association analysis of RNA-seq and WGBS

Correlation of DNA Methylation and Gene Expression in Samples

To determine whether gene expression influences DNA methylation in a sample, genes were categorized into four classes based on their expression levels, including a non-expression group (RPKM (reads per kilobase per million reads mapped) ≤ 1), a low-expression group (1 < RPKM ≤ 10), a middle-expression group (10 < RPKM ≤ 100), and a high-expressin group (RPKM > 100).

To analyze whether DNA methylation influences gene expression in a sample, genes were categorized into four classes according to their methylation level, including a non-methylation group, a low-methylation group, a middle-methylation group, and a high-methylation group (non-methylation was not considered, and genes were divided into 3 groups on average).

The relationships between DNA methylation and gene expression within the ±2kb flanking regions and gene body region were statistically discerned using Spearman correlation analysis. rho > 0 means a positive correlation, while rho < 0 means a negative correlation.

Correlation of DNA Methylation and Gene Expression Between Groups

To analyze whether differentially expressed genes (DEGs) influence DNA methylation between groups, DEGs were categorized into four classes based on their differential expression pattern, including a special-down group (genes specifically expressed in the control group), and a special-up group (genes specifically expressed in the treatment group), an other-down group (genes down-regulated expressed in the treatment group), an other-up group (genes up-regulated expressed in the treatment group).

To determine whether DNA methylation levels in differentially methylated regions (DMRs) influence gene expression between groups, genes were classified according to genomic location. including the ±2kb flanking regions and the gene body region.

Correlation of DMRs and DEGs Between Groups

To explore the potential functions of DNA methylation responsible for DEGs, common genes between DMRs related genes and DEGs were analyzed. Gene ontology (GO) enrichment analysis and KEGG pathway enrichment analysis were also conducted for DEGs with DMRs.

RESULTS

The Identification of Chickens that Produced Dark-blue and Light-Blue Eggshells

To distinguish the Xuefeng Black-bone chickens that produced dark-blue and light-blue eggshells, we used a colorimeter to measure the color of eggshells (obtuse, medium, and acute ends) of 20 chickens from group S and group Q, respectively. The color of the blue eggshell was measured as shown in Table 1. The L*(lightness) and a*(redness) values in the eggshell color of Group S were extremely significantly lower than those of Group Q (p<0.01), which indicates that the selection of hens in the current study was successful.

Table 1
Eggshell color in laying hens with dark blue eggshells or light blue eggshells.

Items	Group		SEM	pp-value
	S	Q
L*	63.1	72.3	1.02	0.00
a*	-3.59	-0.75	0.22	0.00
b*	6.47	7.30	0.49	0.95

• S, dark blue eggshell; Q, light blue eggshell; SEM, standard error of the mean. p values <0.05 were considered significant and p values <0.01 were considered extremely significant.

Overview of RNA-seq Results

To investigate the molecular mechanism of blue eggshell color variations, transcriptome analysis of shell glands were conducted in group Q and group S. Transcriptome sequence data for all samples can be found in NCBI (https://www.ncbi.nlm.nih.gov/sra/PRJNA1022583, SUB13889596). The important results of RNA-seq are shown in Tables 2 and 3 and Fig. 1. According to the tables, the number of clean reads generated from each library were 45374908, 47888306, 51922410 and 57728652, respectively. The Q20 value was more than 97.83% (Table 2). The data showed that the sequencing results were of good quality and reliable, being suitable for use in subsequent bioinformatics analysis. After the low-quality reads, adaptor sequences and rRNA reads were removed from the total reads. 89.26-91.14% of clean reads were mapped uniquely to the reference (Anas platyrhynchos) genome, and only a small proportion of them was mapped to multiple locations in the genome (Table 3), which showed that the degree of gene coverage was high, the annotation was perfect, and the sequencing result was reliable. Principal component analysis (PCA) also showed global differences between group Q and group S (Fig. 1A). All evidence suggested that our data was repeatable and reproducible. A total of 233 DEGs were downregulated, and 205 DEGs were upregulated in the Q-vs-S comparison (Fig. 1B). A Volcano plot of the DEGs showed their separate profiles between the Q and S group, and the results showed that the |log2FC| value of most differential genes ranged from -2 to 2 (Fig. 1C).

Table 2
Summary of the RNA-seq data collected from Q and S.

Sample name	Clean read num	Adapter (%)	Low-quality (%)	Poly-A (%)	N (%)	Q20%	GC%
Q1	51922410	34008(0.07%)	231288(0.44%)	0(0%)	1354(0%)	7582937593(97.83%)	3859227762(49.79%)
Q2	57728652	36318(0.06%)	218304(0.38%)	0(0%)	1546(0%)	8446124937(97.99%)	4231072105(49.09%)
S1	47888306	27240(0.06%)	177948(0.37%)	0(0%)	1310(0%)	7002172783(97.97%)	3517380655(49.22%)
S2	45374908	23070(0.05%)	158830(0.35%)	0(0%)	1300(0%)	6645637031(98.03%)	3335550893(49.20%)

• N, unknown base rates higher than 10%; Q20, sequencing error rates lower than 1%.

Table 3
Summary of clean reads and genes mapped to the reference genome from Q and S.

Sample name	Total reads	Mapping ratio	Multiple mapped	Unique mapped	Known genes	Know genes in four samples	Novel number	Total genes
Q1	51922410	91.11%	960492 (1.85%)	46236417 (89.26%)	14205(85.74%)	15315 (92.44%)	820	15025
Q2	57728652	92.29%	1130698 (1.96%)	52035617 (90.32%)	14301(86.32%)		831	15132
S1	47888306	92.97%	877125 (1.84%)	43546340 (91.14%)	14237(85.93%)		818	15055
S2	45374908	92.79%	896144 (1.98%)	41136415 (90.81%)	14096(85.08%)		808	14904

Figure 1
Overview of RNA-seq results. (A) PCA of different samples based on the normalized expression levels of all expressed genes. (B) Summary of the number of upregulated and downregulated genes in the Q-vs-S comparisons. (C)Volcano plot showing p-values against fold changes.

DEG Functional Annotation

The DEGs from Q and S were annotated using the Gene Ontology (GO) database. Among the biological processes, the cellular process, biological regulation, metabolic process and signaling were the most frequent terms in both comparisons (Fig. 2A). Moreover, the terms related to eggshell pigmentation were also enriched, though they were not the most frequent terms, which were iron ion transport, transition metal ion transport and negative regulation of bile acid biosynthetic process; some relative genes, namely MELTF, CP, HEPH, HEPHL1, ATP2C2, NRIH4, and FGF19 were involved in them (Table S2). Binding, catalytic activity, and molecular function regulator were found to appear most frequently in the molecular function in both comparisons (Fig.2B). Furthermore, the GO term transporter activity, which may be related to biliverdin transport, was enriched, and the related genes were ASIC2, RBP4A, SLC1A6, SLC16A3, SLC37A2, MCOLN2, GC and KCNQ5 and so on (Table S2). In the cellular component, cell and cell part were the top terms in both comparisons (Fig.2C), and the terms secretory granule and secretory vesicle, which were related with secretory, with target genes ABCA12 and SLC2A14, were also identified (Table S2).

Figure 2
Gene Ontology (GO) classification of DEGs (Q-vs-S). (A)Biological Processes. (B)Molecular Functions. (C)Cellular Components. The results are summarized in three main categories: biological processes, cellular components and molecular functions. The x-axis indicates the second-level GO term, and the y-axis indicates the number of genes.

We used the KEGG database to further understand genes’ biological functions. As shown in Fig. 3, porphyrin and chlorophyll metabolism, drug metabolism-cytochrome P450, and metabolism of xenobiotics by cytochrome P450 were significantly enriched in the Q-vs-S comparison (p<0.05). Through transcriptome sequencing of the liver of blue shell laying hens, Chen et al (2021). found that the differential genes were enriched in drug metabolism-cytochrome P450, which suggested that they are related to biliverdin synthesis. In the porphyrin and chlorophyll metabolism pathway, the target genes were UGT1A1, HEPG, ALAD, and CP (Table S3); while in the drug metabolism-cytochrome P450 pathway and metabolism of xenobiotics by cytochrome P450 pathway, the enriched related genes were GSTA3, ADH1C, HPGDS, and UGT1A1 (Table S3).

Figure 3
Top 20 pathways in KEGG enrichment by QValue (Q-vs-S). Gene ratio is the ratio of the differentially expressed number of genes in the pathway and the total number of genes in the pathway. The higher the gene radio, the higher the degree of enrichment. QValue is the P-value after the multiple hypothesis test correction, in the range of 0 to 1; the closer the QValue is to zero, the more significant the enrichment.

Quantitative Real-time Polymerase Chain Reaction Analysis

To ensure the accuracy of differential expression data obtained from the RNA-Seq, 13 genes were selected for qRT-PCR assay using the same RNA samples used for RNA-seq, comprising 6 down-regulated genes (UGT1A1, ALAD, HEPH, SLC35B4, SLC2A14, CA8), 6 up-regulated genes (CP, ARFGAP3, CD164, FADS2, SLC16A3, AGPAT4), and 1 gene with no significant differential expression (SLCO1B3). Pearson’s correlation coefficient was used to evaluate the correlation between the sequencing data and the quantitative data. The results showed that the correlation coefficient is 0.913, which was very high, indicating that RNA-seq analysis generated dependable data (Fig. 4).

Figure 4
qRT-PCR analysis for Up-regulation and down-regulation of DEGs (Q vs S).

The DNA Methylation Atlas of Q and S

The color of the blue eggshell becomes lighter with the increase of the weeks of age, which suggests that DNA methylation might be involved in the regulation of eggshell color. To ascertain the role of DNA methylation in Xuefeng Black-bone chickens’ shell gland, WGBS was used to characterize their cytosine methylation patterns. The raw datasets of WGBS can be accessed at the NCBI (https://www.ncbi.nlm.nih.gov/sra/PRJNA1022583, SUB13889848). In the present study, the number of raw reads generated from each library ranged from 179409670 to 209269624 (Table 4). After removing the low quality, N (unknown), and connector contamination reads, the number of clean reads finally obtained in the Q and S groups was between 177689330 and 206954352 (Table 4). 85.27%, 84.51%, 82.88%, and 84.86% of chicken genomes were covered with the uniquely mapped reads of the Q1, Q2, S1, and S2 groups, respectively (Table 4). The sequencing depth was more than 21.31. As shown in Table 5, the effective coverage of functional element regions in each sample was high in three contexts (CG, CHG and CHH). The Q20 value was more than 0.96. These results showed that the sequencing results were reliable. Furthermore, the Circos plot showed the DNA methylation levels in the different sequence contexts (mCG, mCHG, and mCHH) (where H is A, C or T) in different chicken chromosomes (1-33 and the Z, W, MT chromosome; Fig. 5).

Table 4
The summary of data generated by WGBS.

Sample name	Raw reads	Clean reads	Clean rate (%)	Mapping rate (%)	Sequence depth	Q20
Q1	179409670	177689330	99.04	85.27	21.64	0.96
Q2	198510066	196622472	99.05	84.51	23.73	0.96
S1	183920368	180062822	97.90	82.88	21.31	0.96
S2	209269624	206954352	98.89	84.86	25.09	0.96

Table 5
Effective coverage of functional element regions of each genome.

Sample name	Contexts	Genebody	Upstream_2k	Downstream_2k	Five_UTR	Three_UTR	Exon	CDS	Intron
Q1	CG_rate(%)	78.29	56.16	64.23	49.00	78.23	67.29	69.46	81.08
	CHG_rate(%)	90.98	74.13	81.35	64.74	88.98	82.12	82.25	91.95
	CHH_rate(%)	91.11	74.57	80.93	67.11	88.16	83.35	83.36	91.84
Q2	CG_rate(%)	72.51	45.46	55.87	37.84	71.80	58.49	60.74	76.08
	CHG_rate(%)	87.47	66.11	75.02	55.74	84.37	75.98	76.00	88.75
	CHH_rate(%)	87.90	67.22	75.35	59.02	84.00	77.83	77.63	88.86
S1	CG_rate(%)	76.84	53.88	62.10	47.63	76.43	65.18	67.06	79.80
	CHG_rate(%)	90.01	72.14	79.51	63.23	87.60	80.47	80.48	91.06
	CHH_rate(%)	90.29	72.96	79.45	65.55	86.95	81.90	81.81	91.09
S2	CG_rate(%)	80.85	60.57	68.17	54.12	81.06	71.14	73.16	83.31
	CHG_rate(%)	92.46	77.46	83.83	68.84	90.60	84.61	84.77	93.32
	CHH_rate(%)	92.50	77.62	83.29	70.88	89.71	85.61	85.70	93.14

Figure 5
Distribution of identified methylation sites on each chromosome. The outer ring represents the chicken genome labeled with chromosome number and position. (A-D) CG Methylation;(E-H) CHG Methylation; (I,J,M,S) CHH Methylation;(A,B,E,F,I,J) Q;(C,D,G,H,K,L) S.

Global DNA Methylation Patterns

In the present study, we analyzed the DNA methylation levels in three contexts: CG, CHG, and CHH (where H is A, C, or T) to analyze the differences in global DNA methylation profiles between the two study groups. As shown in Table 6, the highest proportion (58.61%-66.18%) of cytosines were methylated in the CG context, and only a small proportion (1.45%-1.51%) of cytosines were methylated in the non-CG contexts (CHG and CHH contexts). Meanwhile, we analyzed the proportion of methylation levels in the CG, CHG and CHH contexts to the total methylation levels in the whole genome (Fig. 6A), and we found that the CG context (79.58%-85.29%) accounted for the most of them, as compared to the CHG and CHH contexts (less than 20%). To investigate the patterns of methylated cytosines in the chickens’ shell glands, we analyzed the mC sequence preferences in different sequence contexts. Our results indicated that the methylated cytosines tend to be located in the CG, CHG, and CHH (H = A > T) (Fig. 6B). Methylation levels for specific functional areas (CGI, Promoter, Exon, Intron, Utr5prime, Utr3prime) are shown in Fig. 6C. We found that the methylation levels in the exon, intron, and utr3prime regions were higher than in others, and there was no change in the CHG and CHH contexts.

Figure 6
Global DNA methylation profiles. (A) Comparison of DNA methylation patterns in different samples. (B) Sequence preferences for methylation in various sequence contexts. 9 bp base information around the position of mCG, mCHG, mCHH at the highest or lowest methylation levels, in which the methylated cytosine is in the fourth position. (C) Global DNA methylation levels in different specific functional regions between different samples.

Table 6
Statistical status of c-site methylation in the whole genome (Q vs S).

Content	Q1 (%)	Q2 (%)	S1 (%)	S2 (%)
mC	3.37	3.30	3.06	3.37
mCG	66.18	65.08	58.61	65.20
mCHG	0.70	0.72	0.69	0.70
mCHH	0.76	0.79	0.78	0.75

Analysis of Differentially Methylated Locus (DMCs) and Differentially Methylated Region (DMRs)

In the present study, a total of 5832 DMCs were discovered. The number of DMCs that were upregulated or downregulated in the four contexts (C, CG, CHG, CHH) are shown in Fig. 7. GO analysis and KEGG pathway analysis were performed to look into the DMCs and DMRs related genes’ potential biological roles. Our results showed that the GO analysis of the DMCs-related genes was mainly enriched in the signaling, biological regulation, and regulation of biological process in the biological process. Meanwhile, in the cellular components, cell, cell part and organelle were the most frequent terms. Regarding molecular functions, binding and catalytic activity were the main ones enriched (Fig. 8A, B, C). Moreover, the terms related to biliverdin synthesis and transport, negative regulation of bile acid metabolic process, and ferrous iron transmembrane transporter activity were also enriched, and the related genes were SLC40A1, MALRD1 and PROX1 (Table S2). The terms that were enriched in the three contexts (CG, CHG, CHH) were similar (Fig. 8A, B, C).

Figure 7
Number of DMCs. The number of DMCs that were upregulated or downregulated in four contexts (C, CG, CHG, CHH).

Figure 8
Analysis of DMCs. (A) GO enrichment analysis of DMCs in CG context. (B) GO enrichment analysis of DMCs in CHG context. (C) GO enrichment analysis of DMCs in CHH context. (D) KEGG enrichment analysis of DMCs in CG context. (E) KEGG enrichment analysis of DMCs in CHG context. (F) KEGG enrichment analysis of DMCs in CHH context.

The results of the KEGG pathway analysis showed that 121, 67, and 103 pathways were enriched in the CG, CHG, and CHH contexts, respectively. In CG, adrenergic signaling in cardiomyocytes, AGE-RAGE signaling pathway in diabetic complications, and ubiquitin mediated proteolysis were the most frequent pathways (Fig. 8D); in CHG, the genes were enriched in ECM-receptor interaction, RNA degradation, and insulin signaling pathway (Fig. 8E); in CHH, the genes were enriched in adrenergic signaling in cardiomyocytes, and protein export (Fig. 8F). Furthermore, genes related to eggshell pigmentation were also enriched in relevant terms and pathways, such as ErbB signaling pathway, cell adhesion molecules, ABC transporters, porphyrin and chlorophyll metabolism, neuroactive ligand-receptor interaction, and so on. Meanwhile, porphyrin and chlorophyll metabolisms were also identified in the CG context, and the related genes were EPRS and GLRB (Table S3).

We also identified DMRs to study the differential methylation between differently colored eggshells. We detected 6384 DMRs between group Q and S in four different contexts (C, CG, CHG, CHH), with most DMRs occurring in the CG context (Fig. 9A). Therefore, we carried out the analysis of DMR related genes in the CG context. Among the biological roles, the GO analysis of DMRs related genes showed those involved in cell process, biological regulation, biological regulation process, metabolic process and multicellular organismal process and so on; among cell components, membrane, cell junction, and extracellular region were enriched and may be associated with transmembrane transport of pigments; and among molecular functions, transporter activity, cargo receptor activity, and binding may be related to material transportation (Fig. 9B). Interestingly, glycine-gated chloride ion channel activity targeting genes GLRA2 and GLRA3 was enriched in molecular function (Table S3). The results of KEGG pathway analysis showed that the DMRs were enriched mostly in focal adhesion, adrenergic signaling in cardiomyocytes, and gap junction (Fig. 9C). Among these terms and pathways, ABC transporters were related to eggshell pigment, and genes GLRA2 and GLRA3 were involved with it (Table S3).

Figure 9
Analysis of DMRs. (A) The number of DMRs that were upregulated or downregulated in four contexts (C, CG, CHG, CHH). (B) GO enrichment analysis of DMRs in CG context. (C) KEGG enrichment analysis of DMRs in CG context.

Correlation Analysis Between DNA Methylation and Gene Expression

The data from WGBS and RNA-Seq were jointly analyzed to further the current understanding of the relationships between transcriptome and methylation of molecular mechanisms regulating blue eggshell color differences among Xuefeng Black-bone chickens. Firstly, we calculated the average methylation rates of four types of genes (none, low, middle, and high expression) in the gene body, the upstream 2kb of the transcription start site (TSS), and the downstream 2kb of transcription termination sites (TTS). The results showed that in the context of CG, the methylation level of the genes with high expression levels in the three gene regions (the upstream 2kb of the TSS, gene body, and the downstream 2kb of the TSS) of the four samples was low, and the methylation level of the upstream 2kb of the gene TSS was the lowest. There was no significant difference in DNA methylation levels in the three regions (the upstream 2kb of the TSS, gene body, and the downstream 2kb of the TTS) in both the CHG and CHH contexts (Fig. 10 and Fig. 11). Then we performed Spearman correlation analysis for gene expression levels and DNA methylation levels in the three gene regions. We found that in the context of CG, DNA methylation levels were negatively correlated with gene expression levels in three different gene regions (the upstream 2kb of the TSS, gene body and the downstream 2kb of the TTS) in all samples, while there was no significant correlation between DNA methylation levels and gene expression levels in the contexts of CHG and CHH. Moreover, there was a negative correlation between DNA methylation rates and gene expression levels in the upstream 2kb of the TSS gene in the three contexts (CG, CHG, CHH) (Fig. 12A).

Figure 10
Methylation levels of different expression levels of genes in different regions of Q group. (A) Q1 mCG. (B)Q1 mCHG. (C)Q1 mCHH. (D)Q2 mCG. (E)Q2 mCHG. (F)Q2 mCHH.

Figure 11
Methylation levels of different expression levels of genes in different regions of S group.(A)S1 mCG. (B)S1 mCHG. (C)S1 mCHH. (D)S2 mCG. (E)S2 mCHG. (F)S2 mCHH.

Analysis of Genes Shared by DEGs and DMR-Related Genes

According to the analysis of gene expression level and methylation rate, DNA methylation regulated gene expression in the CG environment. Therefore, an interaction analysis of common genes between DMR and DEG was carried out in the context of CG. Disregarding the location of the coding genes, a total of 55 genes showed differences, not only in expression levels, but also in methylation levels. Among them, the majority (50 genes) was in the gene body region. There were 3 common genes in the upstream region and 6 common genes in the downstream region (Fig.12B), indicating that most of the differential expression regulated by DNA methylation in the eggshell gland of Xuefeng Black-bone chickens occurred in the gene body region. To look into common genes’ potential biological roles, GO analysis and KEGG pathway analysis were performed. GO enrichment analysis results showed that, among the biological processes, cellular process, biological regulation and regulation of biological process were the most enriched terms. Signaling and intracellular bile acid receptor signaling pathway were also enriched in the biological processes, which were related to biliverdin secretion (Table S2). Membrane, membrane part, and cell junction were the terms enriched in among the cellular components. Finally, in terms of the molecular functions, transporter activity, transcription regulator activity, bile acid receptor activity, and ferroxidase activity were enriched, which were related to biliverdin synthesis and porphyrin metabolism (Fig. 13A). The related genes involved in them were NR1H4 and HEPHL1 (Table S2). Additionally, KEGG enrichment analysis results showed that cell molecular adhesion, ABC transporters, ECM-receptor interaction and neuroactive ligand-receptor interaction were the enriched pathways (Fig. 13B), and the ABC transporters pathway was related to the synthesis of biliverdin. The related genes were ABCA4 and ABCA12 (Table S3).

Figure 12
Conjoint analysis of DNA methylation and transcriptome. (A) Correlation between DNA methylation and gene expression. (B) Venn diagram of DMR and DEG associated genes in CG sequence.

Figure 13
Analysis of common genes of DMRs and DEGs. (A) GO enrichment analysis. (B) KEGG enrichment analysis.

DISCUSSION

Consumers currently believe that blue eggshell eggs are produced by local chicken breeds, the feeding method of which is often organic, harmless and cage-free. Therefore, blue eggshell eggs are more popular than eggs of other eggshell colors in China, and there is a positive association between the degree of blue color and its popularity. For consumer purposes, darker blue eggshells mean higher quality eggs. This has forced sellers to select darker blue eggs. However, the eggshell color of eggs laid by the same chicken varies greatly at different weeks of age, which is typical of epigenetic features. Nevertheless, the molecular mechanism of blue eggshell color variations is largely unknown.

In this study, we detected global gene expression and DNA methylation profiles in the shell gland of dark and light blue eggshells, providing large amounts of information for further studies on the regulatory mechanisms underlying poultry eggshell color changes. According to the RNA-seq data, our results revealed that many genes related to biliverdin biosynthesis and transport, such as UDP glucuronosyltransferase family 1 member A1 (UGT1A1), ceruloplasmin (CP), hephaestin (HEPH), aminolevulinate dehydratase (ALAD), nuclear receptor subfamily 1 group H member 4 (NR1H4), albumin (ALB), and so on. Previous studies have made preliminary investigations on the functions of these genes. Among them, ALAD is the rate-limiting enzyme in the process of biliverdin synthesis (Pranawidjaja et al., 2015). Under the catalysis of ALAD, 5-aminolevulinic acid is converted into porphobilinogen, which is the raw material for biliverdin synthesis; CP and HEPH play an important role in the process of the redox of Fe²⁺ (Johannesson et al., 2012; Wierzbicka et al., 2014), which is a necessary step in the oxidation process of protoporphyrin IX to heme; UGT1A1 is responsible for metabolizing bilirubin into bilirubin b-diglucuronide (Qi et al., 2020; Zhang et al., 2021); and NR1H4 and ALB are involved in bile acid and bile salt recycling process. Researches in human medicine have shown that the mutation of NR1H4 can cause obstruction of bile salt transportation, leading to bile stasis (Davit-Spraul et al., 2012; Czubkowski et al., 2021). Therefore, it is speculated that the expression of this gene will promote the transport of bile salts and bile acids, increase the content of biliverdin synthesis raw materials, and thus regulate the formation of the eggshell color. ALB, as a carrier of bile salts, may play a role in the transport of raw material, and further affect the eggshell color. These genes can be considered potential candidate genes for regulating eggshell pigment deposition.

In the present study, functional annotation analysis of the DEGs revealed that these genes play important roles in some biliverdin metabolism- and transport-related pathways, including the porphyrin metabolism pathway, scavenging of heme from plasma pathway, bile acid and salt recycling pathway, cell adhesion molecules (CAMs), and neuroactive ligand-receptor interaction. Porphyrin metabolism and scavenging of heme from the plasma pathway play important roles in the biliverdin synthesis, which is a derivative of heme and results from heme oxygenase-1 (HO-1) activity (Maines, 1997; Elbirt et al., 1999). As a bile pigment, the recycling process of bile acid and salt provides the raw material for biliverdin synthesis. Bai et al. (2019) also screened CAMs and neuroactive ligand-receptor interaction pathways when they screened Putian duck green shell eggs related genes using RNA-Seq technology, which suggested the CAMs may play an important role in the eggshell color. Additionally, many SLC family members and ABC transporter family members, including SLC35F3, SLC1A6 and ABCA12 were enriched. Similar results were discovered in the study of the shell gland of ducks that also produced blue eggs (Zheng et al., 2014), which indicated that the SLC family members and the ABC transporter family may have an important influence on the difference in blue eggshell colors.

DNA methylation is a covalent epigenetic modification of cytosine to 5-methylcytosine, occurring within CpG dinucleotides. DNA methylation usually occurs in CpG regions, and gene expression is usually inhibited when methylation occurs in promoter or enhancer regions. According to our WGBS data, over 79% of mC were found to be located in the CG context, and other sequence contexts (CHG, CHH) accounted for less than 20%. We observed that the methylation level of the exon, and intron regions and 3’UTR of the chicken genome consisted of a large proportion, while the promoter region was relatively low, which is similar to a previous study (Ball et al., 2009). However, the relationship between gene body methylation and gene expression is still not clear. In previous studies, a positive correlation was seen in cancer cell lines (Shann et al., 2008), and a negative correlation was seen in human B cells (Rauch et al., 2009). In the present study, the gene body methylation level had a negative correlation with gene expression in the CG context. The effect of DNA methylation levels in gene body region on gene expression is unclear, so it is worthy of future studies.

In our study, a total of 5832 DMCs and 6484 DMRs were identified. GO and KEGG enrichment analysis found that blue-eggshell egg-related pathways were strongly enriched, including glycine-gated chloride ion channel activity, and ABC transporters. Glycine is the raw material for the synthesis of biliverdin, and co-synthesizes porphobilinogen with CoA succinate. Research has shown that glycine binds to receptors and acts on glycine-gated chloride ion channels, resulting in increased chloride ion transport; at the same time, glycine receptors can enhance glycine transport through the activities of various second messengers in cells (Xu et al., 1996; Lynch et al., 1999), such as Ca²⁺. Meng’s research identified a number of genes related to Ca²⁺ and Cl^- transport through the transcriptome analysis of green duck eggshell glands (Meng, 2017), indicating that these genes can regulate eggshell color synthesis by affecting ion transport. Therefore, it is speculated that the methylation of the genes involved in the glycine-gated chloride ion channel activity pathway, including GRAL2 and GRAL3, may be activated by Ca²⁺, which further increases the transport of glycine and regulates the synthesis of biliverdin. Additionally, many studies have shown that members of the ABC transporter family can regulate blue eggshell color; in this study, it is speculated that it may play a role through the modification of DNA methylation.

Subsequently, we performed a combined analysis of RNA-seq data and WGBS data, and the results showed that DNA methylation levels and gene expression were highly negatively correlated in the CG context, which is identical to previous research. CpG methylation significantly alters the local DNA shape, which can potentially impact the accessibility of transcription factors to target sequences (Rao et al., 2018), and thus negatively regulate gene expression. This relationship was not observed in the CHG and CHH contexts. Venn diagram analysis found that 55 DEGs were regulated by methylation. Most genes were located in the gene body region, but only 3 genes were located in the promoter region.

The GO and KEGG enrichment analysis of these genes found that, among the enriched pathways, bile acid receptor signaling and activity, ferrous oxidase activity, and ABC transport pathway had a correlation with the formation of blue eggshell eggs, and the geneS NR1H4, HEPHL1, ABCA12, and ABCA4, were involved. This indicated that DNA methylation plays an important role in the regulation of the blue eggshell color in the Xuefeng Black-bone chickens. As the bile acid receptor in the bile acid metabolism pathway, NR1H4 plays an important role in bile excretion (Cai et al., 2007; Kovacs et al., 2008). HEPHL1 and HEPH may act synergistically in the oxidation of copper and iron ions in the porphyrin metabolic pathway, thereby participating in the synthesis of biliverdin (Chen et al., 2010). Therefore, it is speculated that DNA methylation modification regulated the expression of these genes, further affecting the synthesis of biliverdin and regulating the differences in blue eggshell colors.

CONCLUSION

In conclusion, our study first combined transcriptome sequencing and whole-genome methylation sequencing to analyze the molecular mechanisms regulating blue eggshell color variations. We uncovered several potential genes (such as NR1H4, HEPH1, ABCA4, and ABCA12) and pathways (such as intracellular bile acid receptor signaling pathway and ABC transporters) related to the blue eggshell. Our study will support the further study of molecular mechanisms in blue eggshell color variations in poultry. However, our work clearly has some limitations. Our mining and analysis of such huge RNA-seq and WGBS data are still not sufficient, thereby further exploration with bioinformatics analysis tools is needed in the future. Moreover, the screened potential genes and pathways need to be further functionally verified by experimental means.

ACKNOWLEDGMENTS

This study was funded by the undergraduate innovation training program (XCX2021017) and the graduate science innovation program (XCX202137) from Hunan Agricultural University, and the Department of Science and Technology of Hunan Province (2021NK4228). All the authors of this article contributed to the experimentation and writing process.

REFERENCE

• Akalin A, Kormaksson M, Li S, et al. MethylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biology 2012;13:R87. https://doi.org/10.1186/gb-2012-13-10-r87
  » https://doi.org/10.1186/gb-2012-13-10-r87
• Bai DP, Lin XY, Wu Y, et al. Isolation of blue-green eggshell pigmentation-related genes from Putian duck through RNA-seq. BMC Genomics 2019;20:66. https://doi.org/10.1186/s12864-019-5436-4
  » https://doi.org/10.1186/s12864-019-5436-4
• Bakken GS, Vanderbilt VC, Buttemer WA, et al. Avian eggs: thermoregulatory value of very high near-infrared reflectance. Science 1978;200(4339):321-3. https://doi.org/ 10.1126/science.200.4339.321.
  » https://doi.org/
• Ball MP, Li JB, Gao Y, et al. Targeted and genome-scale strategies reveal gene-body methylation signatures in human cells. Nature Biotechnology 2009;27(4):361-8. https://doi.org/10.1038/nbt.1533
  » https://doi.org/10.1038/nbt.1533
• Cai SY, Xiong LS, Wray CG, et al. The farnesoid X receptor FXRalpha/NR1H4 acquired ligand specificity for bile salts late in vertebrate evolution. American Journal of Physiology-Regulatory Integrative and Comparative Physiology 2007;293:R1400-1409. https://doi.org/10.1152/ajpregu.00781.2006
  » https://doi.org/10.1152/ajpregu.00781.2006
• Chen HJ, Attieh ZK, Syed BA, et al. Identification of zyklopen, a new member of the vertebrate multicopper ferroxidase family, and characterization in rodents and human cells. Journal of Nutrition. 2010;140:1728-35. https://doi.org/10.3945/jn.109.117531
  » https://doi.org/10.3945/jn.109.117531
• Chen JF, Hua GY, Han DP, et al. An EAV-HP insertion in the promoter region of SLCO1B3 has pleiotropic effects on chicken liver metabolism based on the transcriptome and proteome analysis. Scientific Reports 2021;11:7571. https://doi.org/10.1038/s41598-021-87054-9
  » https://doi.org/10.1038/s41598-021-87054-9
• Chen Q,Wang ZP. A new molecular mechanism supports that blue-greenish egg color evolved independently across chicken breeds. Poultry Science 2022;101(12):102223. https://doi.org/10.1016/j.psj.2022.102223
  » https://doi.org/10.1016/j.psj.2022.102223
• Chen SF, Zhou YQ, Chen Y, et al. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 2018;34:i884-i890. https://doi.org/10.1093/bioinformatics/bty560
  » https://doi.org/10.1093/bioinformatics/bty560
• Czubkowski P, Thompson RJ, Jankowska I, et al. Progressive familial intrahepatic cholestasis - farnesoid X receptor deficiency due to NR1H4 mutation: a case report. World Journal Clinical Cases 2021;9(15):3631-6. https://doi.org/10.12998/wjcc.v9.i15.3631
  » https://doi.org/10.12998/wjcc.v9.i15.3631
• Davit-Spraul A, Gonzales E, Jacquemin E. NR1H4 analysis in patients with progressive familial intrahepatic cholestasis, drug-induced cholestasis or intrahepatic cholestasis of pregnancy unrelated to ATP8B1, ABCB11 and ABCB4 mutations. Clinics and Research in Hepatology Gastroenterology 2012;36(6):569-73. https://doi.org/10.1016/j.clinre.2012.08.008
  » https://doi.org/10.1016/j.clinre.2012.08.008
• Elbirt KK, Bonkovsky HL. Heme oxygenase: recent advances in understanding its regulation and role. Proceedings of The Association of American Physicians 1999;111:438-47. https://onlinelibrary.wiley.com/doi/abs/10.1111/paa.1999.111.5.438
  » https://onlinelibrary.wiley.com/doi/abs/10.1111/paa.1999.111.5.438
• Ishikawa S, Suzuki K, Fukuda E, et al. Photodynamic antimicrobial activity of avian eggshell pigments. Febs Letters 2010;584:770-4. https://doi.org/10.1016/j.febslet.2009.12.041
  » https://doi.org/10.1016/j.febslet.2009.12.041
• Johannesson T, Kristinsson J, Torsdottir G, et al. Ceruloplasmin (Cp) and iron in connection with Parkinson's disease (PD) and Alzheimer's disease (AD). Laeknabladid. 2012;98(10):531-7. https://doi.org/10.17992/lbl.2012.10.457
  » https://doi.org/10.17992/lbl.2012.10.457
• Kilner RM. The evolution of egg colour and patterning in birds. Biological Reviews of the Cambridge Philosophical Society 2006;81:383-406. https://doi.org/10.1017/S1464793106007044
  » https://doi.org/10.1017/S1464793106007044
• Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nature Methods. 2015;12:357-60. https://doi.org/10.1038/nmeth.3317
  » https://doi.org/10.1038/nmeth.3317
• Kovacs P, Kress R, Rocha J, et al. Variation of the gene encoding the nuclear bile salt receptor FXR and gallstone susceptibility in mice and humans. Journal of Hepatology 2008;48:116-24. https://doi.org/10.1016/j.jhep.2007.07.027
  » https://doi.org/10.1016/j.jhep.2007.07.027
• Lang MR,Wells JW. A Review of eggshell pigmentation. World's Poultry Science Journal 1987;43:238-46. https://doi.org/10.1079/WPS19870016
  » https://doi.org/10.1079/WPS19870016
• Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature Methods 2012;9:357-9. https://doi.org/10.1038/nmeth.1923
  » https://doi.org/10.1038/nmeth.1923
• Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 2011;12:323. https://doi.org/10.1186/1471-2105-12-323
  » https://doi.org/10.1186/1471-2105-12-323
• Li ZJ, Ren TH, Li WY, et al. Association between the methylation statuses at CpG sites in the promoter region of the SLCO1B3, RNA expression and color change in blue eggshells in lushi chickens. Frontiers in Genetics 2019;10:161. https://doi.org/10.3389/fgene.2019.00161
  » https://doi.org/10.3389/fgene.2019.00161
• Lister R, Pelizzola M, Dowen RH, et al. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature 2009;462:315-22. https://doi.org/10.1038/nature08514
  » https://doi.org/10.1038/nature08514
• Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 2014;15:550. https://doi.org/10.1186/s13059-014-0550-8
  » https://doi.org/10.1186/s13059-014-0550-8
• Lynch JW, Han NI, Schofield PR. Building new function into glycine receptors: a structural model for the activation of the glycine-gated chloride channel. Clinical and Experimental Pharmacology and Physiology 1999;26:932-4. https://doi.org/10.1046/j.1440-1681.1999.03150.x
  » https://doi.org/10.1046/j.1440-1681.1999.03150.x
• Maines MD. The heme oxygenase system: a regulator of second messenger gases. Annual Review of Pharmacology and Toxicology 1997;37:517-54. https://doi.org/10.1146/annurev.pharmtox.37.1.517
  » https://doi.org/10.1146/annurev.pharmtox.37.1.517
• Meng GH. Exploring potential relationship underlying blue-green eggshell color of duck eggs association with candidate genes and major SLCO2B1 gene [dissertation]. Shaanxi Yangling (MG): Northwest A & F University; 2017. Available from: http://cdmd.cnki.com.cn/Article/CDMD-10712-1017101861.htm.
• Pranawidjaja S, Choi SI, Lay BW, et al. Analysis of heme biosynthetic pathways in a recombinant Escherichia coli. Journal of Microbiology and Biotechnology 2015;25(6):880-6. https://doi.org/10.4014/jmb.1411.11050
  » https://doi.org/10.4014/jmb.1411.11050
• Pertea M, Pertea GM, Antonescu CM, et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nature Biotechnology 2015;33:290-95. https://doi.org/10.1038/nbt.3122
  » https://doi.org/10.1038/nbt.3122
• Rao S, Chiu TP, Kribelbauer JF, et al. Systematic prediction of DNA shape changes due to CpG methylation explains epigenetic effects on protein-DNA binding. Epigenetics & Chromatin 2018;11(1):6. https://doi.org/10.1186/s13072-018-0174-4
  » https://doi.org/10.1186/s13072-018-0174-4
• Rauch TA, Wu X, Zhong X, et al. A human B cell methylome at 100-base pair resolution. Proceedings of The National Academy of Sciences of The United States of America 2009;106(3):671-8. https://doi.org/10.1073/pnas.0812399106
  » https://doi.org/10.1073/pnas.0812399106
• Robinson MD, McCarthy DJ, Smyth GK. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010;26:139-40. https://doi.org/10.1093/bioinformatics/btp616
  » https://doi.org/10.1093/bioinformatics/btp616
• Samiullah S, Roberts J, Chousalkar K. Oviposition time, flock age, and egg position in clutch in relation to brown eggshell color in laying hens. Poultry Science 2016;95:2052-7. https://doi.org/10.3382/ps/pew197
  » https://doi.org/10.3382/ps/pew197
• Samiullah S, Roberts JR, Chousalkar K. Eggshell color in brown-egg laying hens - a review. Poultry Science 2015;94:2566-75. https://doi.org/10.3382/ps/pev202
  » https://doi.org/10.3382/ps/pev202
• Shann YJ, Cheng C, Chiao CH, et al. Genome-wide mapping and characterization of hypomethylated sites in human tissues and breast cancer cell lines. Genome Research 2008;18(5):791-801. https://doi.org/10.1101/gr.070961.107
  » https://doi.org/10.1101/gr.070961.107
• Wang B, Chen XY. Comparative analysis of egg quality with different shell color. Journal of Hebei North University (Natural Science) 2020; 3633-36. http://dx.chinadoi.cn/10.3969/j.issn.1673-1492.2020.01.006
  » http://dx.chinadoi.cn/10.3969/j.issn.1673-1492.2020.01.006
• Wang XT, Zhao CJ, Li JY, et al. Comparison of the total amount of eggshell pigments in Dongxiang brown-shelled eggs and Dongxiang blue-shelled eggs. Poultry Science 2009;88:1735-1739. https://doi.org/10.3382/ps.2008-00434
  » https://doi.org/10.3382/ps.2008-00434
• Wang ZP, Qu LJ, Yao JF, et al. An EAV-HP insertion in 5' Flanking region of SLCO1B3 causes blue eggshell in the chicken. PLoS Genetics. 2013;9:e1003183. https://doi.org/10.1371/journal.pgen.1003183
  » https://doi.org/10.1371/journal.pgen.1003183
• Wierzbicka D, Gromadzka G. Ceruloplasmin, hephaestin and zyklopen: the three multicopper oxidases important for human iron metabolism. Postepy Higieny i Medycyny Doswiadczalnej 2014;68:912-24. https://doi.org/10.5604/17322693.1111136
  » https://doi.org/10.5604/17322693.1111136
• Wragg D, Mwacharo JM, Alcalde JA, et al. Endogenous retrovirus EAV-HP linked to blue egg phenotype in Mapuche fowl. PLoS One 2013;8:e71393. https://doi.org/10.1371/journal.pone.0071393
  » https://doi.org/10.1371/journal.pone.0071393
• Xi YX, Li W. BSMAP: whole genome bisulfite sequence MAPping program. BMC Bioinformatics 2009;10:232. https://doi.org/10.1186/1471-2105-10-232
  » https://doi.org/10.1186/1471-2105-10-232
• Xu FQ, Li A, Lan JJ, et al. Study of formation of green eggshell color in ducks through global gene expression. PLoS One 2018;13:e0191564. https://doi.org/10.1371/journal.pone.0191564
  » https://doi.org/10.1371/journal.pone.0191564
• Xu TL, Nabekura J, Akaike N. Protein kinase C-mediated enhancement of glycine response in rat sacral dorsal commissural neurones by serotonin. Journal of Physiology 1996;496(Pt 2):491-501. https://doi.org/10.1113/jphysiol.1996.sp021701
  » https://doi.org/10.1113/jphysiol.1996.sp021701
• Yuan QY, Lu LZ. Progresses in inheritance of genes for avian eggshell color. Hereditas 2007;29:265-68. https://doi.org/10.1360/yc-007-0265
  » https://doi.org/10.1360/yc-007-0265
• Zhang Y, Cheng RC, Guo R, et al. Clinical significance of serum protoporphyrin IX measurement in patients with HBV-related chronic liver diseases. Journal of Clinical Hepatology 2021;37:1075-80. http://dx.chinadoi.cn/10.3969/j.issn.1001-5256.2021.05.020
  » http://dx.chinadoi.cn/10.3969/j.issn.1001-5256.2021.05.020
• Zheng CW, Li ZS, Yang N, et al. Quantitative expression of candidate genes affecting eggshell color. Animal Science Journal 2014;85(5):506-10. https://doi.org/10.1111/asj.12182
  » https://doi.org/10.1111/asj.12182
• Zhu ZW, Zhang J, Wang Q, et al. Recent advances in the production of high-value products of porphyrin metabolism pathways and their microbial synthesis. Scientia Sinica Vitae 2020;50:1405-17. https://doi.org/10.1360/SSV-2020-0152
  » https://doi.org/10.1360/SSV-2020-0152

Publication Dates

• Publication in this collection
  26 Jan 2024
• Date of issue
  2024

History

• Received
  17 Oct 2022
• Accepted
  15 Nov 2023