Microarray analysis of genes involved with shell strength in layer shell gland at the early stage of active calcification.

The objective of this study was to get a comprehensive understanding of how genes in chicken shell gland modulate eggshell strength at the early stage of active calcification. Four 32-week old of purebred Xianju hens with consistent high or low shell breakage strength were grouped into two pairs. Using Affymetrix Chicken Array, a whole-transcriptome analysis was performed on hen's shell gland at 9 h post oviposition. Gene ontology enrichment analysis for differentially expressed (DE) transcripts was performed using the web-based GOEAST, and the validation of DE-transcripts was tested by qRT-PCR. 1,195 DE-transcripts, corresponding to 941 unique genes were identified in hens with strong eggshell compared to weak shell hens. According to gene ontology annotations, there are 77 DE-transcripts encoding ion transporters and secreted extracellular matrix proteins, and at least 26 DE-transcripts related to carbohydrate metabolism or post-translation glycosylation modification; furthermore, there are 88 signaling DE-transcripts. GO term enrichment analysis suggests that some DE-transcripts mediate reproductive hormones or neurotransmitters to affect eggshell quality through a complex suite of biophysical processes. These results reveal some candidate genes involved with eggshell strength at the early stage of active calcification which may facilitate our understanding of regulating mechanisms of eggshell quality.

INTRODUCTION

The chicken eggshell is a porous bioceramic container which protects the egg against physical damage and microbial contamination. Avian eggshell consists of the innermost bilayered membranes, a calcified layer composed of a mamillary and pallisade layer, and the outermost cuticle. The calcified layer consists of both inorganic minerals and extracellular matrix. It is well known that the shell mineral amount (thickness or density) is the main factor contributing to the mechanical properties of the eggshell (Ahmed et al., 2005). However, the organic matrix, although its content in the calcified layer is only 2 to 3.5%, is of great importance to the deposition of bicarbonate and calcium ions, and to eggshell strength by controlling calcite crystal nucleation, growth, size and orientation (Greenfield et al., 1984).

The organic matrix in the calcified layer is comprised of a complex suite of components. In the acid soluble part of chicken eggshell matrix, 520 proteins have been identified (Mann et al., 2006), including several abundant proteins such as ovalbumin (Hincke, 1995), ovotransferrin (Gautron et al., 2001b), lysozyme (Hincke et al., 2000), osteopontin (Pines et al., 1995), sialoprotein (Solomon, 1999), clusterin (Mann et al., 2003), ovocleidin-17 (Hincke et al., 1995), ovocleidin-23 (Mann, 1999), ovocleidin-116 (Carrino et al., 1997), ovocalyxin-32 (Gautron et al., 2001a) and ovocalyxin-36 (Gautron et al., 2007). Many of the above components have been reported to undergo various post-translation modifications,which allow them to be effective chelators for interacting with the inorganic materials (Veis, 1989; Reyes-Grajeda et al., 2004; Mann et al., 2007), or to mediate protein–protein interactions to facilitate the assembly of the organic matrix (Lakshminarayanan et al., 2002; Ney et al., 2006).

It has been demonstrated that some genes in hen oviduct are associated with eggshell formation, whose expression is dependent on mechanical strain (Pines et al., 1995; Lavelin et al., 1998; Lavelin et al., 2002). It is proposed that some genes may function as crucial modulators for eggshell quality through regulating signal transduction, ion transportation, expression or modification of organic components, and many other processes. However, despite the importance of eggshell strength in the poultry industry, very few transcriptome-wide studies regarding this trait have been published to date (Yang et al., 2007; Dunn et al., 2009; Jonchère et al., 2010).

It is well documented that various parts of the avian eggshell are formed in specific regions of the oviduct as the egg passes through them. During the laying sequence, about 4 h after previous oviposition, the next egg arrives at and will take about 1h to pass through the white isthmus, in which the bilayered shell membranes are built around the egg. Then the egg enters the initial part of the shell gland, the red isthmus (tubular shell gland), and stays there for about 5 h to form mammillary knobs (Reyes-Grajeda et al., 2004). Finally the egg reaches the uterus (the main part of shell gland) and stays for an additional 15 h to form the palisade layer (Creger et al., 1976). It is known that the mamillary layer is the base of calcite crystal nucleation and crystal growth, and the palisade layer is the main part of the calcified shell, both of which affect global eggshell quality (Reyes-Grajeda et al., 2004; Jonchère et al., 2010).

In this study, we focused on the shell gland (uterus tissue near the red isthmus) at about 9 h post oviposition (corresponding to the early stage of active calcification, or to the transition stage from mammillary knob formation to construction of the palisade layer), and identified differentially expressed genes (DE-genes) in the layers with high shell strength compared to those with weak eggshell. Our results provide insight into the candidate genes involved in the mamillary layer formation and calcification that is crucial to the mechanical properties of avian eggshells.

MATERIALS AND METHODS

Animal treatments

Ninety purebred Xianju hens (a widely-bred Chinese indigenous chicken breed) of 28 weeks old were individually housed in laying cages. Birds were maintained under a cycle of 16 h light and 8 h dark. All birds were fed ad libitum with water and a mash layer diet (165 g protein, 35 g Ca, 11.29 MJ ME/kg, as recommended by NRC of China, 2004).

After 10 d of adaptation for hens, the oviposition time of each egg was initiated to be observed and recorded, then egg weight and shape index (length/width) were measured immediately. Following strength testing, the egg content was discarded and the shell was washed, dried at room temperature and weighed. Shell thickness without membranes was measured with a digital micrometer. Shell index (g/100 cm2) (Sauveur, 1988) was calculated as I = (C/S)×100, in which C is the weight of shell with membranes, S is the shell surface (cm2) with S = 4.68×P2/3 where P = egg weight (g). All above measurements were consecutively carried out daily for 16 d.

Finally, 2 groups of 2 hens with consistent high or low shell breakage strength were found. The differences between the eggshell properties of the selected 4 hens were analyzed by One-way ANOVA variance analysis in SPSS statistic software.

The four hens of interest were humanely sacrificed about 9 h after the previous oviposition. It is of note that all of the sacrificed hens had eggs in their uteruses (Figure 1A). The fat was removed from the uterus tissues near red isthmus and the tissues were then frozen in liquid nitrogen immediately and stored at −80°C. The animal treatments were approved by the Commission for Animal Welfare of Zhejiang A&F University.

Figure 1.

