Skip to main content

Advertisement

RNA-sequencing analysis of shell gland shows differences in gene expression profile at two time-points of eggshell formation in laying chickens

Abstract

Background

Eggshell formation takes place in the shell gland of the oviduct of laying hens. The eggshell is rich in calcium and various glycoproteins synthesised in the shell gland. Although studies have identified genes involved in eggshell formation, little is known about the regulation of genes in the shell gland particularly in a temporal manner. The current study investigated the global gene expression profile of the shell gland of laying hens at different time-points of eggshell formation using RNA-Sequencing (RNA-Seq) analysis.

Results

Gene expression profiles of the shell gland tissue at 5 and 15 h time-points were clearly distinct from each other. Out of the 14,334 genes assessed for differential expression in the shell gland tissue, 278 genes were significantly down-regulated (log2 fold change > 1.5; FDR < 0.05) and 413 genes were significantly up-regulated at 15 h relative to the 5 h time-point of eggshell formation. The down-regulated genes annotated to Gene Ontology (GO) terms showed anion transport, synaptic vesicle localisation, organic anion transport, secretion and signal release as the five most enriched terms. The up-regulated gene annotation showed regulation of phospholipase activities, alanine transport, transmembrane receptor protein tyrosine kinase signalling pathway, regulation of blood vessels diameter and 3, 5-cyclic nucleotide phosphodiesterase activity as the five most enriched GO terms. The putative functions of genes identified ranged from calcium binding to receptor activity. Validation of RNA-Seq results through qPCR showed a positive correlation.

Conclusions

The down-regulated genes at 15 h relative to the 5 h time-point were most likely involved in the transport of molecules and synthesis activities, initiating the formation of the eggshell. The up-regulated genes were most likely involved in calcium transportation, as well as synthesis and secretory activities of ions and molecules, reflecting the peak stage of eggshell formation. The findings in the current study improve our understanding of eggshell formation at the molecular level and provide a foundation for further studies of mRNA and possibly microRNA regulation involved in eggshell formation in the shell gland of laying hens.

Background

In laying hens, eggshell formation takes place in the shell gland region of the oviduct over an 18 h period [1]. The eggshell is composed of six distinct layers having calcite as a main component [2]. The calcium and bicarbonate ions contributing to the calcite (calcium carbonate) are secreted from the epithelial cells of the mucosa of the shell gland (reviewed in Hincke et al.) [3]. In addition to other regulators and transporters, the calbindin (CALB1) gene is involved in Ca2+ transportation across the cell membrane for eggshell formation [4]. Calcification within the shell gland is associated with stimuli initiated by ovulation or by neuroendocrine factors that control and coordinate both ovulation and calcium secretion [5]. Calcification of the egg first occurs slowly, increases to a rate of up to 300 mg/hr. over a duration of 15 h, and then again slows during the last 2 h before oviposition [5]. A higher rate of eggshell calcification may be correlated with a significantly higher level of calbindin mRNA expression that peaks at 16 h compared with 0–4.5 h of post-oviposition time [6]. While the expression level of the calbindin gene increases during the ovulatory cycle at a time coincident with eggshell calcification, there is little or no change in the tissue levels of calbindin protein, indicating post-translational control of calbindin levels [4]. Calcium secretion from the shell gland cells increases approximately 7 h after ovulation and reaches a maximum level as the shell is being formed [5]. The hormonal signals affecting changes in the rate of calcium secretion are not fully understood, although estrogen involvement has been suggested [7]. It is suggested that secretion of calcium from the shell gland cells may occur both by active transport and diffusion [8], involving expenditure of metabolic energy [5]. Calcium secretion appears to be linked functionally to luminal HCO3 concentration [8]. It seems that there is the involvement of a number of synthetic pathways in eggshell formation. About 37 ion transport genes have been shown to be involved in eggshell formation [9].

The organic components of the eggshell are shell membranes, mammillary cores, shell matrix, cuticle and pigment [10, 11]. The inorganic components of the eggshell are mammillary layer, palisade layer and surface crystal layer [10, 12]. More than 500 eggshell proteins have been identified in laying hens [13, 14]. All layers except for the shell membranes are formed in the shell gland. The membranes are composed of 10% collagen (types I, V and X) as well as 70–75% of other proteins and glycoprotein containing hexosamine and galactose [15,16,17,18,19] and lipids [20, 21]. The mammillary cores contain protein, carbohydrate and fat [22]. The eggshell matrix is a series of layers of protein and acid mucopolysaccharide [10, 23]. Some of the vital eggshell matrix proteins are ovocalyxin-36 [24, 25], ovocleidin-17 [26], ovocalyxin-32 [27] and ovocalyxin-25 [28] all of which possess antimicrobial functions. The cuticle is composed of glycoprotein (90%), polysaccharides (4%), lipids (3%) and inorganic phosphorus (3%) including hydroxyapatite crystals [10, 19, 29]. The major pigments of avian eggshells are protoporphyrin, zinc porphyrin, biliverdin and zinc biliverdin [30]. The calcified eggshell consists primarily of calcite, the most stable polymorph of calcium carbonate [10]. The metallo-proteinase family of proteins has been shown to play a role in reproductive tract remodelling [31]. Genes such as SPP1 (Secreted phosphoprotein 1), ACP1 (Acid phosphatase 1), PENK (Proenkephalin), RCAN1 (Regulator of calcineurin 1), CALB1 and CYP26A1 (Cytochrome P450 family 26 subfamily A member 1) have been shown to be actively involved in eggshell formation [32, 33]. However, global gene regulation in the eggshell formation of laying hens has not been reported.

We hypothesised that the regulation of genes involved in eggshell formation in the shell gland differs at different stages of eggshell formation on a global scale depending on the shell gland’s molecular and energetic requirements. To test this hypothesis we collected shell gland tissue for analysis when the egg was either forming in the distal magnum or isthmus and in the shell gland regions of the oviduct in brown-egg laying hens. Therefore, the main objective of the current study was to acquire a comprehensive picture of the transcriptional changes in the shell gland of brown-egg laying hens at two different time-points of eggshell formation. We also aimed to identify unknown candidate genes involved in eggshell formation.

Methods

Rearing of laying hens

Day old Isa-Brown laying chickens were obtained from the Baiada Hatchery at Tamworth, NSW, Australia. At the hatchery, day-old chickens received Rispens vaccine against Marek’s disease. The chickens were raised in isolation sheds at the University of New England under strict biosecurity protocols. All chickens were fed commercial starter to 4 weeks of age, pullet grower to 18 weeks of age and layer mash until the termination of the experiment. From the isolation sheds, pullets were moved at 18 weeks of age to individual cages in an isolated poultry house. At 35 weeks of age, eggs were collected and processed for traditional egg quality measurements following the method of Samiullah et al. [34]. Hens were then divided into a 1 × 2 factorial design in such a way that the egg weight and eggshell colour (L*) were not significantly different (P > 0.05) between the groups (Additional file 1: Table S1). Individual hen oviposition times were recorded by video camera, and each hen was processed at a specific post-oviposition time (5 and 15 h). At the time of euthanasia, the egg in individual hens was either in the distal magnum/isthmus (5 h post-oviposition time-point) or in the shell gland (15 h time-point).

