Transcriptome analysis of the uterus of hens laying eggs differing in cuticle deposition

Background Avian eggs have a proteinaceous cuticle. The quantity of cuticle varies and the deposition of a good cuticle in the uterus (Shell-gland) prevents transmission of bacteria to the egg contents. Results To understand cuticle deposition, uterus transcriptomes were compared between hens with i) naturally good and poor cuticle and, ii) where manipulation of the hypothalamo-pituitary-gonadal-oviduct axis produced eggs with or without cuticle. The highest expressed genes encoded eggshell matrix and cuticle proteins, e.g. MEPE (OC-116), BPIFB3 (OVX-36), RARRES1 (OVX-32), WAP (OVX-25), and genes for mitochondrial oxidative phosphorylation, active transport and energy metabolism. Expression of a number of these genes differed between hens laying eggs with or without cuticle. There was also a high expression of clock genes. PER2, CRY2, CRY1, CLOCK and BMAL1 were differentially expressed when cuticle deposition was prevented, and they also changed throughout the egg formation cycle. This suggests an endogenous clock in the uterus may be a component of cuticle deposition control. Cuticle proteins are glycosylated and glycosaminoglycan binding genes had a lower expression when cuticle proteins were deposited on the egg. The immediate early genes, JUN and FOS, were expressed less when the cuticle had not been deposited and changed over the egg formation cycle, suggesting they are important in oviposition and cuticle deposition. The uterus transcriptome of hens with good and poor cuticle deposition did not differ. Conclusions We have gained insights into the factors that can affect the production of the cuticle especially clock genes and immediate early genes. We have demonstrated that these genes change their expression over the period of eggshell formation supporting their importance. The lack of differences in expression between the uterus of hens laying eggs with the best and worse cuticle suggest the genetic basis of the trait may lie outside the oviduct.


Background
The cuticle that covers the eggs of many bird species is invisible [1]. Perhaps for that reason it has been a relatively neglected structure. However, invisibility does not mean it is not important, for birds that live in challenging aquatic environments and where the nests are dirty it is of critical importance for the prevention of bacterial contamination and tends to be much thicker the greater the challenge [1][2][3]. Conversely, species from drier environments seem to have less cuticle or cuticle associated structures [4,5]. The cuticle is a highly glycosylated proteinaceous layer that is specifically deposited on the outside of avian eggs [6], filling the pores in the shell, so preventing aqueous access and bacterial contamination but still allowing gas exchange [7][8][9]. The predominant cuticle proteins in the chicken appear to be the BPI fold containing family B member 3 (BPIFB3 also known as ovocalyxin-36), kunitz-like protease inhibitor, matrix extracellular phosphoglycoprotein (MEPE also known as ovocleidin-116), ovocleidin-17 (OC-17), ovocalyxin 25 (WAP), clusterin (CLU) and retinoic acid receptor responder 1 (RARRES1 also known as ovocalyxin-32) [10]. Lower levels of other proteins have been reported depending on the method used to recover the cuticle [11][12][13]. Despite these proteins having a not dissimilar constitution to the organic matrix of the shell [14], termination of shell formation occurs before the deposition of the cuticle [6]. Therefore, cuticle deposition is not an extension of shell organic matrix production but is a specific event that occurs in the uterus, often referred to as the shell gland pouch, very close to the time of oviposition [6]. This was elucidated by utilizing precise temporal manipulation of the hypothalamo-pituitarygonadal-oviduct axis to manipulate the daily laying cycle and obtain eggs at the same time in the daily cycle which had cuticle deposition or did not have cuticle deposition [6]. However, from this study it became obvious that secretion of the components of the cuticle are not controlled by the proximate endocrine and paracrine factors responsible for the act of oviposition, namely arginine vasotocin (AVT) and prostaglandin [6]. The question of what precise pathways are responsible for the secretion of the cuticle however, remains unresolved.
Previously we speculated that there may even be intrinsic timing events that control the succession of events in the uterus from mineralization through to cuticle deposition and oviposition [6]. We might also expect that factors which are responsible for the overall maintenance of a fully developed oviduct and which act to prevent apoptosis would be important. These factors would affect the overall quality of the cuticle and the expression of genes which control this, these factors might include the receptors for ovarian steroids [15] and growth factors thought to interact with them [16,17].
We also demonstrated on multiple lines of divergent breeds of chicken that the quantity of cuticle is a heritable trait and the hen's level of cuticle deposition on eggs was remarkably constant through its life [10,18]. Most importantly, the quantity of cuticle had a significant effect on bacterial penetration of eggs [8,10,19]. It seems therefore, that the vertical and horizontal transmission of bacteria that threaten the developing embryo [20,21] are prevented by the deposition of a good cuticle. However, it is not known what genetic or environmental factors influence the quantity of cuticle deposited.
Given the state of knowledge about the cuticle and its deposition, we undertook the current study to try and elucidate some of these unanswered questions. Namely, what genes or pathways are involved in the deposition of the cuticle and what genes contribute to the variation in cuticle deposition within a chicken breed? By generating uterus transcriptomes from hens that produced eggs with good or poor cuticle coverage and from hens that were manipulated to produce eggs with or without cuticle, our intention was to increase both our understanding of the uterus transcriptome and to find new avenues to pursue in understanding the biology of this fascinating protective coating found on the surface of eggs.