Eggs and forming eggshells obtained from hen shell glands at the moment of tissue sampling. Figure 1A: Eggs obtained from the shell glands of hens when sampling uterus tissues. A-I is the egg from #35 hen, A-II from #19 hen, A-III from #80 hen, A-IV from #40 hen; these hens were all slaughtered at about 9 h after previous oviposition (P.O.). While A-V and A-VI eggs are from another 2 hens culled at 11.5 to 12 h after P.O., respectively. Figure 1B: The forming “shell” sampled at different stage of eggshell formation. B-I is the “shell” from the hen slaughtered at 4 h after P.O.; B-II is the “shell” of above egg from #40 hen. B-III is the “shell” sampled at 12 h post oviposition.

Measurement of eggshell strength

After egg weight and shape index measurements, the uncracked fresh eggs were individually placed lengthways with its blunt end upward in the FHK testing machine (Fujihara Co., Tokyo, Japan), and the vertical pressure was increasingly loaded upon the eggshell until the eggshell cracked and the eggshell strength was recorded as the maximum load (kgf).

RNA preparation

About 500 mg of the tissue of the uterus near red isthmus, including the mucosa, muscularis and outer serosa was powdered under liquid nitrogen. The total RNA was extracted using the RNAiso Plus Mini Kit (TaKaRa, Dalian, P.R. China) according to the manufacturer’s instructions. RNA concentration and purity were measured by a NanoDrop spectrophotometer (NanoDrop® ND-1000, NanoDrop Technologies, DE).

Microarray hybridization and image acquisition

Microarray analysis was performed by the Bioassay Laboratory of CapitalBio Corporation (CapitalBio Co., Beijing, China). Briefly, the RNA integrity was firstly assessed using a Bioanalyzer (Agilent Technologies, Cheshire, UK), then 2 μg of total RNA was used for reverse transcription and biotin-labeled cRNA synthesis according to the manufactures’ instructions, and finally subjected to microarray hybridization. The Affymetrix GeneChip® Chicken Genome Array (Affymetrix, Santa Clara, CA, USA) was used in this study, which contains 38,535 probesets corresponding to >28,000 chicken genes. Following 16 h of hybridization, the arrays were immediately washed, stained and scanned using Affymetrix® GeneChip® scanner 3000 (Affymetrix, Santa Clara, CA, USA), and the image files were processed into raw CEL intensity files using GeneChip Operating Software (GCOS version 1.2).

Pre-processing and normalization of microarray data

The raw intensity files generated by GCOS were imported and processed by R with Bioconductor packages. The total RNA quality was firstly verified statistically again by plotting the 5′-3′ hybridization signal trends across all target transcripts. Then the microarray intensity was processed into transcript expression by the Affymetrix MAS5.0 method implemented in the R package, a procedure including background normalization, PM/MM probe correction, expression summarization and constant normalization on probeset level.

Identification of DE-transcripts

To identify DE-transcripts, the 4 array samples were first grouped into two pairs of high vs. low eggshell strength according to eggshell property differences of the hens (see results). According to Cheuk and Cheng (2011), Affymetrix platform is relatively precise and sensitive in detecting signals, the DE-transcripts were identified as those with fold-change >= 2 in either of the two pairs of comparison and a statistical significant difference between high strength and low strength samples (p<0.05, Welch t-test). It is of note that the log-odds values (Lods) of expression fold-change were used in the analysis; therefore, the DE-transcripts always have an absolute Lods value no less than 1 (|Lods|≥1).

Gene ontology enrichment analysis

Gene ontology enrichment analysis for DE-transcripts was performed using the web-based GOEAST (Zheng and Wang, 2008) Affymetrix analysis tool, with FDR cut-off of 0.05 using Yekutieli’s FDR adjustment method.

Validation of differential expression by qRT-PCR experiments

Twenty-one DE-transcripts, with a fold-change ranging from low to high, were selected for further validation with qRT-PCR experiments; and all two groups of microarray samples were tested. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the internal reference in all the PCR experiments. The primer sequences for qRT-PCR experiments can be found in Table 1.

Table 1.

Descriptions of specific primers used for real-time RT-PCR

Gene symbol	Accession no.	Forward primer (5′-3′)	Reverse primer (5′-3′)	Amplicon (bp)
ACYP2	XM_419292	CGGCTCGCTCAAGTCGGTGG	GGCCCTGAACTTGGCCCGTC	152
AMDHD1	XM_416158	GCACTGGGAAGTGCGTATTGCCA	TCTTCCGTGGCCTTCCTGGTGT	175
ATP6V1A	NM_204974	TGCAACATGGCAGGTGCTGCT	TGCCAGGCCCCAGTTCCACT	187
CA5B	XM_414195	CAGCTTGGCCACCTGCACTCC	ACACGTCGCTGGGTCGTAGCT	175
CHST3	NM_205121	TGATGGCCACCACACGCACC	CTGCAGCACGTCGCGGTACA	170
COL12A1	NM_205021	AGGCGAGTCTTCCCCGACGG	GCGCTGTCCTCATGTCTGCCC	171
CRABP1	NM_001030539	CGCCCCGCCATGCCTAACTT	AACTGGTCCCCGTCCTGGCG	161
CRYBB1	NM_204180	ACCTGGCGGACTGCGGGTT	CGGTAGCTGCTGGACCAGGTG	151
EXOC6B	XM_420892	AACCCCACCACAGCCCTCGT	TGGCTGTTGATGAGGCCGCG	149
FGB	XM_420369.2	GCTGCTCCTGCTGCTCCTGC	GTGCCACGGGCCTGAGTGTG	155
GAS2L3	XM_416172	GGAGTAGTGCTGGCAGTCCTGC	CCTGGGCCGTGTCTGGGAGT	193
GIT2	NM_204206	TCGCTTGCCATGCCGTGAGG	GCAACGTGGAGCGGGGTGTT	168
MAN1A2	XM_416490	ACGTGGACACCAGCAAGGGGG	TCCTTTGCCTCTTCCAGGGCCTTT	148
NDST4	XM_420638	CGAGCAGCTTCCCTCATCCCCAA	TGCCCAGGGGCTTGACGTAA	156
NPY	NM_204587	GAGGACGCTCCCGCAGAGGA	TCGAAGGGTCTTCAAACCGGGA	175
OC416916	XM_415207	TGGAGGTGGAGCACAAACATCTGC CCACCGAGCACACAGCCAGAAA		200
PLCXD1	NM_001128637	CCTGGCCTGCAGGAATTTTGATGG	AGCCACGCTGCCACATGGTC	137
RCJMB04_34k20	NM_001031112	GGACAGGCGGGCGAGAGAGT	TGGTGGTAACACGCACGCTGA	126
SLC8A1	XM_415002	CGTGTTTGTGGCACTGGGGACA	ATGGCCGCGATGGACCAAGC	159
TBXAS1	XM_416334	TGTGTGGTGCTGGGACAGCGT	ATACAGCCACGGGGTCCTGCT	188
WDR72	XM_425069	GGCTGTTATCAGGGGGCCAGGA	GCACACGCAGCACACTACGC	161
GAPDH	NM_204305	GGGCTGCTAAGGCTGTGGGG	TCAGGGGCCCATCAGCAGCA	177