Tissue collection

A total of forty hens were euthanised with CO2 gas and the shell gland was aseptically retracted through an abdominal incision. An approximately 500 mg sample tissue was cut from the centre of the shell gland and transferred directly to RNALater (Sigma Aldrich, Sydney, Australia). The samples were stored at − 20 °C and were processed for total RNA extraction within one day of collection. For total RNA extraction, a whole piece of shell gland tissue (all tissue layers) was processed.

Total RNA extraction and purification

Total RNA was extracted using TRIsure (Bioline, Australia), according to the manufacturer’s instructions. Briefly, an approximately 100 mg of tissue (wet weight) was homogenized in 1 mL of TRIsure using an IKA T10 basic Homogenizer (Wilmington, NC, USA). After the RNA pellets were washed with 1.5 mL ethanol (75%), 50 μL of UltraPure™ DEPC-treated water (Ambion, USA) was used to dissolve RNA pellets. The dissolved RNA was further purified using an RNeasy Mini Kit (Qiagen, GmbH, Hilden, Germany) as per the manufacturer’s instructions. A DNase-I step was performed to remove traces of genomic DNA from the extracted total RNA. The elution of RNA from the spin column with 50 μL of RNase-free water was repeated twice and the eluted RNA solutions were mixed thoroughly. The purified RNA was analysed in a NANODROP-8000 spectrophotometer (ThermoFisher Scientific, Wilmington, DE, USA) to measure its quantity and purity. The absorbance measurements of the spectrophotometer 260/280 and 260/230 ratios were in the range of 1.8–2.1 and 1.9–2.2, respectively. RNA integrity and purity were further examined in an Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany) as per the manufacturer’s instructions for an Agilent RNA 6000 Nano Kit. All the RNA showed distinct 18S and 28S bands with an average RNA integrity number (RIN) of > 9.1. Representative shell gland purified total RNA samples (Fig. 1) were processed by the Australian Genome Research Facility (AGRF) for RNA-Seq analysis.

Fig. 1
figure1

Schematic diagram explaining the selection of samples for RNA-Seq and differentially expressed genes (DEGs) data analysis. To validate the results of RNA-Seq, qPCR was performed on all 40 RNA samples and the data were analysed in qbase+ software for gene expression study. For bioinformatics analysis, the 5 h time-point was taken as the reference control as compared to the 15 h time-point

cDNA libraries preparation

Illumina’s TruSeq Stranded mRNA Prep Kit was used for processing the RNA samples. The process included mRNA purification via oligo (dT) beads, fragmentation of mRNA with divalent cations and heat, and 1st strand cDNA and 2nd strand cDNA syntheses. cDNA libraries were prepared by DNA fragment end repair, 3′ adenylation of DNA fragments, sequence adaptor ligation and amplification of library via PCR. In total, 12 cDNA libraries (i.e. one library for each sample) were constructed for sequencing - 6 samples for each of the time-points at 5 h and 15 h post-oviposition. Sequencing of libraries using 100 bp single read was performed on an Illumina HiSeq 2000 sequencing system.

Sequence quality control and sequence data evaluation

The primary sequence data were generated using the Illumina bcl2fastq 2.17.1.14 pipeline.