Results
Zoo-technical data and cuticle measurement in hens from experiments 1 and 2 Experiment 1: comparison of the top and tail of the cuticle deposition distribution At post mortem, there was no difference in body weight between the hens sampled from the top (High) or the bottom (Low) of the distribution for cuticle quality (1968 ± 142 g v 1925 ± 215 g; n = 8; F1,14 0.65) nor the oviduct weight (82.5 ± 6.7 v 81.3 ± 8.8 g; n = 8; F1,14 0.76). As it was the basis of selection of the samples, the cuticle (ΔAbs @640 nm) was more abundant in the High group versus the Low group as indicated by a larger difference in absorbance (High 0.374 ± 0.009 v Low 0.153 ± 0.16 ΔAbs @640 nm; n = 8; F1,14 < 0.001). However, there was less pigment on the High group eggs as indicated by the lower pre-test absorbance measurement (High 0.201 ± 0.011 v Low 0.286 ± 0.016 Abs @ 640 nm; n = 8; F1,14 < 0.001).
Analysis of uterus RNA-sequencing data Experiment 1: comparison of the top and tail of the cuticle deposition distribution Eight biological replicate samples of chicken uterus from each of the extremes of the distribution of cuticle deposition, i.e. High and Low cuticle quality, were analysed. Quality control and data processing steps rejected one replicate from the High cuticle deposition group resulting in 7 High and 8 Low biological replicates of "clean" data totalling 136,267,433 and 185,605,401 aligned reads respectively. Of a total of 18,346 genes present in the chicken genome, our sequencing reads were assembled onto 12,253 genes (RPKM > 1) (Additional file 1: Table  S1). To identify the DEGs between the two groups, we used EdgeR with adjusted P < 0.05 for all comparisons. No DEGs were detected between the samples from the High and Low groups.

Experiment 2: manipulation of cuticle deposition by endocrinological intervention
Eight biological replicate samples of chicken uterus in each of two well-studied experimental manipulations of oviposition timing; i.e. the injection of GNRH1 and AVT were available. Quality control and data processing steps rejected one replicate from the AVT treatment group resulting in 8 GNRH1 and 7 AVT biological replicates of "clean" data totalling 129,731,248 and 197,066,236 aligned reads respectively. Of a total of 18,346 genes present in the chicken genome, our sequencing reads were assembled onto 12,248 genes (RPKM > 1) (Additional file 2: Table S2). To identify the DEGs between the two groups, we again used EdgeR with adjusted P < 0.05 for all comparisons. A total of 3431 DEGs were detected between the treatment groups ( Fig. 1a; Additional file 3: Table S3). Of these, 1736 were significantly up-regulated and 1695 were significantly downregulated in GNRH1 treated hens compared to AVT hens (Fig. 1b), indicating that uterine genes are differentially expressed in these two treatment groups which because one group deposited cuticle whilst the other did not may be related to this process. Genes from each group with a fold-change > 1.5 or < − 1.5 (46 and 115, respectively) were selected for further analysis ( Fig. 1b; Additional file 4: Table S4 and Additional file 5: Table  S5).
Expression profile of genes in the uterus close to oviposition in experiment 1 and 2 In the uterus, the 30 genes with the highest expression are detailed in Table 1. There logCPM values are between 11 and 15.3, where the genes in the data set as a whole had a mean expression of 4.4 and an upper quartile value of 6.1 (Table 1).