To begin, total RNA was individually reverse transcribed with the SYBR® PrimeScript™ RT-PCR kit II (TaKaRa, Dalian, China) according to the manufacturer’s instructions. Then above RT-PCR kit was further used for fluorescence detection on an ABI Prism 7500 Sequence Detection System (Applied Biosystems, USA). All samples were analyzed in triplicates.

Dissociation curve analysis was conducted to ensure that a single PCR product with appropriate size was amplified in each reaction. On the other hand, the examination of PCR efficiency was performed based on LinRegPCR program (12.X) (Ramakers et al., 2003; Ruijter et al., 2009) to ensure internal and target transcript primers were amplified with similar efficiency.

The differential expression levels (Log2 units) were calculated using the equation Log2 units (high versus low) = −ΔΔCt, where ΔΔCt = (Ctth-Ctih)-(Cttl-Ctil). Ct is the threshold cycle number when the amount of amplified product reaches a stable threshold. Ctth and Ctih represented the Ct of target transcript and internal reference transcript of “high eggshell strength sample”, respectively. Correspondingly, Cttl and Ctil represented the Ct of target transcript and inner-reference transcript of “low eggshell strength sample”, respectively.

RESULTS

Eggshell quality of hens under study

Among 90 tested hens, only 6 of them laid eggs at a similar laying rate with consistent high eggshell strength (defined as ≥4.5 kgf) or low strength (defined as ≤3.5 kgf). Two of these 6 hens were sacrificed at about 11.5 to 12 h post oviposition, but the eggshells collected from the shell glands showed more calcification extent than expected (Figure 1A and 1B). To focus on the initial stage of active calcification, we decided to use the uterus tissues near red isthmus from the remaining 4 hens, namely #19, #35, #40 and #80, at about 9 hr post oviposition. The eggs and eggshells harvested from these 4 hens at the moment of tissue sampling are shown in Figure 1A and Figure 1B. The eggshell breaking strength is also shown in Table 2 for these 4 hens, with consistent high (#19 and #40) or low (#35 and #80) eggshell quality.

Table 2.

Parameters related to eggshell quality of hens in this study

Hen	shell strength (kgf)	shell thickness (mm)	shell index (g/100 cm2)	shell weight (g)	egg weight (g)	shape index
#40	5.17±0.40 A	0.367±0.016 A	8.07±0.28 A	4.297±0.144 a	38.46±1.52 C	1.274±0.024 B
#19	4.75±0.21 A	0.328±0.012 B	7.18±0.25 B	3.968±0.195 b	40.72±1.26 B	1.322±0.030 A
#35	2.99±0.71 B	0.324±0.023 B	6.97±0.68 B	4.363±0.426 a	48.60±1.38 A	1.278±0.032 AB
#80	2.54±0.69 B	0.272±0.022 C	6.25±0.71 C	3.451±0.417 c	40.58±1.06 B	1.313±0.052 AB

Values are from eggs laid by each hen of interest during the period of observation. Distinct capital letters in the same column indicate parameters between hens with a significant difference (p< 0.01), and distinct small letters indicate the significant difference is at level p<0.05.

To eggs from #19 and #35 hens, the differences of both shell strength and shell weight were very significant (p<0.01, t-test), but there was no significant difference (p>0.05) for other eggshell quality metrics, such as shell thickness or shell index (Table 2). On the other hand, the differences of all of above eggshell metrics between eggs from #40 and #80 hens were very significant (p<0.01) (Table 2). To get rigorous microarray data, therefore, we grouped the #19 with #35 hens due to the similarity of some of eggshell mechanical properties of the paired individuals; while #40 and #80 hens were also grouped as another pair.

Differentially expressed transcripts

The expression level of all probesets in 4 array samples were analyzed, and 1,195 DE-transcripts between uterus samples with high shell strength and low shell strength were identified. These DE-transcripts correspond to 941 unique genes. Among them, 407 genes were up-regulated in high strength samples comparing to low strength samples, and the other 534 genes were down-regulated. The expression profile of all 1,195 DE-transcripts is shown in the heatmap in Figure 2. As shown in the heatmap, samples #19 with #40 and #35 with #80 were grouped as clusters among different samples, consistent with the similarity of their eggshell quality.

Figure 2.

Heatmap and dendrogram of differentially expressed transcripts (DE-transcripts). Each cell represents the (normalized) gene expression value for given DE-transcript (row) in the specified sample (column). Cell colors indicate gene expression level: red: highly expressed; yellow: medium expression; green: lowly expressed. Row-side and column-side dendrogram represent the hierarchical clustering of DE-transcript expression for different transcripts or samples, respectively. Clustering is based on “complete-linkage” method using Euclidean-distance.

According to gene ontology annotations, the DE-transcripts are involved in a variety of biological processes. The most prominent DE-transcripts were found related to the following processes: signal transduction (88 DE-transcripts), ion transport and extracellular matrix organization (77 DE-transcripts), carbohydrate metabolism and protein modification (26 DE-transcripts) (Table 3).

Furthermore, avian calcified eggshell is a biomaterial composed of calcium salt and special ECM. The ECM is mainly comprised of collagens, glycoproteins and proteoglycans. Among the DE-transcripts, COL8A2, COL12A1, COL13A1, LOC424798, LAMA2, LAMA4, LAMB4, and LAMC1 may be related to extracellular matrix formation; while CHST3, GALNTL1, NDST4, LARGE, POFUT2, RCJMB04_28l23, and MAN1A2 are all localized in the endoplasmic reticulum or Golgi apparatus, and likely mediate the processes of carbohydrate metabolism, or posttranslation glycosylation modification.

Gene ontology (GO) term enrichment of DE-transcripts

It is of note that although many DE-transcripts were found related to various biological processes according to their ontology annotations, they are not necessarily correlated to the eggshell quality, due to random noise or other non-specific confounding factors commonly existing in microarray or other high-throughput experiments. Therefore, using web-based GOEAST (Zheng and Wang, 2008) we further identified significantly enriched GO terms among all the DE-transcripts. According to biology processes or molecular functions, the enriched GO terms can be roughly classified into several groups (Tables 4 and 5).