Initial quality control of the RNA sequences was evaluated by FastQC v0.11.5 [35]. The raw reads were also screened for the presence of any Illumina adaptor/overrepresented sequences, low quality sequences, empty reads and cross species contamination. Illumina adaptors and contaminated sequences were removed through trim_galore and Fastq_clean (https://ieeexplore.ieee.org/document/6999309/).

Reads mapping and raw gene counts

Tophat aligner (v2.0.14) [36] was used to map reads to the genomic sequences. The counts of reads mapping to each known gene were summarised at gene level using the featureCounts v1.4.6-p5 utility of the subread package (http://subread.sourceforge.net/). The cleaned sequence reads were aligned against the Gallus gallus genome (Built version 5 Ensembl release 86) [37].

Reference guided transcript assembly

The transcripts were assembled with Stringtie tool v1.2.4 (http://ccb.jhu.edu/software/stringtie/) utilising the reads alignment and reference annotation based assembly option (RABT). This option generated assembly for known and potentially novel transcripts. The Ensembl annotations (Gallus_gallus.Gallus_gallus-5.0.86.gtf) for genome were used as a guide. Common gene names were converted to Entrez IDs using Ensembl of chicken genome assembly.

Differential gene expression data analysis

Gene expression was calculated in counts-per-million (CPM) with a hard filter of 0.5 in edgeR (v3.14.0). Trimmed mean of M values (TMM) normalisation was applied to estimate gene expression and identify differentially expressed genes (DEGs) using R packages (R version 3.3.1) ‘edgeR’ [38] and ‘limma’ (3.28.21) [39]. During differential gene data analysis, false discovery rate/adjusted p-value was used for multiple test comparison (BH-adjustment). DEGs obtained at the 15 h time-point were compared to the 5 h time-point using the later as reference. To obtain further insight on the functions of the DEGs encoding hypothetical proteins, Ensembl BLAST/BLAT searches were performed with nucleotide and protein sequences queries using a cut-off e-value of < 10− 20.

Hierarchical clustering analysis

Average gene counts for the top 50 significantly down- or up-regulated genes at the 15 h relative to the 5 h time-point of eggshell formation were considered when performing the hierarchical clustering. The clustering was performed in gplots (version 3.0.1) of R packages (version 3.3.1), and the results were presented as heatmaps.

Functional annotation of DEGs

The DEGs (log2 fold change > 1.5; FDR < 0.05) were subjected to functional analysis using ClueGO version 2.2.6 [40, 41] + CluePedia version 1.2.6 [42] plugins in Cytoscape version 3.4.0 [43] as has previously been used in a similar study [44]. The DEGs were enriched for terms specific for biological process (BP), molecular functions (MF) and cellular component (CC). The annotation enrichment of the DEGs was performed with the 5 h time-point being considered as a reference control.

To create the annotation network, ClueGO investigates the distribution of the target genes across the Gene Ontology (GO) terms and pathways. CluePedia is a Cytoscape plugin for pathway insights using integrated experimental and in silico data [42]. CluePedia extends the functionality of ClueGO down to gene level [42]. In ClueGO analysis, the P value was calculated using the right-sided hypergeometric tests with Benjamini-Hochberg adjustment for multiple test correction [45]. An adjusted P ≤ 0.001 indicated a statistically significant deviation from the expected distribution, and that the corresponding GO terms and pathways were enriched for the target genes. The association strength between the terms was calculated using a corrected kappa statistic score of 0.4, in ClueGO [41, 46]. The relationship between the selected terms was defined based on their shared genes in a similar way. The created network showed the GO terms as nodes and size of the nodes reflected the enrichment significance. The network was automatically laid out using the organic layout algorithm supported by the Cytoscape software [43]. Functional groups represented by their most significant term were visualized in the network providing an insightful view of their interrelations [40].

Primer design, specificity and amplification efficiency for qPCR

Primers for the candidate target genes were designed in NCBI software by choosing an option for exon-intron spanning (Table 1). Primers for the reference and CALB1 genes were sourced from published literature. Specific amplifications of the primers were confirmed by a single peak of melting curve analysis and a single amplicon band of appropriate size using Agilent 2100 Bioanalyzer gel electrophoresis and a DNA 1000 Kit (Agilent Technologies, Waldbronn, Germany). PCR amplification efficiencies and correlation coefficients (R2) were determined with the amplifications of a series of six 10-fold RNA dilutions [47,48,49]. Amplification efficiency was calculated based on the following equation: \( E\ \left(\%\right)=\left({10}^{\frac{1}{- slope}}-1\right)\times 100 \) %. The primer pairs were used for qPCR analysis only if the amplification efficiency was in a range of 90 to 105%, and linear correlation coefficient R2 > 0.980 [49,50,51].

Table 1 Candidate target and reference* genes in expression studies by qPCR in the shell gland of laying hens

Quantitative PCR validation of RNA-Seq results

Quantitative PCR was performed on 40 samples of shell gland tissue RNA with the SensiFAST SYBR® Lo-ROX One-Step RT-PCR Kit (Bioline, Sydney, Australia). Master mix was prepared as per the manufacturer’s protocol and 4 μL of RNA template with 1:100 dilutions was added to the reaction wells using a QIAgility robotic (Qiagen, Hilden, Germany). The reaction was run in duplicates of 20 μL in a Rotor-gene Disc 100 (Qiagen, Hilden, Germany) with a Rotor-Gene Q thermocycler (Qiagen, Hilden, Germany). No template control (NTC) and no reverse transcriptase (−RT) control were also included to detect possible contamination. Thermocycling conditions for a 2-step PCR were: reverse transcription at 45 °C for 10 min, first denaturation at 95 °C for 2 min, then 40 cycles of denaturation at 95 °C for 5 s and annealing at appropriate temperatures (shown in Table 1) for 20s. The fluorescent data were acquired at the end of each annealing step during PCR cycles. A melting step was conducted to assess the specificity of PCR amplification.

Statistical analysis

Egg quality data were analysed by Statview software (SAS Institute Inc., Version 5.0.1.0). To calculate the relative expression of the candidate target genes, Cq values were analysed in qbase+ software version 3.0 [52]. The Cq values of target genes were normalized against previously optimised reference genes (TBP: TATA-Box binding protein and YWHAZ: Tyrosine 3-monooxygenase/Tryptophan 5-monooxygenase activation protein zeta) [53] to obtain normalized relative quantities of individual genes. Candidate target gene specific amplification efficiencies were used based on the method of Pfaffl [54]. The normalized relative quantities were further analysed in Statview software to compare the means from the time-points of 5 and 15 h. Significant differences were separated by the Tukey-Kramer test at probability < 0.05.

Results

Differential gene expression in shell gland tissue

A total of 261,684,549 (26.17 Gb of data bulk) clean reads with an average length of 100 bp were generated from the twelve libraries divided into two groups (G1 and G2; Fig. 1). The reads feature summary is depicted in Table 2. The feature summary shows that the percentage of reads mapped to Gallus gallus genome was ≥80%. Multi-dimensional scaling (MDS) plot showed that there was a significant effect of time-point on the expressed genes (Additional file 2: Figure S1).

Table 2 Sequence quality and alignment information of 12 shell gland samples in two groups (G1 and G2)

A total of 14,334 gene transcripts were assessed for differential expression after filtering was applied. Differential gene expression analysis showed 691 (log2 fold change > 1.5; FDR < 0.05) differentially expressed genes (DEGs) between the 5 h time-point and the 15 h time-point of eggshell formation. Among the 691 DEGs, there were 278 significantly down-regulated and 413 significantly up-regulated genes at the 15 h time-point relative to the 5 h time-point of eggshell formation. Among the DEGs at the 15 h relative to the 5 h time-point, SLC13A5 (Solute carrier family 13 member 5), KLB (Klotho beta), XAF1 (XIAP associated factor 1), FIBIN (Fin bud initiation factor homolog (zebrafish)) and POMGNT1 (Protein O-linked mannose N-acetylglucosaminyltransferase 1 (Beta 1,2-) were the top five most down-regulated genes. A full list of the top 50 significantly down-regulated genes at the 15 h relative to the 5 h time-point is shown in Table 3. Among the DEGs that were significantly up-regulated at the 15 h relative to the 5 h time-point, the top five genes were POMC (Proopiomelanocortin5), CALB1 (Calbindin), SPP1 (secreted phosphoprotein 1), NEU4 (Neuraminidase 4) and CEMIP (Cell migration inducing hyaluronan binding protein). A full list of the top 50 significantly up-regulated genes at the 15 h relative to the 5 h time-point is shown in Table 4.

Table 3 Top 50 down-regulated DEGs at 15 h relative to 5 h time-point
Table 4 Top 50 up-regulated DEGs at 15 h relative to 5 h time-point

DEGs analysis for hypothetical functions

Most of the DEGs with hypothetical functions appeared to possess domains that function in diverse cellular activities (Table 5). The associated GO terms showed that the functions of the unknown genes may be correlated with the synthesis and secretory activities in the shell gland during an eggshell formation. In addition, there were 6.11 and 6.31% of lincRNA significantly (log2 fold change > 1.5; FDR < 0.05) down- and up-regulated, respectively, at the 15 h relative to the 5 h time-point of eggshell formation.

Table 5 Putative functions of mRNA sequences that did not annotate to any known gene IDs during Ensembl annotations of RNA-Seq transcripts

Functional annotation of DEGs down- or up-regulated at 15 h relative to 5 h time-point

An enrichment gene set analysis was performed to identify the associated Gene Ontology (GO) terms specific to Biological Process (BP), Cellular Component (CC) and Molecular Functions (MM). A total of 278 genes (log2 fold change > 1.5; FDR < 0.05) significantly down-regulated at the 15 h relative to the 5 h time-point were mapped to the GO terms specific for BP, CC and MF pathways. The most enriched GO terms associated with DEGs are depicted in Fig. 2. Out of the 14 GO terms revealed, the five major terms associated with the down-regulated genes were anion transport (GO:0006820), synaptic vesicle localization (GO:0097479), organic anion transport (GO:0015711), secretion (GO:0046903) and signal release (GO:0023061) (Fig. 2a). It should be noted that all of the 14 GO terms were significantly enriched (enrichment pathway P value < 0.05) (Fig. 2a). For the functional analysis of genes significantly up-regulated at the 15 h relative to the 5 h time-point, a total of 413 genes were mapped to the GO terms specific for BP, CC and MF pathways. All of the GO terms enriched were significant at an enrichment pathway P value < 0.05 (Fig. 2b). Out of the total 10 GO terms, the five major terms were: regulation of phospholipase activity (GO:0010517), alanine transport (GO:0032328), transmembrane receptor protein tyrosine kinase signaling pathway (GO:0007169), regulation of blood vessel diameter (GO:0097746) and 3′,5′-cyclic-nucleotide phosphodiesterase activity (GO:0004114). Network representation of the enriched GO terms and their associated genes obtained from the mapped genes down-regulated at the 15 h relative to the 5 h time-point is depicted in Fig. 3. Network representation of the enriched GO terms and their associated genes obtained from the mapped genes up-regulated at the 15 h relative to the 5 h time-point is depicted in Fig. 4.

Fig. 2
figure2

Functional map of DEGs (log2 fold change > 1.5; FDR < 0.05) enriched for GO terms specific for biological process, cellular component and molecular function. The chart fragments represent the number of genes associated with the terms as a proportion with the total number of genes within the GO term. a GO terms associated with genes significantly down-regulated at the 15 h relative to the 5 h time-point of eggshell formation. b GO terms associated with genes significantly up-regulated at the 15 h relative to the 5 h time-point of eggshell formation. **P < 0.001 and *P < 0.01 show the level of significance of the enriched terms

Fig. 3
figure3

Network representation of the enriched GO terms and their associated genes obtained from the mapping of down-regulated genes at the 15 h relative to the 5 h time-point. The GO terms were identified as nodes and linked based on their p-value < 0.05 and kappa score level (> 0.4). Functionally related groups partially overlapped. The terms are labelled in colours according to hierarchical clustering of GO terms. Terms which have not been grouped are shown in grey

Fig. 4
figure4

Network representation of the enriched GO terms and their associated genes obtained from the mapping of up-regulated genes at the 15 h relative to the 5 h time-point. The GO terms were identified as nodes and linked based on their p-value < 0.05 and kappa score level (> 0.4). Functionally related groups partially overlapped. The terms are labelled in colours according to hierarchical clustering of GO terms. Terms which have not been grouped are shown in grey

Hierarchical clustering analysis

Hierarchical clustering analysis (HCA) was performed using the top 50 DEGs down- or up-regulated genes at the 15 h relative to the 5 h time-point of eggshell formation. The pattern of expression for the top 50 DEGs is presented in Fig. 5a, b. A clear difference for the pattern of DEGs at the two time-points has been visualised.

Fig. 5
figure5

Hierarchical clustering analysis of the top 50 DEGs significantly down- or up-regulated at the 15 h relative to the 5 h time-point. a Top 50 DEGs significantly down-regulated at the 15 h relative to the 5 h time-point. b Top 50 DEGs significantly up-regulated at the 15 h relative to the 5 h time-point. G1a-f represent six samples taken at the 5 h time-point, while G1a-f represent six samples taken at the 15 h time-point of eggshell formation

Validation of RNA-Seq data by qPCR

Quantitative PCR was performed to validate the significantly down- or up-regulated genes at the 15 h relative to the 5 h time-point obtained in RNA-Seq analysis. All primers used for RNA-Seq data validation by qPCR were specific in amplifications (Fig. 6a, b). The amplification efficiencies of individual primers have been depicted in Table 1. The expression levels of nineteen genes selected for validation of RNA-Seq data showed a positive linear relationship (Table 6). The results suggested that the RNA-Seq is a good reference for expression profiling study and the assembly quality of the sequences was desirable. Although the magnitude of fold change obtained by qPCR and RNA-Seq was slightly different, the qPCR results demonstrated a similar trend (positive correlation) compared with the RNA-Seq for the 19 genes being tested (Table 6).

Fig. 6
figure6

DNA gel electrophoresis of the qPCR products showing that the primers were specific in amplification. Panel a: L) DNA ladder; 1) GAL3ST2 (177 bp); 2) CYP7A1 (171 bp); 3) CALB1 (116 bp); 4) TYRO3 (186 bp); 5) CNG4 (243 bp); 6) POMC (115 bp); 7) RHOBTB3 (171 bp); 8) KCNH1 (160 bp); 9) MMP13 (128 bp); 10) KLB (224 bp); 11) GATA3 (157 bp); 12) SPP1 (134 bp). Panel b: L) DNA ladder; 1) YWHAZ (61 bp); 2) TBP (147 bp); 3) IBV T-as positive control (181 bp); 4) CA9 (80 bp); 5) OTOP2 (202 bp); 6) CGA (166 bp); 7) CLDN16 (186 bp); 8) GJA8 (186 bp); 9) SS2 (160 bp); 10) BPIFB3 (234 bp); 11) PPARGCIB-Positive control (82 bp); 12) TLR3-Positive control (203 bp). The upper (purple) and lower markers (green) act as internal standards and are used to align the ladder analysis with the individual DNA sample analysis. The DNA gel in Agilent 2100 Bioanalyzer was performed as per the manufacturer’s instructions of Agilent DNA 1000 Kit. The size of the individual amplicons are consistent with the expected size