GO functional annotation and pathway enrichment analysis of the DEGs in experiment 2
DEGs in experiment 2 were enriched in biological process, cellular component, and molecular function categories by GO analysis (http://www.geneontology.org/). GO terms with P < 0.05 were considered significantly enriched in DEGs. The GO term enrichment analysis showed that there were no enriched categories of GO functional annotation or pathway enrichment for upregulated genes in GNRH1 treated hens compared to AVT hens with a fold change ≥1.5 when subject to analysis by the DAVID Functional Annotation Tool.
The top 10 significant GO functional annotation terms for down-regulated genes in GNRH1 treated hens compared to AVT hens are listed in Table 2. In the Biological process ontology, the highest enrichment of DEGs included signal transduction (26 genes), system development (19 genes) and cell differentiation (19 genes). In the Molecular function ontology, only the glycosaminoglycan binding category (GO:0005539) satisfied the cut-off criteria (P < 0.05) and included 6 genes (HBEGF, CYR61, THBS1, CEMIP, SPP1, REG4). There were no over-represented categories in the cellular component ontology.
Besides the GO analysis, DEGs were also plotted to KEGG reference pathways (https://www.kegg.jp/kegg/ pathway.html) for down-regulated genes in GNRH1 treated hens compared to AVT hens with FDR < 0.05 (Table 3) to collect the molecular interaction, reaction and relation between them. Forty-two KEGG pathways of chicken were assigned. The significantly enriched KEGG pathways (P < 0.05) included pathways such as MAPK signalling (8 genes) and Jak-STAT signalling (5 genes).

Confirmation of gene expression differences by RT-qPCR
Differences in expression between GNRH1 and AVT treatments were confirmed for 13 genes (Fig. 2). These included 10 highly DEGs that were both up and down regulated relative to the GNRH1 treatment. In addition, 3 genes (PER, CRY1 and CRY2) that were not indicated to be in the top rank of DEGs but did differ between treatments, were selected for comparison because they were highly expressed and of biological interest in relation to timing events. PER2 showed differential expression but neither of the CRY genes were differentially expressed in the RT-qPCR (Fig. 2). Although the magnitude of changes in gene expression varied greatly, the results of the 13 genes expression correlated to the RPKM values estimated by RNA sequencing.

Uterus expression of differentially expressed genes at different stages of egg formation
In order to study further the expression of genes that were apparently regulated between the situations were cuticle production was changed, their expression was measured in the uterus at different stages during egg formation including when there was no egg in the uterus (Fig. 3).
FOS expression varied significantly over the egg formation period (F 4,78 = 7.21, P < .001). FOS was high during shell formation; the greatest expression was measured at the middle stage of shell formation; the expression was significantly lower when there was no egg in the uterus and the expression at the late period of shell formation was also significantly lower than either mid or early shell formation (Fig. 3a).
JUN expression varied significantly over the egg formation period (F 4,78 = 5.15, P < .001). The lowest level of expression was measured at the end of shell formation (late) with both the mid and late, pause and white formation expression being significantly lower than the level at the early stage of shell formation (Fig. 3b). HSPB9 had very low levels of expression during all the egg formation stages and there was no significant difference in expression across the egg formation stages or even when no egg was present (pause) (F 4,76 = 0.77; Fig. 3c). SPP1 expression increased over the period of shell formation (F 4,77 = 16.28, P < .001), being significantly higher at the end of it. Expression was very low when no egg was present in the uterus (Fig. 3d). TNFSF10 expression had an opposite trend to SPP1 (F 4,80 = 13.16, P < .001), with significantly higher expression when no egg was in the uterus but no differences between the stages when an egg was present in it (Fig. 3e). PER2 varied over the stages categorised (F 4,79 = 13.96, P < .001) with a significant peak of expression at the end of shell formation, whilst the expression was low during the rest of the stages (Fig. 3f). CRY1 variation was significant (F 4,80 = 2.57, P = 0.045), with a significantly higher expression at the end of shell formation (Fig.  3g). Expression over all the egg formation stages for CRY2 was similar (F 4,78 = 0.72; Fig. 3h). Correlation of expression between genes during stages of egg shell formation Significant correlation existed between the expression of genes across the stages of egg formation, with good correlation between CRY1, CRY2 and PER2 which show correlations in excess 0.6 (P < 0.001; Fig. 4). Other genes showing evidence of correlation around 0.5 (P < 0.001) are TNSF10 with CRY2, FOS, JUN and SPP1; SPP1 with PER2; CRY2 with JUN (Fig. 4).

Discussion
To the best of our knowledge, the present study is the first deep-sequencing analysis of the laying hens' uterus focused on cuticle deposition, which occurs just before oviposition. The differentially expressed genes identified in this study differ considerably from the results of previous studies that compared when an egg was in the uterus and when there was not [22][23][24]. We therefore believe that the experimental regimen applied here captures some of the unique events that are important at the very end of the egg forming process including those that result in cuticle deposition.
We present the analysis of the expression profile of genes in the uterus, close to oviposition. There are 32 individuals sampled at this stage of egg formation from 2 experiments, so we have a highly-replicated data set to identify those genes most highly expressed in the uterus at oviposition. Most studies of the transcriptome focus only on the differences between defined states or tissues and only present that data, however there is much information to be gathered from the study of the most abundant genes in a tissue. When the 30 genes with the   highest expression in the uterus were annotated as might be expected, several of the genes with extremely high expression at oviposition either encode proteins that are known to be expressed and secreted to form the cuticle [10,25] or the organic matrix of the eggshell [14] e.g. MEPE, BPIFB3, RARRES1 and WAP. Other genes are related to mitochondrial oxidative phosphorylation e.g. NADH dehydrogenase subunits, ND4, ND1, ND2, ND5, cytochrome c and b related genes and COX3 that are required for the large amount of energy needed to drive the active secretion of proteins and minerals [26]. These are accompanied by genes important in active transport and energy metabolism, ACTB, GAPDH, ATP6 and ATP1A. A number of these highly-expressed genes were also differentially expressed in the GNRH1 versus AVT experiment described in this study (Experiment 2); MEPE, COX1, COX3, COX2, BPIFB3, Cytochrome b, ATP6, ND5, ATP1A1, ND4, ND2, EIF4A2, UBB, Novel mitochondrial gene, IGLL1, HSPA8, RASD1.
There are however, a number of genes that do not fall immediately into these categories. IGLL1, for example, has previously been reported in the uterus in a microarray experiment [27]. This protein was found in the protein component of the deeper layers of the shell [28] and appears to be a major secreted protein responsible for pigeon fanciers lung [29], but other than the strong indication that it is an epithelial secreted protein and its function remains obscure.
Whey acidic protein (ENSGALG00000020972; WAP) was also very abundant and belongs to the first group of genes encoding proteins secreted as part of the egg shell matrix and in particular the cuticle. It appears to be related to genes found in molluscs and may be involved in the termination of mineral accretion on the shell [30]. It seems likely that this gene produces the protein that has been annotated as similar to Kunitz-like protease inhibitor in proteomic studies of the cuticle [10,25]. It is one of the most abundant proteins in the cuticle [10,25,28] and is synonymous with the protein once named as ovocalyxin 25. Ovocalyxin 25 has previously been attributed with antimicrobial properties [31]. However, the possibility that ovocalyxin 25 is critical in termination of mineral deposition prior to cuticle deposition now seems a very strong hypothesis, given similar molecules can change the formation of mollusc shells [30]. The role of matrix proteins in modulating CaCO 3 crystal formation is well established [32], but both proposed roles are not mutually exclusive.
Overall, the most abundant genes in the uterus reported here are broadly in agreement with the 5 top expressed genes quoted in a recent NGS study that also featured the uterus with the unidentified transcript ENSGALG00000044239 likely being homologous to ENSGALG00000018361, NADH-ubiquinone oxidoreductase chain 5 in this study. The absence of whey acidic protein (ENSGALG00000020972) in that study, which was the 3rd highest in the list reported here, may reflect the time of sampling in this study which was later in shell formation and would reinforce the role of this gene at the end of calcification and in cuticle formation.
Given the focus of the studies on the cuticle, we were also interested in genes that were highly expressed and that might be related to the timely secretion of the cuticle just before oviposition [6]. A notable gene of interest was ENSGALG00000005521, which encodes Period Circadian Regulator 2 (PER2) with an expression of 10.1 log cpm. There were also high levels of other components of the circadian clock mechanism, Cryptochrome Circadian Regulator 2 (CRY2) and Cryptochrome Circadian Regulator 1 (CRY1), at 9.4 log cpm and 7.6 log cpm respectively. We also found CLOCK (3.9 log cpm) and the gene encoding BMAL1, (ARNTL, 5.5 log cpm) at what were average levels of expression for this study. The presence of a functioning clock [33] in peripheral tissues is not new [34] and the elements of the clock have been observed before in the chicken uterus where they were seen to cycle [35]. PER2 was also differentially expressed in a NGS study comparing hen uteri 15 h apart [24]. Although not in the list of genes with the most significant differences in expression, we also found significant differences between the AVT and GNRH1 treated hens in the expression of PER2 (P = 0.0028), CRY2 (ENSGALG00000008436; P = 0.0001), CRY1 (ENSGALG00000012638; P = 0.0002), CLOCK (ENSGALG00000013793; P = 0.001) and ARNTL (ENSG ALG00000005378; P = 0.00002). These differences combined with the high level of expression of PER and CRY genes, confirmed by RT-qPCR, suggests that Clock genes may be important in driving the cyclic expression of genes in the uterus. When examined over the egg laying cycle the expression of PER2 showed strong cyclicity related to the time in the egg formation cycle. Although not as distinct a pattern, CRY1 but not CRY2, expression was also highest at the end of egg formation. Cyclicity of the clock related genes has previously been observed in the oviduct of mice which could be entrained by external stimuli including endocrine signals [36]. In that study clock time appears to have been shifted by signals that alter the timing of the ovulatory surge [36]. Such a phase shift could be critical in the timing of cuticle deposition in the hen's uterus or other key events during shell formation and therefore provides a testable hypothesis.
A further aspect of this study was to compare the uterus transcriptome from the extremes of the population distribution for the trait. A photograph of what the extremes of the cuticle deposition distribution look like is available in the supplemental material of a recent paper along with the associated the genetic parameters of the trait [18]. Hens at the low end of the distribution have virtually no cuticle whilst those at the top have a strong even cuticle deposition resulting in a large numerical difference in cuticle reflectance between the eggs laid by the hens that were used for transcriptome analysis. Interestingly we did not find any detectable differences in gene expression between these two groups. One explanation could be that the genetic effects that result in large phenotypic differences in this trait may not be expressed in the uterus, but may reside in neural structures that control release of the cuticle proteins in the uterus or in a timing mechanism related to ovarian-uterus interaction or other unknown mechanisms. Yet we found no evidence for the AVT prostaglandin system that is thought to mediate the brain to ovary signalling of oviposition timing being involved in cuticle deposition [6].
Alternative hypothesis would be that undetectable differences in expression might be responsible for major Corrected expression on the X-axis was corrected using the geometric mean of LBR and NDUFA1 expression to normalize for any differences between tissues. To better visualize the large differences in expression, the data are presented on two different graphs, the genes in the righthand panel having greater than 10-fold higher expression than those in the left-hand panel changes in cuticle deposition. Given that in the endocrinological intervention using AVT/GNRH we observed clear effects on expression, it seems likely that if there were even small changes in expression we would have detected these. Both experiments had the same statistical power. Given the prominence of clock genes in the Fig. 3 Corrected expression of a FOS, b JUN, c HSPB9, d SPP1, e TNFSF10, f PER2, g CRY1 and h CRY2 mRNA in uterus tissue at different stages of egg formation measured by RT-qPCR (n = 7-11) are represented on the y axis. All genes expression was corrected by the geometric mean of LBR and NDUFA1 expression. The y-axis scale therefore varies depending on the relative level of expression of the genes in the tissue, the larger the number the higher the relative expression of the gene. On the x-axis are the physiological states from which the uterus samples were collected. Early, mid-and late describes the stage of shell formation in the uterus and indicates that the egg was not in the magnum when the samples was taken. Pause represents no egg in the uterus and magnum is when the egg is in that region of the oviduct. The respective F and P values from ANOVA are inset in each graph AVT/GNRH experiment the possibility that variation in their protein expression rather than gene expression is a possibility, indeed in many clock mechanisms there are observations of differences in the profile of the different components of the clock [37]. These differences might be genetic differences in, for example, the half-life of clock proteins that would not be apparent in the transcriptome [37].
In experiment 2 we compared gene expression in birds where we had manipulated the hypothalamo-pituitarygonadal-oviduct axis pharmacologically using either GNRH1 or AVT. Birds on each treatment laid eggs at the same time but either with (GNRH1) or without (AVT) a cuticle. The most differentially expressed genes in this case fell into the GO molecular function glycosaminoglycan binding (0005539): HBEGF, CYR61, THBS1, CEMIP, SPP1, REG4. For those genes tested, all were confirmed to have a lower level of expression in the GNRH1 group who had completed a normal albeit advanced egg formation cycle than the AVT group by RT-qPCR. In experiment 3 where we examined gene expression at different stages of the egg production cycle these same genes expression varied quite significantly. For example SPP1, widely known as osteopontin in relation to mineralisation, was clearly higher at the end of the shell formation period, a fact noted previously [38] and recently in a transcriptome study [22]. These genes must presumably reduce their expression rapidly after oviposition which would make sense as there is no longer an egg in the oviduct and no requirement for events associated with glycosylation of matrix or cuticle proteins.
A further category of genes identified from experiment 2 was a group of genes that had both relatively high expression and a large significant difference in expression, where expression in the GNRH1 group was lower than in the AVT group. This category included the immediate early genes JUN and FOS as well as NR4A3, and HSPB9. There were also several genes associated with the KEGG pathways for both MAPK and Jak-STAT signalling pathways. The immediate early gene products, JUN and FOS, classically act as transcription factors and can form dimers to activate transcription [39]. As the name suggests these genes respond rapidly to stimuli which could explain why the AVT group expression was much higher uterus responded to the AVT injection. This seems plausible because in the natural cycle (experiment 3) expression was lower in the late phase of shell formation. In the chick oviduct JUN and FOS have also been shown to be responsive to the steroid environment, increasing with oestrogen and decreasing in response to progesterone [40][41][42]. In birds the secretion of progesterone peaks at ovulation and is critical to maintaining the ovulatory surge [43]. Previously, it had been proposed that progesterone might be important to cuticle formation [44]. However, in a previous publication we utilised the lack of a progesterone on a pause day to demonstrate there was no direct effect of the presence or absence of a progesterone ovulatory surge on cuticle deposition [6] although arguably the timing effect of the previous surge could still be critical.
Two genes had both high levels of expression and a large significant difference between the AVT and GNRH group. For these genes, TNFSF10 and GKN2 (Ovocalyxin 21), expression was lower in the AVT group which produced eggs without cuticles versus the GNRH group that did have cuticles on the eggs produced. However, it should be noted, that the differences in GKN2 expression was not validated by qPCR. Tumor Necrosis Factor Superfamily Member 10, TNFRSF10 is known as TRAIL and is the ligand for receptors known as TRAILR, but can also be bound by the decoy receptor TNFRSF11B, also known as osteoprotegerin (OPG). OPG is present in the uterus and is highly upregulated as the steroid environment changes with sexual maturity [45]. OPG, under the influence of oestrogen, is thought to prevent the activation of death domain receptors within the oviduct [46] by TRAIL related apoptosis. In the current study, TNFRSF10 had higher expression in the GNRH1 group, which may be related to subsequent fluctuation in ovarian steroids which maintain the oviduct and the extremely rapid regression and apoptosis which occurs when steroid support is removed. This observation was confirmed by a transcriptome study comparing the uterus during and before egg shell calcification where TNFSF10 was higher when no egg was in the uterus [22].
Although relatively few of the top differentially expressed genes identified in this paper are in common with the study of Kahn [22], those that are common to both studies change very dynamically during the egg formation cycle. These include BPIFB3 (ovocalyxin 36), BPIL3, GAL3ST2, GNG4, HABP2, NPTX1, POMC, SLC13A5, SPP1, TNFRSF6B and TNFSF10. Using a microarray approach, Brionne [34] profiled the expression of genes in the uterus in the presence of a hard shelled egg and identified a number of genes some of which are common to both studies; BPIFB3 (ovocalyxin 36), GKN2, BPIL3, POMC, PENK, SPP1 and TNFRSF6B [23] whilst in a similar transcriptome study using RNA sequencing there was shared differential expression observed in the uterus with and without an egg of BPIL3, DUSP4, GAL3ST2, GKN2, MXD4, PENK, RCAN1 and again POMC and SPP1 [24]. The genes in common between our study and these studies are more likely to be involved in shell formation and therefore those genes which are not in common, such as the clock genes and the immediate early genes as prime candidates for the control of cuticle deposition.

Conclusions
By manipulating the production of eggs with or without cuticle endocrinologically we have demonstrated two categories of genes that may be responsible for controlling the production of the cuticle. Clock genes, which determine the timing of cellular events and ensure that the cells respond appropriately to their environment, are possible mediators of cuticle deposition as the timing of deposition is very tightly controlled. Clock gene expression during shell formation show dynamic changes that support the hypothesis. Immediate early genes, which are critical in the activation of cellular processes by external factors, have large fold differences between the physiological states where cuticle has been secreted and those where it has not. The lack of differentially expressed genes when the top and tail of the cuticle deposition were compared suggests that the uterus itself may not be the site of genetic differences in cuticle deposition, therefore other tissues may be responsible for genetic variability of this trait.

Animals and sampling details
All hens used in experiments were sourced commercially as day old layer hen chicks Lohmann Brown (Lohmann GB, Worcester, England) and were housed in our facilities at Roslin. Hens from the top and tail of the distribution were collected from a Lohmann Tierzucht flock with permission. In total 58 animals were used to provide material for the studies in this paper and material from a further 36 animals which had been collected for at the same time as material utilised in previous publications [47].
Hens were reared and maintained as described previously; for experiment 1, the study investigating the Rhode Island Red hens that laid eggs from the top and tail of the distribution i.e. good or poor cuticle, were as described in Bain et al. [10]. This is a commercial pure line that is used in the breeding of Lohmann brown hens and meant we had around 2000 70-week-old hens to select from.
Hens were reared in cages in an environmentally controlled poultry house, receiving 16 h of light per day and were fed following standard commercial practice. Hens were selected based on measurements made at 31 and 50 weeks of age [18]. Genetic correlations across age for the trait of cuticle deposition was near 1 [18] suggesting hens that were laying good or poor cuticle covered eggs earlier in life would continue to do so and indeed that was the case. Nonetheless, two eggs from each hen were measured to ensure the hens were still laying eggs with good or poor cuticle. Twenty-Five High and 17 Low individuals from the distribution were available with tissue and eight of the most extreme individuals were selected from each group for RNA-Seq analysis. Hens were killed immediately after oviposition as specified in Schedule 1 of the UK Animals (Scientific Procedures) Act (1986) using an intra-venous injection of 200 mg/kg of Pentobarbitone Sodium (Euthatal, Merial Animal Health Limited, Harlow, Essex, UK).
For experiment 2, the studies utilizing treatments to create the presence or absence of cuticle were as detailed in Wilson et al. [6]. Briefly the Lohmann Brown (Lohmann GB, Worcester, England) layer hens (Gallus gallus domesticus) were maintained on a 28 h cycle to synchronise their ovulatory cycle.
Uterine tissue was collected from hens that laid eggs with (n = 8) or without (n = 8) cuticle at the same time. This was achieved by causing premature oviposition on a 28-h daily cycle with either GNRH1 peptide where the cuticle was deposited as normal or by using AVT where eggs are produced with no cuticle. This approach and all the peptides used were exactly as used and validated previously to study the physiology of cuticle deposition [6]. Briefly, GNRH1 (28.9 × 10 − 6 mol/kg body weight) was administered intravenously 10 h before a normal oviposition on the 28-h ahemeral light cycle and intravenously administered arginine vasotocin (AVT) (1.05 × 10 − 6 mol/ kg body weight) 4.5 h prematurely. Hens were reared on the floor, to peak-of-lay (22-28 weeks), following commercial management practice, except for the lighting, which was 14 L:10D. Hens were transferred to individual cages prior to the experimental period on a 28-h ahemeral light/dark cycle (14 L: 14D).
Hens were killed immediately after oviposition as specified in Schedule 1 of the UK Animals (Scientific Procedures) Act (1986) using an intra-venous injection of 200 mg/kg of Pentobarbitone Sodium (Euthatal, Merial Animal Health Limited, Harlow, Essex, UK).
For experiment 3 uterus samples were collected at the same time as the magnum samples featured in Gong et al. [47], this allowed uterus samples across the entire spectrum of 24 h egg formation cycle to be obtained and the stage of egg formation was determined by microscopical examination of the shell. Briefly, uterus tissue was collected either when the egg was in the magnum where the albumen is secreted (n = 7), in the uterus (shell gland) and recorded as early (n = 9), mid-(n = 9), and late (n = 9) calcification based on electron microscopy of the shell or during a pause day (n = 11) when there was no ovulation and therefore, no egg in the oviduct.
Selection of sample size was based on previous studies on cuticle deposition [6] to supply data for variance and effect size using for a power calculation (Genstat, 18th edition, VSN International Ltd., Hemel Hempstead, England). The number of animals was the minimum number required to get the a statistically significant effect on cuticle deposition with a power of 1.
Animals in experiment 1 were assigned to experimental treatments by ranking and randomization on body weight. The individual cage was the experimental unit in all experiments.

Zootechnical data, tissue collection and storage
Hens were weighed at cull and the whole oviduct removed and weighed, minus any egg it contained. The uterus pouch was laid out with the pouch facing away from the dissector. A small incision was made in the centre of the pouch and the uterine fluid removed. A longer incision (~5 mm) was then made through the pouch away from the dissector and an approximately square section through the pouch was taken (< 100 mg). The uterine tissues in experiment 1-3 were collected into RNAlater (Life Technologies, Paisley, Scotland), stored overnight at 4°C and subsequently at − 80°C prior to total RNA extraction.

Measurements of egg cuticle deposition and shell pigmentation
Cuticle and pigment were measured spectrophotometrically as described previously [6,10]. Briefly, absorbance of the eggshell at 640 nm was measured before and after staining with a cuticle-specific dye consisting of tartrazine/lissamine green B (Sigma-Aldrich, Poole, Dorset, England). Pre-measurement; absorbance at 640 nm on the unstained egg (Abs @ 640 nm) is reflects the intensity of the pigmentation of the egg, post measurement; difference in absorbance at 640 nm between the unstained egg and the same egg stained with cuticle dye (ΔAbs @640 nm). This represents the amount of cuticle deposition.
RNA isolation, sequencing and quality control analysis RNA was isolated from each tissue sample for sequencing using TRIzol (Invitrogen Ltd., Renfrewshire, Scotland) as described previously. RNA was treated with DNAse I and purified using a RNeasy Mini Kit (Qiagen, Manchester, England). The quality and concentration of RNA was checked with a Nanodrop ND-1000 Spectrophotometer (Thermo Scientific; Waltham, MA USA). 2-4μg of total RNA from each sample were used to make individual bar-coded libraries using Illumina TruSeq RNA Sample Preparation kit (Illumina, Cambridgeshire, England). These libraries were sequenced on an Illumina HiSeq 2000 sequencer (HiSeq high output v4 50PE; Illumina, San Diego, CA, USA). raw data was quality controlled using the FastQC package (Babraham bioinformatics, Cambridgeshire, England). Adapters were trimmed from reads using cutadapt version 1.3 with the parameters -q 30 -m 50 -a AGATCGGAAGAGC.

Transcriptomic data analysis
Trimmed reads were aligned to the Gallus gallus reference genome (Ensembl v82, Galgal 4) with the TopHat2 Aligner (http://ccb.jhu.edu/software/tophat; version 2.0.13) using default parameters [48]. Expression levels of the transcripts were quantified as RPKM (reads per kilobase per million). Read counts were generated using HTSeq-count (version 0.6.0) [49] and Python (version 2.7.3) by Edinburgh Genomics, the University of Edinburgh. Subsequent data filtering and analysis was performed with EdgeR software [50,51]. Genes with RPKM < 1 in at least two samples were not included in the analysis. After normalization, gene-wise generalized linear models (GLM) were fitted on the samples followed by likelihood ratio test to determine differentially expressed genes (DEGs) with a P-value of 0.05 as a threshold. Ensembl gene IDs from each group with a fold-change > 1.5 or < − 1.5 were uploaded to the DAVID Functional Annotation Tool [52,53] and analysed for gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment.

RT-qPCR
mRNA expression differences were verified by qPCR using cDNA from 8 uterine tissue samples from the two treatments (experiment 2). For the tissue expression experiment, cDNA from between 7 and 11 uterus samples was measured for each egg formation stage. Primers were designed using Primer 3 software [54] on exon-exon spans when possible for some of the differentially expressed genes (Additional file 6: Table S6). The qPCR was performed using the Agilent Brilliant III Ultra-Fast SYBR Green qPCR Mix and MX3005P real time system (Agilent; Santa Clara, CA USA). The measurements were carried out in duplicate, essentially as described previously [45,55]. To confirm that only products of the correct length, free from primer dimers, were amplified by each primer pair; the products were assessed by electrophoresis and sequenced directly to confirm identity.
Concentrations were normalized using a geometric mean of the expression of two 'housekeeping genes', LBR and NDUFA1 [56]. One-way ANOVA was used for statistical analysis of data using Minitab Express V.1.3.0. In experiment 3, the Pearson correlation matrix for gene expression was produced using GenStat 18th edition (VSN International Ltd., Hemel Hempstead, England).
Conference proceedings; An brief account of this work has been presented at the UK branch meeting of the Worlds Poultry Science Association [57].