Table 4.

Enriched gene ontology (GO) terms revealed from identified DE-transcripts according to biological_process ontology

Group	GOID	Term	p	Gene symbol or representative public ID
Reproductive hormone synthesis and regulation	GO:0060126	Somatotropin secreting cell differentiation	0.013	OTX2, WNT4
	GO:0021984	Adenohypophysis development	0.044	OTX2, WNT4
	GO:0032355	Response to estradiol stimulus	0.044	SOCS2, AREGB
Signal transduction	GO:0030514	Negative regulation of BMP signaling pathway	0.030	TOB1, GREM1
Biophysical processes	GO:0048741	Skeletal muscle fiber development	0.012	SLC23A2, CHAT
	GO:0015074	DNA integration	0.012	LOC770294, LOC770705, ENS-3
	GO:0055117	Regulation of cardiac muscle contraction	0.013	P2RX4, NKX2–5
	GO:0009409	Response to cold	0.018	IL4, SLC27A1
	GO:0048747	Muscle fiber development	0.022	SLC23A2, CHAT
	GO:0046209	Nitric oxide metabolic process	0.024	P2RX4, CPS1
	GO:0007586	Digestion	0.030	PGA5, PRSS2, LOC396365
	GO:0015849	Organic acid transport	0.033	SLC23A2, OCA2, LOC770309, SLC7A14, SLC27A1
	GO:0046942	Carboxylic acid transport	0.033	SLC23A2, OCA2, LOC770309, SLC7A14, SLC27A1
	GO:0055002	Striated muscle cell development	0.034	SLC23A2, CHAT, TTN, NKX2–5
	GO:0006942	Regulation of striated muscle contraction	0.037	P2RX4, NKX2–5
	GO:0002028	Regulation of sodium ion transport	0.044	NKX2–5, NEDD4L
Reproductive biophysical processes	GO:0060748	Tertiary branching involved in mammary gland duct morphogenesis	0.009	WNT4, AR
	GO:0060745	Mammary gland branching involved in pregnancy	0.013	WNT4, AR
	GO:0060562	Epithelial tube morphogenesis	0.019	DEAF1, WNT3, GREM1, WNT4, NKX2–5, HOXA11, AR, AREGB
	GO:0060444	Branching involved in mammary gland duct morphogenesis	0.020	WNT4, AR, AREGB
	GO:0009994	Oocyte differentiation	0.024	WNT4, GDF9
	GO:0048599	Oocyte development	0.024	WNT4, GDF9
	GO:0060603	Mammary gland duct morphogenesis	0.033	WNT4, AR, AREGB
	GO:0060135	Maternal process involved in female pregnancy	0.037	WNT4, AR

GOID represents the identifiers, and Term represents term definitions for Gene Ontology term entities. p: p-value of significance (Welch t-test).

Table 5.

Enriched gene ontology (GO) terms revealed from identified DE-transcripts according to molecular function ontology

Group	GOID	Term	p	Gene symbol or representative public ID
Signal transduction	GO:0005030	Neurotrophin receptor activity	0.013	NTRK1, NTRK2
	GO:0001614	Purinergic nucleotide receptor activity	0.017	P2RX4, P2RY2, ENSGALG00000012080
	GO:0016502	Nucleotide receptor activity	0.017	P2RX4, P2RY2, ENSGALG00000012080
	GO:0043121	Neurotrophin binding	0.024	NTRK1, NTRK2
	GO:0035586	Purinergic receptor activity	0.026	P2RX4, P2RY2, ENSGALG00000012080
	GO:0004888	Transmembrane signaling receptor activity	0.049	OXTR, LOC431251, SSTR3, CHRM2, ADRA2B, P2RX4, P2RY2, EDNRB2, GABRB2, GABRG2, LOC428961, NPFFR2, GRIN2B, GRIN3A, NTRK1, NTRK2, EPHB6, DDR2, TMPRSS6, PCSK5, CCKAR, IFNAR2, CSF1R, TLR5, OR10A7, LOC768958, LOC769317, LOC777484, GPR39, GPR97, ENSGALG00000017405, ENSGALG00000017093, ENSGALG00000012080.
Biophysical processes	GO:0003951	NAD+ kinase activity	0.013	C5orf33, NADK
	GO:0005319	Lipid transporter activity	0.049	ATP11C, ATP8A2, ATP8B3, APOB, LOC769564, SLC27A1

GOID represents the identifiers, and Term represents term definitions for gene ontology term entities. p: p-value of significance (Welch t-test).

A group of processes are involved in reproductive hormone regulation, which contain Somatotropin secreting cell differentiation (GO:0060126), adenohypophysis development (GO:0021984), and response to estradiol stimulus (GO:0032355) (Table 4).

As shown in Table 4 and Table 5, many DE-transcripts are involved in signal transduction, such as GO terms purinergic nucleotide receptor activity (GO:0001614), nucleotide receptor activity (GO:0016502), purinergic receptor activity (GO:0035586), transmembrane signaling receptor activity (GO:0004888), and negative regulation of BMP signaling pathway (GO:0030514). Among them, GO:0004888 dominantly contains 33 transcripts encoding signal receptors, and these receptors could be further classified into several subgroups: OXTR, LOC431251 and SSTR3 belong to reproductive hormone receptors; CHRM2, ADRA2B, P2RX4, P2RY2, EDNRB2, GABRB2, GABRG2, LOC428961 and NPFFR2 function as receptors mediating neurotransmitters or neuropeptide; GRIN2B and GRIN3A could modulate the efficiency of synaptic transmission; NTRK1 and NTRK2 belong to the receptor tyrosine kinase (RTK) family, and are involved with neurotrophin (GO:0005030 - neurotrophin receptor activity; and GO:0043121 - neurotrophin binding) (Table 5).

Besides various enriched molecular function shown above, many biophysical processes are also found to be enriched among the DE-transcripts, including a series of processes and subgroups (Tables 4 and 5). GO:0003951 (NAD+ kinase activity) modulate the metabolism or redox in cell (Table 5). Enrichment of GO:0009409 (response to cold) may reflect the fact the rearing condition of experimental hens was in the winter at room temperature about 2 to 10°C. GO:0046209 (nitric oxide metabolic process) may regulate vascular or smooth muscle relaxation or other functions. GO:0002028 is involved in ion transportation. While the subgroup processes of muscular development and activity include skeletal muscle fiber development (enrichments of GO:0048741, GO:0048747 and GO:0055002) and striated muscle contraction regulation (enrichments of GO:0055117 and GO:0006942). It is of note that there is almost no striated muscle in avian uterus except smooth muscle. However, the chicken genome project was completed in 2004, and the functional gene database of G. gallus remains incomplete, some ontology annotations of DE-genes may refer to mammalian homologs, which may account for our results. The genes related to muscular cell contraction are likely to modulate the mobility of uterus to facilitate egg rotation and calcification (Johnson, 1986; Jonchère et al., 2010). Similarly, there is no digestion process in the uterus, three genes in GO:0007586 (digestion process), PGA5 (an aspartic acid protease, which is involved in ovulation (Peluffo et al., 2011)), PRSS2 and LOC396365 (preprogastrin), are likely to promote the maturation of secretary extracellular proteins or regulate the secretion of uterus glands and mobility of uterus.

