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

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. Electronic supplementary material The online version of this article (10.1186/s12864-019-5460-4) contains supplementary material, which is available to authorized users.


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 Ca 2+ 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 HCO 3 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].
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.

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 CO 2 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 Ultra-Pure™ 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.

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 Fig. 1 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 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 fea-tureCounts 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 (log 2 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     [49][50][51].

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.

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).
A total of 14,334 gene transcripts were assessed for differential expression after filtering was applied. Differential gene expression analysis showed 691 (log 2 Table 4.

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 (log 2 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.

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 (log 2 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.

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.

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). 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 (R 2 = 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   To retrieve the best homology hit, the sequence IDs were blasted against chicken, duck, turkey and human reference genomes in Ensembl BLAT database. The cut off criterion was established as e-value <10E − 20 . a Represents genes significantly (log 2 fold change > 1.  [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 log 2 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 (log 2 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   [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 (log 2 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 (log 2 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 (log 2 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 (log 2 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 (log 2 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  [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 (log 2 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 log 2 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 For qPCR, the relative expression level of genes at the 5 h and 15 h time-points was calculated in qbase+ software based on 2^-ΔΔCq with genes amplification specific efficiency. For gene expression data normalisation, TBP and YWHAZ were used as reference genes. Plus and minus signs show down-or up-regulation of genes at the 15 h relative to the 5 h time-point