Table 6 Comparison of the gene expression data between RNA-Seq and qPCR

All the genes tested were either significantly down or up-regulated (P < 0.05) at the 15 h relative to the 5 h time-point of eggshell formation. Regression analysis showed a weak positive correlation (R2 = 0.526; P value 0.004) between the qPCR and RNA-Seq data.

Discussion

Significant advances have been made in understanding the morphological and biochemical aspects of eggshell biogenesis. However, the molecular mechanisms underpinning the formation of various layers of the eggshell formation are still not well understood. The present study focused on how the regulation of genes was related to eggshell formation by the study of DEGs between the time points when the egg was either in the distal magnum or isthmus (5 h time-point, post oviposition time) or in the shell gland (15 h time-point, post oviposition time). For simplicity of data presentation, the 5 h time-point was taken as reference control to examine the expression changes of the genes at 15 h time-point when the eggshell formation was already initiated. Quantitative PCR results validated RNA-Seq data; therefore, RNA-Seq was used for genome-wide exploration of the gene expression profile of the shell gland. The RNA-Seq analysis revealed many DEGs down- or up-regulated at the 15 h relative to the 5 h time-point of eggshell formation. Some of the genes identified in the current study have been previously implicated in eggshell formation [55]; however, we have also identified multiple new genes that potentially play vital roles during active stages of eggshell formation. In addition, the current study has picturised the expression profile of shell gland when the egg was either in the distal magnum or the distal isthmus, reflecting the preparatory molecular mechanisms occurring in the shell gland. Various layers of eggshell result from the deposition of organic matrix and inorganic minerals secreted to the lumen of shell gland. In the current study, elaborating on molecular mechanisms occurring during eggshell formation; significantly up-regulated genes, such as CALB1, POMC, SPP1, BPIFB3 and down-regulated genes, such as SLC13A5, KLB, XAF1 and MMP13 reflect the differential expression profile of the shell gland during eggshell formation.