The final group of reproductive biophysical processes also includes several subgroups of processes (Table 4). Epithelial tube morphogenesis (GO:0060562) may regulate the development of uterus glands (tubular epithelial glands). Oocyte development subgroup contains oocyte differentiation (GO:0009994) and oocyte development (GO:0048599). Female pregnancy subgroup contains enrichments of GO:0060135, GO:0060745, GO:0060748, GO:0060444 and GO:0060603.

Overall, laying is an avian reproductive behavior, and eggshell calcification is regulated by relative reproductive hormones and neurotransmitters, which may finally affect eggshell quality through a complex suite of biophysical reactions.

Confirmation of DE-transcripts by qRT-PCR

21 DE-transcripts (9 up-regulated and 12 down-regulated) were chosen for validation using qRT-RCR experiments, and the four microarray samples were tested in pairs for #19 vs #35 and #40 vs #80, respectively. As shown in Figure 3, 16 out of the 21 tested transcripts (76%) were confirmed by qRT-PCR experiments, though the absolute fold-change values are slightly different. The remaining transcripts, CRYBB1, EXOC6B, LOC416916, MAN1A2 and CHST3, showed inconsistent differential expression between qRT-PCR and microarray experiments.

Figure 3.

Real-time RT-PCR validation of microarray data. Expression levels of the first 12 transcripts (AMDHD1, ATP6V1A, CRYBB1 FGB, GAS2L3, GIT2, MAN1A2, NDST4, PLCXD1, RCJMB04-34k20, SLC8A1 and TBXAS1) were down-regulated in microarray experiment, and the last 9 transcripts (ACYP2, CA5B, CHST3, COL12A1, CRABP1, EXOC6B, LOC416916, NPY and WDR72) were up-regulated in microarray experiment. Among the above transcripts, five transcripts (CRYBB1, MAN1A2, CHST3, EXOC6B and LOC416916) failed to be verified by real-time RT-PCR.

GAPDH, CHST3, GALNTL1, NDST4, LARGE, SP1, RHOBTB2, and WDR72 were selected to examine the PCR efficiency. The results showed the PCR efficiency of these genes ranged from 86.8% to 94.2%, and the PCR efficiency of inner reference (GAPDH) and other genes seemed nearly similar.

DISCUSSION

Laying is regarded as avian reproductive behavior, which is regulated by reproductive hormones and neurotransmitters. The chicken oviduct has been extensively used as a model to study hormonal induction of protein synthesis (Khuong and Jeong. 2011). Under the control of steroid hormones or neurotransmitters, the tubular gland epithelial cells synthesize and secrete a great variety of proteins to form egg white and eggshell when egg passes through the oviduct (Mann et al., 2006). In this paper, 1,195 DE-transcripts have been identified to be related with eggshell strength. GOEAST analysis further identify some significantly enriched GO terms, and the enriched GO terms suggest that some DE-transcripts mediate reproductive hormones or neurotransmitters to affect eggshell quality (Tables 4 and 5).

Both terms GO:0060126 and GO:0021984 are involved in reproductive hormone regulation, and share two genes, OTX2 and WNT4.

Otx2 is a paired-like homeodomain transcription factor, which can mediate GnRH (gonadotropin releasing hormone) signaling (Kelley et al., 2000). Functional studies revealed that Otx2 is required as early as gastrulation for neural induction, and even for brain development (Rhinn et al., 1998). However, Otx2 is also of importance for neurogenesis and cellular proliferation in multiple other tissues (Layman et al., 2011).

As a member of the WNT family, Wnt4 is a secreted glycoprotein signaling molecule and involved in paracrine signaling (Diaz et al., 2011). Wnt4 is critical for female sex determination and differentiation (Chen et al., 2011). In the female, Wnt4 is positively involved in ovarian development; while in the male mutated WNT4 will result in aberrant testis development (Diaz et al., 2011; Barrionuevo et al., 2012). On the other hand, Wnt4 is also potent to regulate the development of the female reproductive tract (Franco et al., 2011). Furthermore, WNT4 is expressed postnatally in ovarian follicles and corpora lutea, and its expression increases in response to gonadotropin (Hsieh et al., 2002). Wnt4 mediates follicle development and fertility by regulating the expression of genes involved in steroidogenesis, prostaglandin biosynthesis, tissue remodeling, and angiogenesis (Hsieh et al., 2002; Boyer et al., 2010). Moreover, excluding its reproductive contributions, WNT4 is also tightly associated with bone strength (Zmuda et al., 2011).

Our results also show that some DE-transcripts are involved in signal transduction (Tables 4 and 5), among which, NTRK1, NTRK2, P2RX4, and P2RY2 are overlapped in multiple enriched GO terms (Table 5).

Ntrk1, also named TrkA, and Ntrk2 TrkB, are two members of the neurotrophic tyrosine kinase receptor (NTKR) family. These kinases are membrane-bound receptors mediating various functions of neurotrophins, such as cell survival, migration, outgrowth of axons and dendrites, synaptogenesis, remodeling of synapses, and synaptic transmission (Ohira1 and Hayashi, 2009). So far, several neurotrophins have been well studied, such as nerve growth factor (NGF), brain-derived neurotrophic factor (BDNF), neurotrophin 3 (NT-3), and NT-4/5. NTKRs are high affinity receptors of neurotrophins. TrkA mediates the biological response of NGF, while BDNF and NT-4/5 are the preferred ligands for TrkB (Ohira1 and Hayashi, 2009). Additionally, NTKRs also play roles in some biomaterials. NT-4 may modulate proliferation and differentiation of the dental epithelium and promote production of the enamel matrix via the TrkB-MAPK pathway (Yoshizaki et al., 2008).

Both P2RX4 and P2RY2 are purinergic receptors. Purinergic receptors are subdivided into metabotropic P2Y receptors and ionotropic P2X receptors. P2Y receptors are coupled to G-protein and trigger inositol 1,4,5-triphosphate (IP3)-induced intracellular Ca2+ release following activation of phospholipase C, while P2X receptors are ligand-gated ion channels. P2RX4 will be discussed later, while P2RY2 performs a dominant role in calcium signaling during osteoblast differentiation (Nishii et al., 2009). It is known that extracellular ATP, UTP, and PPi can strongly block the mineralization of matrix nodules, while this potent inhibition of bone formation is mediated by P2RY2 (Orriss et al., 2007). Furthermore, P2RY2 is also involved in inhibition of intercellular communication between osteoblasts (Hoebertz et al., 2003).

At present, there are at least three cDNA microarray studies globally investigating the gene expression in chicken shell gland (Yang et al., 2007; Dunn et al., 2009; Jonchère et al., 2010), but the overlap among the DE-genes from these studies is not plentiful. Different animals, tissue samples or treatment methods may partially account for this problem.

Yang et al. (2007) harvested uterus tissues at 2 h post oviposition, and screened out 34 known genes in the shell glands between hens with low and high egg productive rates. This study and our data share a single gene CALD1 (caldesmon 1) (Figure 4). CALD1 is a ubiquitous actin and calmodulin binding protein, and functions as a substrate for mitogen-activated protein kinase (Childs et al., 1992) or as serine and threonine kinases (Sutherland et al., 1994).

Figure 4.

MA-plot of all microarray tested chicken-genes. X-axis: the average normalized expression values across all 4 (strong+weak) eggshell samples (in Log2 scale); Y-axis: the log-odds ratio between the average expression values of strong vs. weak eggshell samples; grey dots: nondifferentially expressed genes (non-DE genes); red dots: differentially expressed genes (DE-genes) identified in this study; big blue dots: DE-genes (CALD1, NADK, LOC422993, LAMP3, P2RX4, FSTL1, TUBGCP4, WDR77, RCJMB04_6g16, PWP1 and SGK1) reported in previous studies.

Dunn et al. (2009) identified 266 DE-genes in shell glands from 25-week old mature hens comparing to 12-week old juveniles from high and low bone quality lines, respectively. The tissues sampled when eggs passed through the oviducts but not in shell glands. Three DEgenes are also found in our data: NADK (NAD kinase), LOC422993 (Similar to interferon-induced membrane protein Leu-13/9–27), and LAMP3 (lysosomal-associated membrane protein 3) (Figure 4), suggesting potential crucial function of these genes in not only early stage of eggshell calcification but also other stages of eggshell formation.

Jonchère et al. (2010) used the 40-week old hens at 18 h post oviposition (corresponding to the rapid phase of calcification), and identified 469 DE-known genes in uterus versus both white isthmus, and magnum. There are 7 genes consistently identified in their study and our data, such as P2RX4 (purinergic receptor P2X, ligand-gated ion channel, 4), FSTL1 (follistatin-like 1), TUBGCP4 (Tubulin, gamma complex associated protein 4), WDR77 (WD repeat domain 77), RCJMB04_6g16 (microtubule-associated protein 1 light chain 3 beta), PWP1 (PWP1 homolog in S. cerevisiae) and SGK1 (serum/glucocorticoid regulated kinase 1) (Figure 4).

On the other hand, three additional DE-genes in our data were previously found in the acid soluble part of chicken eggshell organic matrix (Mann et al., 2006). These three genes, FSTL1 (follistatin-like 1), CAMK2D (calcium/calmodulin-dependent protein kinase (CaM kinase) II delta) and KRT75 (keratin 75) (Figure 4), could reflect potential interaction of eggshell calcification and organic matrix formation.

Among these overlapping DE-genes, both P2RX4 (Jonchère et al., 2010) and NADK (Yang et al., 2007) are also present in our enriched GO terms (Tables 4 and 5), and FSTL1 (Mann et al., 2006; Jonchère et al., 2010) occurs in more than three relative studies.

P2RX4 is one member of the P2X receptors (P2RX). P2RX are ionotropic ATP-gated ion channels conducting Ca2+ inflow (Fodor et al., 2009), with high capability of Ca2+ permeabilities corresponding to at least 100-fold those of Na+ (Burnashev, 1998). In chondrogenic mesenchymal cells, P2X4 receptors could conduct Ca2+ inflow to elevate intracellular Ca2+ levels, and finally promoting extracellular matrix production (Fodor et al., 2009). Eggshell calcification requires considerable ion transportation, especially Ca2+, and various extracellular matrix synthesis and secretion, whether and how P2RX4 channels regulate these processes requires further studies.

NAD kinases (NADKs) are a family of enzymes transferring a phosphate group from ATP to NAD to generate and maintain the cellular NADP pool (Pollak et al., 2007). It is reported, that during development of placenta, the expression level of NADK appears drastically elevated (Lerner et al., 2001).

Fstl1 is a secreted glycoprotein belonging to the BM-40/SPARC/osteonectin family containing both calcium-binding domain and Follistatin-like domain (Hambrock et al., 2004). As a mesenchymal factor, Fstl1 is critical for oviduct development, and determines the differentiation of secretary epithelial cells and ciliated epithelial cells in the oviduct (Umezu et al., 2010). This means Fstl1 may modulate chicken endometrium development during eggshell formation. However, Fstl1 is also present in the organic part of eggshells (Mann et al., 2006), and Jonchère et al. (2010) propose it may be a uterine antiprotease

IMPLICATIONS

Above all, using Affymetrix Chicken Array, 1,195 DE-transcripts were identified in the shell gland between “high shell strength” and “low shell strength” hens, which represent 941 unique known genes. According to gene ontology annotations, these transcripts are involved in a wide range of biological processes; the most prominent DE-transcripts relate to signal transduction, metabolism, extracellular matrix, or ion transport and homeostasis, and so on. Furthermore, Gene Ontology (GO) term enrichment of DE-transcripts suggests that avian eggshell calcification is likely to be regulated by relative reproductive hormones and neurotransmitters, which may finally affect eggshell quality through a complex suite of biophysical processes.

Table 3. i)

DE-transcripts related with signaling, ion transportation, extracellular matrix protein, and carbohydrate metabolism or post-translation glycosylation modification