DEGs that were significantly up-regulated at the 15 h relative to the 5 h time-point and enriched for GO term pathway analysis showed active stages of eggshell formation. The GO term regulation of phospholipase activity (GO:0010517) shows that the hydrolysis of lipids was higher in order to produce energy for the synthetic processes of eggshell formation. Among the DEGs in phospholipase activity, CEMIP, ANGPTL3, WNT11, EREG, MAP 3 K15 and SLC20A1 were significantly up-regulated with log2 fold changes of 7.555, 6.198, 4.799, 3.867, 3.736 and 3.302, respectively. CEMIP is mainly involved in metabolism, glycosaminoglycan and calcium release metabolism pathways. CEMIP interacts with BIP/HSPA5 for the release of calcium from endoplasmic reticulum [56]. The higher expression levels of CEMIP and HSPA5 (log2 fold change 2.403) might indicate their role in calcium release for peak stages of eggshell formation.

ANGPTL3 is a member of angiopoietin-like (ANGPTL) genes that have diverse functions in various pathophysiological [57] and developmental [58] conditions in mammals. The N terminal chain of ANGPTL3 is also important for lipid metabolism. A higher mRNA expression of ANGPTL3 was observed in mouse uterus on day 6.5 of pregnancy [59]. In the chicken oviduct, a higher expression of ANGPTL3 was linked with molecular mechanisms involved in tissue development and remodelling [60]. A significantly higher expression of ANGPTL3 at the 15 h time-point shows its direct role in eggshell formation. It seems that ANGPTL3 might have been up-regulated by the release of endocrine hormones involved in molecular mechanisms of eggshell formation and oviposition. PTN is among the estrogen stimulating genes, possesses antimicrobial properties [55] and expresses in chicken oviduct [61]. The current study confirms a significant up-regulation of PTN during active stages of eggshell formation. The WNT11 gene functions in developmental processes and its up-regulation at the 15 h time-point (log2 fold change 4.799) compared with the 5 h time-point reflects its role in the peak/active stages of eggshell formation. WNT11 was up-regulated during eggshell formation in laying hens observed in other study [55]. In sheep uterus, WNT family encodes signalling regulator molecules vital for cell growth, differentiation and cell-cell interactions [62].

At the 15 h time-point, among the other significantly up-regulated genes were CALB1, POMC, SPP1, BPIFB3/OCX-36, LOC415478, KCNH1, BPIL3 and OTOP3 that have previously been implicated in eggshell formation [55]. A significant up-regulation of CALB1 at the 15 h (log2 fold change 8.081) relative to the 5 h time-point confirms a higher rate of calcium transportation across the cell membrane during the peak stages of eggshell formation. A higher expression of CALB1 during eggshell calcification in the shell gland and in the intestine of chickens has been reported [4, 33, 55, 63]. Low free Ca+ in cells is maintained by calcium uptake in the endoplasmic reticulum through ATP dependent calcium pumps [55]. ATP2A3 appears to play a role in this Ca+ balance, which is confirmed by its up-regulation (log2 fold change 3.484) at the 15 h relative to the 5 h time-point. SPP1 is another important gene involved in eggshell calcification [55, 64]. The peak stages of eggshell formation can be further linked with the significant higher expression of SPP1 (log2 fold change 7.993). A significantly higher expression of SPP1 was observed between 3 and 20 h post-oviposition times in the shell gland of laying hens by Jeong et al [33]. SPP1 is involved in bone mineralisation and is present in chicken eggshell [65, 66]. The expression of SPP1 in chicken uterine tissue is stimulated by the mechanical presence of the forming egg [33, 67]. The gene POMC functions in many biological pathways including the stimulation of the release of cortisol hormone. A significant up-regulation (log2 fold change 9.179) of the POMC at the 15 h time-point highlights its role in the release of hormones necessary for formation of eggshell. A higher expression of POMC was observed when a hard shell egg was forming in hen uterine tissue [55].

BPIFB3 (OCX-36) is a lipopolysaccharide-binding protein/bactericidal-permeability increasing protein (LBP/BPI) that is present in various layers of the eggshell and possesses antibacterial activity [25, 55, 68]. In the oviduct of laying chickens, OCX-36 only expresses in the shell gland [25]. In the current study, a significantly higher expression of OCX-36 (log2 fold change 6.413) at the 15 h vs the 5 h time-point indicates the importance of OCX-36 protein in the shell matrix. A higher expression level of OCX-36 mRNA has been shown in chicken shell gland in the presence of an egg [25, 55].

The second most enriched GO term at the 15 h time-point, alanine transport (GO:0032328), indicates that the alanine and 2-aminopropanoic acid transport across the shell gland cells was higher, which might be involved in energy (ATP) production during eggshell formation. Significantly up-regulated SLC6A17 (log2 fold change 6.642) at the 15 h relative to the 5 h time-point indicates the importance of the alanine transport pathway during eggshell formation. The transmembrane receptor protein tyrosine kinase signaling pathway (GO:0007169) is usually initiated by the binding of an extracellular ligand to a receptor on the surface of the target cells where the receptor possesses tyrosine kinase activity to regulate transcription. The most enriched GO:0007169 indicates higher transcriptional activities of the cells involved in the synthesis and secretion of macromolecules needed for eggshell formation. The GO enriched term regulation of blood vessel diameter (GO:0097746) suggests that the blood flow to the shell gland at the 15 h time-point was significantly affected by the eggshell formation as has been shown previously [69]. The genes involved in GO term cyclic 3′,5′-phosphodiesterase activity (GO:0004114) encode enzymes that degrade the phosphodiester bond in cAMP and cGMP molecules. The up-regulation of these genes in the shell gland at the 15 h relative to the 5 h time-point indicate their role in energy production during the synthetic activities of shell gland for eggshell formation.