Gene symbol	Transcript ID	Log2 units (strong VS weak)	p-value	Category
LOC771699	XM_001234946	3.853	0.0004	signaling
SH3PXD2A	XM_421741	3.410	0.0397	signaling
LOC429955	XM_427511	3.296	0.0034	signaling
PDCL2	XM_420702	3.100	0.0115	signaling
LOC430487	XM_428042	3.055	0.0074	signaling
SEMA3G	XM_414289	2.933	0.0213	signaling
RHOBTB2	XM_001232709	2.905	0.0069	signaling
RXFP1	XM_420385	2.795	0.0205	signaling
PIK3C2B	XM_417956	2.777	0.0194	signaling
OR10A7	XM_425093	2.768	0.0067	signaling
NPY	NM_205473	2.401	0.0329	signaling
PDE8B	XM_425218	2.401	0.0498	signaling
GREM2	XM_419552	2.284	0.0423	signaling

Table 3. ii)

DE-transcripts related with signaling, ion transportation, extracellular matrix protein, and carbohydrate metabolism or post-translation glycosylation modification

Gene symbol	Transcript ID	Log2 units (strong VS weak)	p-value	Category
SPAG9	XM_420098	2.215	0.0423	signaling
LOC396365	NM_205400	2.166	0.0009	signaling
MAPKBP1	XR_026772	2.157	0.0187	signaling
OXTR	NM_001031569	2.128	0.0134	signaling
STC2	XM_414534	2.120	0.0135	signaling
MPP3	XM_418108	1.921	0.0491	signaling
GREM1	NM_204978	1.903	0.0220	signaling
C20orf32	XM_417499	1.889	0.0117	signaling
CRABP1	NM_001030539	1.778	0.0115	signaling
C1orf107	NM_001031051	1.770	0.0411	signaling
RGS9	XM_415685	1.657	0.0177	signaling
ARHGEF12	XM_417890	1.437	0.0014	signaling
LTBP3	XM_426444	1.433	0.0409	signaling
NGEF	NM_001010841	1.390	0.0307	signaling
PDE1A	XM_421969	1.363	0.0127	signaling
CRHBP	XM_424801	1.355	0.0159	signaling
SRGAP1	NM_001080101	1.335	0.0426	signaling
PDE9A	XM_416748	1.196	0.0475	signaling
NPFFR2	NM_001034825	1.124	0.0283	signaling
WNT4	NM_204783	1.114	0.0180	signaling
FGD4	XM_416365	1.099	0.0121	signaling
SOCS2	NM_204540	1.033	0.0377	signaling
TOB1	NM_001001467	1.015	0.0347	signaling
RND3	XM_422158	0.887	0.0094	signaling
ARHGAP28	XM_419140	−4.485	0.0231	signaling
WNT3	NM_204675	−3.902	0.0422	signaling
LOC428961	NM_001142671	−3.469	0.0058	signaling
VAV2	NM_204142	−3.397	0.0232	signaling
TBC1D20	XM_001235014	−3.154	0.0389	signaling
HTT	XM_420822	−3.086	0.0065	signaling
P2RY2	XM_425667	−3.079	0.0294	signaling
GRAP2	XM_001234081	−2.967	0.0402	signaling
INPP4A	XM_416886	−2.749	0.0130	signaling
ITSN1	XM_416715	−2.744	0.0139	signaling
EDNRB2	NM_204120	−2.644	0.0317	signaling
FGB	XM_420369	−2.631	0.0194	signaling
RXFP3	XM_429217	−2.616	0.0048	signaling
LOC420403	XM_418509	−2.531	0.0049	signaling
SSTR3	NM_001024583	−2.501	0.0203	signaling
GPR39	NM_001080105	−2.376	0.0087	signaling
GPR97	XM_413998	−2.356	0.0230	signaling
FGF14	NM_204777	−2.322	0.0175	signaling
SFRP4	XM_418831	−2.204	0.0162	signaling
LOC430333	XM_001235474	−1.880	0.0153	signaling
GARNL1	XM_421244	−1.850	0.0477	signaling
RAPH1	XM_421961	−1.842	0.0035	signaling
C14orf138	XM_421460	−1.756	0.0012	signaling
LOC421876	XM_419893	−1.750	0.0333	signaling
NLE1	XM_415857	−1.735	0.0463	signaling
CHRM2	NM_001030765	−1.721	0.0310	signaling

Table 3. iii)

DE-transcripts related with signaling, ion transportation, extracellular matrix protein, and carbohydrate metabolism or post-translation glycosylation modification

Gene symbol	Transcript ID	Log2 units (strong VS weak)	p-value	Category
LOC768958	XM_001232128	−1.622	0.0278	signaling
RASL10B	XM_001233673	−1.546	0.0473	signaling
SIPA1L2	XM_419564	−1.545	0.0020	signaling
LOC769317	XM_001231944	−1.531	0.0287	signaling
PLXDC2	XM_418613	−1.503	0.0155	signaling
SPOCK1	XM_414622	−1.491	0.0375	signaling
CSF2RB	XM_001234608	−1.468	0.0498	signaling
LOC429163	XM_426718	−1.434	0.0220	signaling
PLXNC1	XM_416143	−1.398	0.0103	signaling
RCJMB04_19g9	XM_419989	−1.339	0.0019	signaling
PLXNA1	XM_414370	−1.305	0.0467	signaling
TSPAN5	XM_420654	−1.277	0.0142	signaling
LOC431251	NM_001127171	−1.226	0.0414	signaling
ANXA10	XM_001233661	−1.176	0.0441	signaling
RCJMB04_18c11	NM_001012909	−1.153	0.0436	signaling
SPRED2	XM_419341	−1.124	0.0146	signaling
MOBKL1A	XM_420601	−1.116	0.0026	signaling
ALS2CL	XR_026875	−1.038	0.0119	signaling
MPP1	NM_001007917	−1.036	0.0424	signaling
FGF12	NM_204888	−1.016	0.0006	signaling
GNA13	XM_415686	−1.007	0.0195	signaling
ARL10	XM_414552	−0.975	0.0492	signaling
ADRA2B	XM_425203	−0.885	0.0415	signaling
CCKAR	NM_001081501	−0.868	0.0352	signaling
RCJMB04_3n15	NM_001030902	−0.815	0.0373	signaling
TBC1D24	XM_001232296	4.629	0.0023	IT
SCN9A	XM_422021	3.745	0.0047	IT
KCNT2	XM_426614	3.077	0.0015	IT
LOC395893		3.030	0.0392	IT
NIPAL4	XM_414566	2.796	0.0163	IT
ATP6V0A4	NM_001080102	2.790	0.0397	IT
GRIN2B	XM_416204	2.781	0.0291	IT
KCNJ1	XM_425795	2.363	0.0011	IT
POR	XM_415768	2.294	0.0047	IT
KCNK2	XM_001234269	2.270	0.0447	IT
KIRREL3	XR_026874	2.143	0.0193	IT
GABRG2	NM_205345	1.985	0.0176	IT
SLC4A1	NM_205522	1.794	0.0041	IT
NDUFA7	XM_418185	0.995	0.0464	IT
JPH3	XM_414192	0.925	0.0193	IT
CACNA2D1	XM_001231265	0.852	0.0429	IT
EFCAB5	XM_415833	−3.593	0.0342	IT
SPATA22	XM_001235167	−3.590	0.0176	IT
RCJMB04_1f1	NM_001031133	−3.197	0.0328	IT
LOC428404	XM_425965	−3.177	0.0442	IT
ATP13A3	XM_422709	−3.059	0.0065	IT
GABRB2	XM_001232377	−2.963	0.0472	IT
SLC8A1	NM_001079473	−2.694	0.0342	IT
CACNA2D3	XM_414338	−2.622	0.0100	IT
SERINC5	XM_424762	−2.615	0.0129	IT

Table 3. iv)

DE-transcripts related with signaling, ion transportation, extracellular matrix protein, and carbohydrate metabolism or post-translation glycosylation modification

Gene symbol	Transcript ID	Log2 units (strong VS weak)	p-value	Category
LOC421866	XR_027148	−2.613	0.0307	IT
LOC425295	XM_423073	−2.586	0.0212	IT
LOC772391	XM_001235535	−2.573	0.0163	IT
KCNK17	XM_419477	−2.506	0.0047	IT
KCTD16	XM_425217	−2.363	0.0076	IT
CNNM1	XM_421703	−2.309	0.0446	IT
KCNJ5	XM_417864	−2.138	0.0324	IT
RCJMB04_11e10	NM_001030630	−2.102	0.0226	IT
GRIN3A	XM_001232181	−1.680	0.0355	IT
ATP6V1A	NM_204974	−1.202	0.0065	IT
RCJMB04_16a12	NM_001031305	−1.030	0.0261	IT
CNGA3	NM_205221	−0.962	0.0479	IT
P2RX4	NM_204291	−0.882	0.0252	IT
MEGF10	XM_424719	3.910	0.0398	EM
FAT2	XM_414584	3.892	0.0052	EM
SDK2	NM_204538	2.905	0.0097	EM
NRXN3	XM_421297	2.643	0.0139	EM
LAMA4	XM_419780	2.569	0.0322	EM
NTNG1	XM_001231446	2.004	0.0328	EM
CRTAC1	NM_001080211	1.930	0.0132	EM
LAMC1	NM_204166	1.680	0.0257	EM
LAMB4	XM_001232877	1.642	0.0161	EM
PPFIA1	XM_421074	1.583	0.0166	EM
CLDN20	XM_001232002	1.438	0.0101	EM
CDH9	XM_001231501	1.296	0.0311	EM
CHAD	XM_416236	1.268	0.0240	EM
LOC396026	NM_205128	1.239	0.0403	EM
PCDH21	NM_001001759	1.174	0.0187	EM
EPDR1	XM_418830	1.158	0.0469	EM
CPNE8	XM_001231388	1.121	0.0219	EM
COL12A1	NM_205021	0.996	0.0003	EM
PKP2	XM_416362	0.983	0.0447	EM
CD72	NM_205052	0.855	0.0263	EM
NINJ2	XM_416382	−4.039	0.0409	EM
OTOF	XM_420015	−3.980	0.0016	EM
COL13A1	XM_001232218	−3.260	0.0006	EM
SVEP1	XM_424917	−2.967	0.0250	EM
OTOP1	XM_420790	−2.830	0.0263	EM
GPNMB	XM_425991	−2.771	0.0327	EM
PKP1	XM_419240	−2.349	0.0324	EM
LAMA2	XM_419746	−2.336	0.0000	EM
COL8A2	XM_425780	−2.295	0.0170	EM
FNBP4	XM_424260	−2.279	0.0486	EM
EGFL6	XM_416835	−2.268	0.0392	EM
CDH18	XM_426046	−1.955	0.0090	EM
CLDN8	XM_425544	−1.598	0.0288	EM
RCJMB04_34k20	NM_001031112	−1.214	0.0085	EM
SRPX	XM_416781	−1.123	0.0244	EM
DLG1	XM_422701	−1.084	0.0171	EM
FBLN1	NM_204165	−0.981	0.0206	EM

Table 3. v)

DE-transcripts related with signaling, ion transportation, extracellular matrix protein, and carbohydrate metabolism or post-translation glycosylation modification

Gene symbol	Transcript ID	Log2 units (strong VS weak)	p-value	Category
F13A1	NM_204685	−0.949	0.0234	EM
FREM1	XM_424932	−0.910	0.0057	EM
MEGF10	XM_424719	3.910	0.0398	EM
MGAT4C	XM_425447	3.007	0.0497	GM or CM
CHST3	NM_205121	2.839	0.0372	GM or CM
EDEM3	XM_422293	2.539	0.0179	GM or CM
LARGE	NM_001004404	2.071	0.0429	GM or CM
GFPT2	XM_424573	1.919	0.0173	GM or CM
GALNTL1	XM_001231964	1.895	0.0452	GM or CM
WDR77	NM_001030916	1.805	0.0445	GM or CM
NDST3	XM_426325	1.403	0.0121	GM or CM
B3GALT1	XM_426584	1.254	0.0483	GM or CM
OGDHL	XM_421503	1.144	0.0022	GM or CM
MAN1A2	XM_416490	−4.700	0.0087	GM or CM
POFUT2	XM_421892	−2.804	0.0156	GM or CM
LOC772154	XM_001235329	−2.798	0.0006	GM or CM
KLB	XM_423224	−2.335	0.0380	GM or CM
TRIM7.2	NM_001099354	−1.986	0.0041	GM or CM
LOC771361	XM_001234647	−1.433	0.0369	GM or CM
RCJMB04_28l23	NM_001039316	−1.324	0.0034	GM or CM
NDST4	XM_420638	−1.239	0.0261	GM or CM
PFKM	NM_204223	−1.170	0.0460	GM or CM
NUP153	XM_418937	−1.109	0.0378	GM or CM
GPD1L	XM_418763	−0.947	0.0383	GM or CM
PHKA2	XM_416811	−0.926	0.0326	GM or CM
B3GNTL1	XM_415599	−0.883	0.0109	GM or CM
PMM1	XM_416228	−0.787	0.0490	GM or CM
MMP11	XM_001232776	2.209	0.0391	GM or CM
ST3GAL4	XM_417860	−1.043	0.0094	GM or CM

IT represents ion/proton transporter, EM represents extracellular matrix, GM represents post-translation glycosylation modification, and CM represents carbohydrate metabolism.