Genes that were significantly down-regulated at the 15 h relative to the 5 h time-point reflect the activities of the shell gland when the egg was forming either in distal magnum or isthmus and was ready to enter to shell gland in the next hour or so. The down-regulated genes annotated to the most enriched GO term anion transport (GO:0006820) indicate that the genes involved in transport of ions across cell membrane were significantly down-regulated in the shell gland. This indicates that the synthesis and secretory activities in the shell gland cells were already initiated, while the egg was still forming in the distal magnum or isthmus. SLC13A5 (also known as Na+/citrate cotransporter) plays an important role in transporting ions and/or molecules across cell membranes. A significantly lower log2 fold change (− 6.516) of SLC13A5 at the 15 h relative to the 5 h time-point might reflect its role in transportation of ions for the initiation of synthesis of molecules necessary for the initiation of eggshell formation. The genes SLC13A5 and SLC13A2 belong to solute carrier family 13 group of proteins and are sodium-dependent citrate cotransporters in regulating metabolic processes. Among its related pathways are transport of various sugars, bile salts and organic acids, metal ions and amine compounds. In mammalian cells, SLC13A5 mediates Na+-coupled transport of citrate and succinate for tricarboxylic acid cycle [70]. In the GO term synaptic vesicle localisation, most of the genes involved function in transportation of synaptic vesicles across cell membrane. It seems that the genes in this pathway mainly perform activities in neurotransmission necessary for the transport and synthesis of various molecules including hormones in the shell gland as shown in the present study. Some of the genes that were annotated to the third most enriched GO term, organic anion transport, also served as transporters for organic anions across cell membrane. Organic anions contain molecules that are negatively charged and contain carbon in covalent linkage. The significantly enriched GO term secretion indicates the synthesis of substances that were either directly involved in eggshell formation or served a role in transportation of other molecules such as hormones. The enriched GO term signal release indicates that signal secretion to the extracellular medium from a cellular source was occurring around the 5 h time-point. This may indicate that the shell gland cells were actively involved in the synthesis of molecules necessary for either cellular function or initiation of eggshell formation.

The alignment of the sequences with unknown gene/protein functions suggests that these genes are vital to shell gland function in laying chickens. The majority of the significantly up-regulated DEGs with unknown functions were from the 15 h time-point. It seems that these DEGs were involved in the molecular mechanisms necessary for eggshell formation. We suggest further investigation of their roles in the shell gland relative to egg formation. The associated GO terms with the unknown function genes ranged from calcium ion binding to receptor activity. A large number of novel lincRNA in the current study might indicate their role as regulators in the shell gland of laying hens. Further studies should be performed to investigate the spatio-temporal expression of genes involved in the synthesis of various eggshell layers and the role of microRNA and lincRNA in the regulation of genes involved in eggshell formation.

Conclusions

Transcriptome analysis revealed thousands of DEGs in shell gland of laying chickens at the 15 h relative to the 5 h time-point of eggshell formation. The significantly down-regulated DEGs indicate that the synthesis activities were already initiated in the shell gland when the egg was still forming in the distal magnum or isthmus regions of the oviduct. The DEGs significantly up-regulated at the 15 h relative to the 5 h time-point reflect the phospholipid activities and synthesis or transport of molecules for the peak period of eggshell formation. The findings in the current study improve our understanding of eggshell formation at molecular level.

Abbreviations

ACP1:

Acid phosphatase 1

CALB1:

Calbindin

CYP26A1:

Cytochrome P450 family 26 subfamily A member 1

DEGs:

Differentially expressed genes

GO:

Gene ontology

PENK:

Proenkephalin

RCAN1:

Regulator of calcineurin 1

RIN:

RNA integrity number

SPP1:

Secreted phosphoprotein 1

TBP:

TATA-Box binding protein

YWHAZ:

Tyrosine 3-monooxygenase/Tryptophan 5-monooxygenase activation protein zeta

References

  1. 1.

    Nys Y. Relationships between age, shell quality and individual rate and duration of shell formation in domestic hens. Br Poult Sci. 1986;27(2):253–9.

  2. 2.

    Bain M. Egg composition and chemistry. In: Roberts J, editor. Achieving sustainable production of eggs. UK: Burleigh Dodds sScience pPublishing lLimited; 2017. p. 3–23.

  3. 3.

    Hincke MT, Nys Y, Gautron J, Mann K, Rodriguez-Navarro AB, McKee MD. The eggshell: structure, composition and mineralization. Front Biosci. 2012;17:1266–80.

  4. 4.

    Nys Y, Mayel-Afshar S, Bouillon R, Van Baelen H, Lawson DEM. Increases in calbindin D 28K mRNA in the uterus of the domestic fowl induced by sexual maturity and shell formation. Gen Comp Endocrinol. 1989;76(2):322–9.

  5. 5.

    Johnson AL. Reproduction in the female. In: Scanes CG, editor. Avian physiology. UK: Academic; 2015. p. 635–65.

  6. 6.

    Yang JH, Zhao ZH, Hou JF, Zhou ZL, Deng YF, Dai JJ. Expression of TRPV6 and CaBP-D28k in the egg shell gland (uterus) during the oviposition cycle of the laying hen. Br Poult Sci. 2013;54(3):398–406.

  7. 7.

    Eastin WC Jr, Spaziani E. On the control of calcium secretion in the avian shell gland (uterus). Biol Reprod. 1978;19(3):493–504.

  8. 8.

    Eastin WC Jr, Spaziani E. On the mechanism of calcium secretion in the avian shell gland (uterus). Biol Reprod. 1978;19(3):505–18.

  9. 9.

    Jonchère V, Brionne A, Gautron J, Nys Y. Identification of uterine ion transporters for mineralisation precursors of the avian eggshell. BMC Physiol. 2012;12(1):10.

  10. 10.

    Dennis JE, Xiao SQ, Agarwal M, Fink DJ, Heuer AH, Caplan AI. Microstructure of matrix and mineral components of eggshells from white Leghorn chickens (Gallus gallus). J Morphol. 1996;228(3):287–306.

  11. 11.

    Nys Y, Zawadzki J, Gautron J, Mills A. Whitening of brown-shelled eggs: mineral composition of uterine fluid and rate of protoporphyrin deposition. Poult Sci. 1991;70(5):1236–45.

  12. 12.

    Mikšík I, Charvátová J, Eckhardt A, Zk D. Insoluble eggshell matrix proteins-their peptide mapping and partial characterization by capillary electrophoresis and high-performance liquid chromatography. Electrophoresis. 2003;24(5):843–52.

  13. 13.

    Mikšík I, Eckhardt A, Sedláková P, Mikulikova K. Proteins of insoluble matrix of avian (Gallus gallus) eggshell. Connect Tissue Res. 2007;48(1):1–8.

  14. 14.

    Mann K, Maček B, Olsen JV. Proteomic analysis of the acid‐soluble organic matrix of the chicken calcified eggshell layer. Proteomics. 2006;6(13):3801–10.

  15. 15.

    Harris ED, Blount JE, Leach RM. Localization of lysyl oxidase in hen oviduct: implications in egg shell membrane formation and composition. Science. 1980;208(4439):55–6.

  16. 16.

    Wong M, Hendrix MJ, von der Mark K, Little C, Stern R. Collagen in the egg shell membranes of the hen. Dev Biol. 1984;104(1):28–36.

  17. 17.

    Arias JL, Fernandez MS, Dennis JE, Caplan AI. Collagens of the chicken eggshell membranes. Connect Tissue Res. 1991;26(1–2):37–45.

  18. 18.

    Fernandez MS, Araya M, Arias JL. Eggshells are shaped by a precise spatio-temporal arrangement of sequentially deposited macromolecules. Matrix Biol. 1997;16(1):13–20.

  19. 19.

    Baker JR, Balch DA. A study of the organic material of hen’s-egg shell. Biochem J. 1962;82(2):352–61.

  20. 20.

    Hasiak R, Vadehra D, Baker R. Lipid composition of the egg exteriors of the chicken (Gallus gallus). Comp Biochem Physiol. 1970;37(3):429–35.

  21. 21.

    Hasiak R, Vadehra D, Baker R. Fatty acid composition of the egg exterior structures of Gallus gallus. Comp Biochem Physiol. 1970;35(3):751–5.

  22. 22.

    Simkiss K. Tyler C. A histochemical study of the organic matrix of hen egg-shells. J Cell Sci. 1957;98(3):19–28.

  23. 23.

    Nys Y, Gautron J, Garcia-Ruiz JM, Hincke MT. Avian eggshell mineralization: biochemical and functional characterization of matrix proteins. Comptes Rendus Palevol. 2004;3(6–7):549–62.

  24. 24.

    Kovacs-Nolan J, Cordeiro C, Young D, Mine Y, Hincke M. Ovocalyxin-36 is an effector protein modulating the production of proinflammatory mediators. Vet Immunol Immunopathol. 2014;160(1):1–11.

  25. 25.

    Gautron J, Murayama E, Vignal A, Morisson M, McKee MD, Réhault S, Labas V, Belghazi M, Vidal M-L, Nys Y, et al. Cloning of ovocalyxin-36, a novel chicken eggshell protein related to lipopolysaccharide-binding proteins, bactericidal permeability-increasing proteins, and plunc family proteins. J Biol Chem. 2007;282(8):5273–86.

  26. 26.

    Mann K, Siedler F. The amino acid sequence of ovocleidin 17, a major protein of the avian eggshell calcified layer. IUBMB Life. 1999;47(6):997–1007.

  27. 27.

    Gautron J, Hincke MT, Mann K, Panhéleux M, Bain M, McKee MD, Solomon SE, Nys Y. Ovocalyxin-32, a novel chicken eggshell matrix protein isolation, amino acid sequencing, cloning, and immunocytochemical localization. J Biol Chem. 2001;276(42):39243–52.

  28. 28.

    Marie P, Labas V, Brionne A, Harichaux G, Hennequet-Antier C, Nys Y, Gautron J. Quantitative proteomics and bioinformatic analysis provide new insight into protein function during avian eggshell biomineralization. J Proteome. 2015;113:178–93.

  29. 29.

    Wedral EM, Vadehra DV, Baker RC. Chemical composition of the cuticle, and the inner and outer shell membranes from eggs of Gallus gallus. Comp Biochem Physiol B. 1974;47(3):631–40.

  30. 30.

    Kennedy G, Vevers H. A survey of avian eggshell pigments. Comp Biochem Physiol B. 1976;55(1):117–23.

  31. 31.

    Huh M-I, Jung J-C. Expression of matrix metalloproteinase-13 (MMP-13) in the testes of growing and adult chicken. Acta Histochem. 2013;115(5):475–80.

  32. 32.

    Lim W, Jeong W, Kim J, Ka H, Bazer FW, Han JY, Song G. Differential expression of secreted phosphoprotein 1 in response to estradiol-17β and in ovarian tumors in chickens. Biochem Biophys Res Commun. 2012;422(3):494–500.

  33. 33.

    Jeong W, Lim W, Kim J, Ahn SE, Lee HC, Jeong J-W, Han JY, Song G, Bazer FW. Cell-specific and temporal aspects of gene expression in the chicken oviduct at different stages of the laying cycle. Biol Reprod. 2012;86(6):172.

  34. 34.

    Samiullah S, Omar AS, Roberts J, Chousalkar K. Effect of production system and flock age on eggshell and egg internal quality measurements. Poult Sci. 2016;96(1):246–58.

  35. 35.

    Andrews S. FastQC: a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc. Accessed Nov 20 2016.

  36. 36.

    Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

  37. 37.

    Hillier LW, Miller W, Birney E, Warren W, Hardison RC, Ponting CP, Bork P, Burt DW, Groenen MAM, Delany ME. Sequence and comparative analysis of the chicken genome provide unique perspectives on vertebrate evolution. Nature. 2004;432(7018):695–716.

  38. 38.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

  39. 39.

    Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(4):e47.

  40. 40.

    Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman W-H, Pagès F, Trajanoski Z, Galon J. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25(8):1091–3.

  41. 41.

    Mlecnik B, Galon J, Bindea G. Comprehensive functional analysis of large lists of genes and proteins. J Proteome. 2018;171:2–10.

  42. 42.

    Bindea G, Galon J, Mlecnik B. CluePedia Cytoscape plugin: pathway insights using integrated experimental and in silico data. Bioinformatics. 2013;29(5):661–3.

  43. 43.

    Cline MS, Smoot M, Cerami E, Kuchinsky A, Landys N, Workman C, Christmas R, Avila-Campilo I, Creech M, Gross B. Integration of biological networks and gene expression data using Cytoscape. Nat Protoc. 2007;2(10):2366–82.

  44. 44.

    Hamzić E, Kjærup RB, Mach N, Minozzi G, Strozzi F, Gualdi V, Williams JL, Chen J, Wattrang E, Buitenhuis B. RNA sequencing-based analysis of the spleen transcriptome following infectious bronchitis virus infection of chickens selected for different mannose-binding lectin serum concentrations. BMC Genomics. 2016;17(1):82.

  45. 45.

    Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57(1):289–300.

  46. 46.

    Bindea G, Mlecnik B. ClueGO documentation. http://www.ici.upmc.fr/cluego/ClueGODocumentation.pdf. Accessed Jan 20 2017.

  47. 47.

    Rasmussen R. Quantification on the LightCycler. In: Rapid cycle real-time PCR. Berlin: Springer; 2001. p. 21–34.

  48. 48.

    Hellemans J, Mortier G, De Paepe A, Speleman F, Vandesompele J. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 2007;8(2):R19.

  49. 49.

    Bio-Rad Application User Guide. http://www.bio-rad.com/webroot/web/pdf/lsr/literature/Bulletin_5279.pdf. Accessed May 20 2016.

  50. 50.

    Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, Mueller R, Nolan T, Pfaffl MW, Shipley GL. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009;55(4):611–22.

  51. 51.

    Ginzinger DG. Gene quantification using real-time quantitative PCR: An emerging technology hits the mainstream. Exp Hematol. 2002;30(6):503–12.

  52. 52.

    qbase+ qPCR data analysis software version 3.0: https://www.qbaseplus.com/. Accessed 1 May 2017.

  53. 53.

    Khan S, Roberts J, Wu S-B. Reference gene selection for gene expression study in shell gland and spleen of laying hens challenged with infectious bronchitis virus. Sci Rep. 2017;7(1):14271.

  54. 54.

    Pfaffl MW. A new mathematical model for relative quantification in real-time RT–PCR. Nucleic Acids Res. 2001;29(9):e45.

  55. 55.

    Brionne A, Nys Y, Hennequet-Antier C, Gautron J. Hen uterine gene expression profiling during eggshell formation reveals putative proteins involved in the supply of minerals or in the shell mineralization process. BMC Genomics. 2014;15(1):220.

  56. 56.

    Evensen NA, Kuscu C, Nguyen H-L, Zarrabi K, Dufour A, Kadam P, Y-j H, Pulkoski-Gross A, Bahou WF, Zucker S. Unraveling the role of KIAA1199, a novel endoplasmic reticulum protein, in cancer cell migration. J Natl Cancer Inst. 2013;105(18):1402–16.

  57. 57.

    Kikuchi R, Tsuda H, Kozaki K-i, Kanai Y, Kasamatsu T, Sengoku K, Hirohashi S, Inazawa J, Imoto I. Frequent inactivation of a putative tumor suppressor, angiopoietin-like protein 2, in ovarian cancer. Cancer Res. 2008;68(13):5067–75.

  58. 58.

    Lai DM, Tu YK, Hsieh YH, Hsu WM, Lee CC, Cheng WC, Hsieh FJ, Li H. Angiopoietin-like protein 1 expression is related to intermuscular connective tissue and cartilage development. Dev Dyn. 2007;236(9):2643–52.

  59. 59.

    Scott CA, van Huyen D, Bany BM. Angiopoietin-like gene expression in the mouse uterus during implantation and in response to steroids. Cell Tissue Res. 2012;348(1):199–211.

  60. 60.

    Jeong W, Lim W, Ahn SE, Lim C-H, Lee J-Y, Bae S-M, Kim J, Bazer FW, Song G. Recrudescence mechanisms and gene expression profile of the reproductive tracts from chickens during the molting period. PLoS One. 2013;8(10):e76784.

  61. 61.

    Lee J-Y, Jeong W, Lim W, Kim J, Bazer FW, Han JY, Song G. Chicken pleiotrophin: regulation of tissue specific expression by estrogen in the oviduct and distinct expression pattern in the ovarian carcinomas. PLoS One. 2012;7(4):e34215.

  62. 62.

    Hayashi K, Burghardt RC, Bazer FW, Spencer TE. WNTs in the ovine uterus: potential regulation of periimplantation ovine conceptus development. Endocrinol. 2007;148(7):3496–506.

  63. 63.

    Bar A, Striem S, Mayel-Afshar S, Lawson DEM. Differential regulation of calbindin-D28K mRNA in the intestine and eggshell gland of the laying hen. J Mol Endocrinol. 1990;4(2):93–9.

  64. 64.

    Pines M, Knopov V, Bar A. Involvement of osteopontin in egg shell formation in the laying chicken. Matrix Biol. 1994;14(9):765–71.

  65. 65.

    Hincke MT, Chien Y-C, Gerstenfeld LC, McKee MD. Colloidal-gold immunocytochemical localization of osteopontin in avian eggshell gland and eggshell. J Histochem Cytochem. 2008;56(5):467–76.

  66. 66.

    Pines M, Knopov V, Bar A. Involvement of osteopontin in egg shell formation in the laying chicken. Matrix Biol. 1995;14(9):765–71.

  67. 67.

    Lavelin I, Yarden N, Ben-Bassat S, Bar A, Pines M. Regulation of osteopontin gene expression during egg shell formation in the laying hen by mechanical strain. Matrix Biol. 1998;17(8–9):615–23.

  68. 68.

    Jonchère V, Réhault-Godbert S, Hennequet-Antier C, Cabau C, Sibut V, Cogburn LA, Nys Y, Gautron J. Gene expression profiling to identify eggshell proteins involved in physical defense of the chicken egg. BMC Genomics. 2010;11(1):57.

  69. 69.

    Wolfenson D, Frei YF, Berman A. Responses of the reproductive vascular system during the egg-formation cycle of unanaesthetised laying hens. Br Poult Sci. 1982;23(5):425–31.

  70. 70.

    Inoue K, You-Jun F, Zhuang L, Gopal E, Miyauchi S, Ganapathy V. Functional features and genomic organization of mouse NaCT, a sodium-coupled transporter for tricarboxylic acid cycle intermediates. Biochem J. 2004;378(3):949–57.

  71. 71.

    Qi X, Tan D, Wu C, Tang C, Li T, Han X, Wang J, Liu C, Li R, Wang J. Deterioration of eggshell quality in laying hens experimentally infected with H9N2 avian influenza virus. Vet Res. 2016;47(1):35.

  72. 72.

    Li YP, Bang DD, Handberg KJ, Jorgensen PH, Zhang MF. Evaluation of the suitability of six host genes as internal control in real-time RT-PCR assays in chicken embryo cell cultures infected with infectious bursal disease virus. Vet Microbiol. 2005;110(3):155–65.

  73. 73.

    Bagés S, Estany J, Tor M, Pena RN. Investigating reference genes for quantitative real-time PCR analysis across four chicken tissues. Gene. 2015;561(1):82–7.

Download references

Acknowledgements

We thank Kapil Chousalkar at The University of Adelaide, SA, Australia for his technical support in the design of this experiment. The assistance of Lihong Qin (from Animal and Veterinary Institute, Jilin Academy of Agricultural Sciences, China) in functional annotations of RNA-Seq data is gratefully acknowledged.

Funding

This study was funded by the Australian Egg Corporation Limited, under grant number AECL 1UN121 to Juliet Roberts. The funding body did not involve in the collection, analysis and interpretation of data and in writing the manuscript.

Availability of data and materials

All the data obtained in the current study have been presented in this article. The RNA-Seq sequence raw data-set supporting the results of this study have been deposited at the National Center for Biotechnology Information (NCBI), Sequence Read Archive (SRA) under the Accession Number SAMN10461749.

Author information

SK developed the hypothesis, designed and performed the experiment, analysed and interpreted the data, and drafted the manuscript; JR oversaw the animal trials, administrated the overall research project, assisted with the experiment, analysis and interpretation of the data and critically revised the manuscript; S-BW designed gene expression experiment, analysed and interpreted the data, and drafted the manuscript. All authors reviewed and approved the manuscript for publication.

Correspondence to Shu-Biao Wu.

Ethics declarations

Ethics approval and consent to participate

The experimental setup was approved by the University of New England, Animal Ethics Approval Committee under Authority No. AEC15-118. The protocol was carried out in accordance with the guidelines specified in Australian Code for the Care and Use of Animals for Scientific Purposes 8th edition 2013.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Egg quality variables measured for dividing experimental hens into two different groups. (DOCX 12 kb)

Additional file 2:

Figure S1. Multi-dimensional scaling (MDS) plot showing the expression level of genes in 12 different samples. (PDF 27 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Chicken oviduct
  • Eggshell formation
  • Ions transport
  • Gene regulation
  • Transcriptome profiling
  • Matrix proteins