Paradigm of Time-sequence Development of the Intestine of Suckling Piglets with Microarray.

The interaction of the genes involved in intestinal development is the molecular basis of the regulatory mechanisms of intestinal development. The objective of this study was to identify the significant pathways and key genes that regulate intestinal development in Landrace piglets, and elucidate their rules of operation. The differential expression of genes related to intestinal development during suckling time was investigated using a porcine genome array. Time sequence profiles were analyzed for the differentially expressed genes to obtain significant expression profiles. Subsequently, the most significant profiles were assayed using Gene Ontology categories, pathway analysis, network analysis, and analysis of gene co-expression to unveil the main biological processes, the significant pathways, and the effective genes, respectively. In addition, quantitative real-time PCR was carried out to verify the reliability of the results of the analysis of the array. The results showed that more than 8000 differential expression transcripts were identified using microarray technology. Among the 30 significant obtained model profiles, profiles 66 and 13 were the most significant. Analysis of profiles 66 and 13 indicated that they were mainly involved in immunity, metabolism, and cell division or proliferation. Among the most effective genes in these two profiles, CN161469, which is similar to methylcrotonoyl-Coenzyme A carboxylase 2 (beta), and U89949.1, which encodes a folate binding protein, had a crucial influence on the co-expression network.

INTRODUCTION

The pig (Sus scrofa) has been widely used as a model for studying the development of human organs and disease progression, owing to the similarity in size and physiology of the two species (Butler et al., 2007; Lunney, 2007; Ojeda et al., 2008). The Landrace is one of the pig breeds that are used world-wide in research and production because of its genetic stability and high lean percentage.

The alimentary tract is essential for absorption of nutrients and the function of the immune system. A notable feature of the intestine is the complex changes that occur over the lifetime of an individual (Weaver et al., 1991; Thompson et al., 2008). The development of the intestine is regulated by a large number of factors, including nutrients, the micro flora, the epithelium, innate and adaptive immunity (Pacha, 2000; Caicedo et al., 2005; Donovan, 2006; Commare et al., 2007; Ojeda et al., 2008), which constitute the intestinal ecosystem. It is known that the adaptive change involved in enteric development is highly conserved and is a complex process (Gilbert and Lloyd, 2000). The underlying molecular mechanism of intestinal development is associated with genes and the corresponding proteins that participate in this change, and with the interactions and signal transduction processes of the gene and/or protein networks. It is almost impossible. However, to understand the complex interactions of genes and proteins by studying them in an individual and static manner. The study of gene regulatory networks makes it feasible to analyze the interactions of genes in an integral and dynamic way (Tomita et al., 1999; Gilbert and Lloyd, 2000). Gene regulatory networks have attracted considerable research interest owing to the rapid accumulation of genomic information (Kitano, 202). Comprehensive assessment of the gene expression profiles needs to be carried out in order to identify gene networks (Kitano, 2002), and the high-throughput techniques of molecular biology can be used to provide insight into such networks.

Microarray technologies have enabled researchers to explore the dynamics of transcription on a genome-wide scale (Young, 2000; Jiang and Deyholos, 2006; Zhang et al., 2008), and they provide a systematic insight into the crucial factors that control development. They can be used to elucidate the mechanisms of intestinal development at a molecular level. Studies of gene expression profiling using microarray technologies in the development of plants and the human fetal pancreas have made significant progress (Sarkar et al., 2008; Wechter et al., 2008), but there have been no reports of such studies on the development of the porcine intestine. In the present study, we aimed using microarray technology to identify the significant pathways and key genes that regulate the intestinal development of Landrace piglets.

MATERIALS AND METHODS

Animals and tissue preparation

All animals used in this study were managed humanely according to the established guidelines of the Animal Care Committee, Sichuan Agricultural University, P. R. China. The four Landrace sows were mated with a boar of the same breed. Two male piglets were sacrificed at 0, 3, 8, 14 and 21 d after birth, respectively. The mid-section of the jejunum of each piglet was dissected rapidly. The samples were snap-frozen in liquid nitrogen and stored at −80°C until further use.

Microarray analysis

Isolation of RNA was performed with the QIAGEN RNeasy mini kit (QIAGEN, Valencia, CA, USA) following the manufacturer’s instructions. The RNA samples were then analyzed using the Affymetrix GeneChip® Porcine Genome Array (Affymetrix, Santa Clara, CA, USA), which contains 23,937 probes sets that interrogate approximately 23,256 transcripts from 20,201 genes of S. scrofa. Statistics were computed using the GeneChip Operating Software (GCOS). The data discussed in this article have been deposited in the NCBIs Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) and are accessible through GEO Series accession number GSE13456.

Time sequence profile analysis of gene expression

The differentially expressed genes from four time points (d 3, 8, 14, and 21) were compared with those of d 0 using Affymetrix GCOS software. The raw values for gene expression were converted to log2ratio. According to the analysis strategy for an experiment representing a cluster of short time series (Ernst et al., 2005), 80 profiles were defined in the experiment. The profiles were analyzed that had the most significant probability assigned in comparison with that expected by the Fisher exact test and the multiple comparison test. If the p-value of the profiles was less than 0.01, the gene expression profiles were identified to be significantly different with respect to time.

Functional analysis of significant profiles of co-expressed genes

According to previous studies (Yu et al., 2003; Pan et al., 2008; Rawat et al., 2008), co-expressed genes have similar functions. The important functions were analyzed for those profiles found to be significant using Gene Ontology (GO), which organizes genes into hierarchical categories based on biological processes, molecular functions, and cellular components (Ashburner et al., 2000). To undertake GO analysis, the p-values from the Fisher’s exact test for over-representation of the selected genes in all GO biological categories are computed, as described previously. Significance of the GO was accepted when the p-value was less than 0.05. Functional annotation can provide a detailed description of the roles of genes located in deeper categories of GO, according to the hierarchical organization of GO. Therefore, the enrichment value was defined to filter the deeper categories of GO: the higher the enrichment value of GO, the more target GO were expected. Within a given category, the enrichment Re is given by equation,

where nf is the number of differentially expressed genes within the particular category, n is the total number of genes within the same category, Nf is the number of different genes on the entire microarray, and N is the total number of genes on the microarray.

Pathway analysis of differentially expressed genes

A pathway is defined as a functional unit in which genes interact with each other and achieve a process of signal transduction. To explain test phenotypic features with respect to development processes, significant pathways are identified from differentially expressed genes. This is based on the p-values of the Fisher’s exact test for over-representation of the selected genes in all pathway categories, as described perviously (Shalgi et al., 2007). A significant pathway was specified when the p-value was less than 0.05. To assess the significance of a particular category with respect to random chance, the false discovery rate (FDR) was estimated for a set of categories. After resampling 5,000 times, the FDR was defined by equation, where Nk refers to the number of Fisher’s test p-values less

than ×2 test p-values.

The p-values obtained using Fisher’s exact test are lower than ×2 test p-values. Pathway enrichment analysis was performed according to a method reported previously (Kanehisa and Goto, 2000) to investigate whether particular functions were enriched among the genes that control distinctive characters in groups of genes that were differentially expressed on d 3, 8, 14, and 21, in comparison with d 0 (the control time point). With increasing enrichment, the corresponding function becomes more significant. Within a given category, the enrichment Re is given by equation, where nf is the number of different genes within the particular pathway category, and n is the total number of genes within the same category, Nf is the number of different genes on the entire microarray, and N is the total number of genes on the microarray.

Construction and topological analysis of the gene co-expression network

We transformed the normalized expression values of profile 66 and profile 13 from the individual correlation analysis into measures of pair-wise connection strengths for the construction of a co-expression network (Carlson et al., 2006). The network edges were specified to be more than 0.9, which is the value of the coefficient correlation that incidates a strong relationship, such as co-expression of genes. The centrality of a network is represented by the central degree (Barabasi and Oltvai, 2004). It is possible to determine the characteristic variables that are related to the distances among genes. The core of maximum order is referred to as the main core, or the highest k-core of the graph (Huber et al., 2007). A k-core subgraph of a graph can be generated by deleting the vertices recursively from the graph whose degree is less than k (Huber et al., 2007). Larger values of “coreness” clearly correspond to vertices of larger degree and a more central position in the structure of the network. In this work, we applied the notion of the k-core subgraph to predict similarity of gene function (Altaf-Ul-Amin et al., 2006). The nodes labeled by the same color may represent similar gene ontology terms. The highlight ontology of nodes of the same color was assessed by counting the genes of the same GO that had the same color.

Real-time quantitative PCR

Fourteen differentially expressed genes (AREG (Amphiregulin), ATRN (Attractin), CD59 (CD59 molecule), CLU (Clusterin), CST3 (Cystatin C), DPEP1 (Dipeptidase), EPHX1 (Epoxide hydrolase), LEAP2 (Liver expressed antimicrobial peptide 2), LOC396871 (Arginine rich antibacterial peptides), MS4A2 (Membrane-spanning 4-domains), PIAP (Putative inhibitor of apoptosis), RETN (Resistin), UBP (Ubiquitin-specific protease), FBP (Folate binding protein) dentified in the array experiment were selected and analyzed by RT-qPCR. Aliquots of the same RNA samples used for the microarray analysis were used for the real-time PCR. The reverse transcriptions (RT) were carried out with the total cellar RNA using TRIZOL Reagent® (Invitrogen California, USA). Each sample of cDNA was diluted 1:4 fold in sterile ddH2O, and 1 μl of this dilution was used as the template for qPCR. Primers for the PCR reactions were designed by DNAStar (DNASTAR Inc) to have a Tm of 62 to 65°C and an optimal annealing temperature of 63°C, with a length of 100 and 250 bp. Real-time PCR was performed with SYBR® Premix Ex TaqTM (TaKaRa, Dalian, China) in 20 μl reaction volumes using the iCycler System (BioRad, Richmond, CA, USA) according to the manufacturer’s instructions. Each PCR reaction contained 1 μl of cDNA, 0.2 μM of each of the primers (Table 1) and 10 μl of PCR master mix. The initial denaturing time was 95°C per minute, followed by 40 cycles consisting of 94°C for 15 s, 63°C for 20 s, and 72°C for 30 s with a single measurement of fluorescence.

Table 1

Primer sequences and PCR product sizes of genes selected for validation by RT

Gene	GenBank	ID	Primer sequence	Annealing tm(°C)	Product size(bp)
β-actin	SSU07786		Forward5’-GCTGGCCGGGACCTGACCGACT-3’Reverse5’-TTGGCATAATTGTTACACGTTTGG-3’	63.5	154
Amphiregulin AREG	AY028311.1		Forward5’-GGCGCCCGTGGTGCTGTCACTC-3’Reverse5’-CGCTAGCAGGGGGAGCCTCACTT-3’	61.5	179
Clusterin CLU	NM_213971.1		Forward5’-AAGTCCCGCTTCGCCCGGAACATC-3’Reverse5’-CGCTAGCAGGGGGAGCCTCACTT-3’	61.4	225
Dipeptidase DPEP1	NM_214108		Forward5’-CAACACGCCCTGGGCCGACAACTG-3’Reverse5’-ACGGGGGCCTGGGACAGCTTCAGA-3’	64.1	183
Epoxide hydrolase EPHX1	NM_214355		Forward5’-CCCGTCCCCCTGGCCTACAT-3’Reverse5’-AAAGCCGGTGGGCACGTGGACCTT-3’	62.3	243
Liver expressed antimicrobial peptide 2 LEAP2	NM_213788		Forward5’-ACGGCTCCCCGATACCGGAACGAG-3’Reverse5’-TCCTGGGCCACACTTAGGGAACAGC-3’	59.9	169
Arginine rich antibacterial peptides LOC396871	NM_214355.1		Forward5’-TGCCGGAACATCTATGAGAATGACT-3Reverse 5’-AATCAGGGGCAGAAGCTATCTCCTT-3’	57	250
Membrane-spanning 4-domains MS4A2	NM_214092.1		Forward5’-TGCCGGAACATCTATGAGAATGACT-3’Reverse5’-AATCAGGGGCAGAAGCTATCTCCTT-3’	54.8	244
Putative inhibitor of apoptosis PIAP	NM_214181		Forward5’-GCAGCCCGCTTTAAAACATTCTGTA-3’Reverse 5’-TTGAACGCGACTGATGAACTCCTGT-3’	56.1	240
Resistin RETN	NM_213783.1		Forward5’-GGGGCTGCTGGTGTGGGGCAAGT-3’Reverse5’-CGCGCACGTCCCAAGAGCCACAG-3’Forward5’-GGGCCCGGTCGGTTTGCACAACAT-3’	63	218
Ubiquitin-specific protease UBP	NM_213826.1		Reverse5’-CACCGCCTTCTGCCGGCTGTCCT-3’Forward5’-AAGCCGTGCCACCCAAAACTTACTA-3’	61	199
CD59 molecule CD59	NM_214170		Reverse5’-CAGCAGAGCGGTTTTCCCTGATGA-3’Forward 5’-TGGCCAGTCGTCCGTTTGCCTCTG-3’	56.2	167
Attractin ATRN	AB271953		Reverse 5’-AGGGCGCTGGCCACTGCGAGAC-3’Forward: 5’-ACGGGGACCACCTGCGACACATT-3’	61.5	220
Folate binding protein FBP	NM_213853		Reverse: 5’-CGCACTCAAGGCCAAGGCGTAGA-3’Forward: 5’-CCCAACAGGCAAGGCAGCTAAAGC-3’	66	1,320
Cystatin C CST3	NM_001044602		Reverse: 5’-CTGCCAAGACCGCGTGGTTGAGATT-3’	59.3	203

RESULTS

To characterize the genes involved in intestinal development during the first 21 d in piglets’ life, total RNA was extracted from the mid-section of the jejunum on postnatal d 0, 3, 8, 14 and 21. Gene expression was analyzed using a porcine genome array and statistics were computed using GCOS. In total, 8,427 significant transcripts were identified. Among these statistically significant expression ratios, probes representing 4,789 unique genes increased in transcript abundance at least 2.0 fold at one or more time points, while probes representing 3,618 unique genes decreased in abundance at least 2.0 fold. Thus, at least 24% of the total number of transcripts from S. scrofa were strongly (i.e. >2.0 fold) induced by intestinal development; however, at least 18% of transcripts were repressed by an equivalent amount (Table 2).

Table 2

Summary of data obtained from the genechip analysis

	d 3 vs d 0	d 8 vs d 0	d 14 vs d 0	d 21 vs d 0
Total transcripts	20,201	20,201	20,201	20,201
Differentially expressed genes (DEG)	1,881	1,918	2,204	2,404
DEG/total (%)	9.31	9.49	10.91	11.9
Increased genes (IG)	1,117	994	1340	1338
IG/DEG (%)	59.38	51.82	60.8	55.66
Decreased genes (DG)	764	924	864	1,066
DG/DEG (%)	40.62	48.18	39.2	44.34

Gene expression profiles significant to the development process

During the course of intestinal development from d 0 to 21, 80 model profiles were identified in which we defined the value of maximum unit change to be 1. For these model profiles, 30 were identified as significant (Figure 1). Among them, the two most significant profiles were profile 66 and profile 13, according to ascending order of the p-values (Figure 2). If the p-value of a profile was less than 0.01, the gene profile was identified as one that changed significantly with time. According to this principle, the lowest p-value profiles may represent the regulatory genes that affect the development of the porcine 539 genes, the expression of which increased significantly at d 3 and maintained a similar level of expression until d 21. In contrast, profile 13 had a significance level of 8.8×10−87 and contained 351 genes, the expression of which decreased significantly at d 3 and maintained a similar level of expression until d 21.

Figure 1

Model expression profiles of differentially expressed genes. The number in the top left hand corner of a profile box is the model profile ID number, which has no special meaning; the lower left hand corner contains the p-value. Model expression profiles are ordered from left to right and from top to bottom according to gradually increasing significance. Significant profiles have statistical significance (corrected with the multiple permutation test) less than 0.01.

Figure 2

Gene clusters of significant expression profiles 66 and 13. The horizontal axis represents time, and the vertical axis shows the time series of the level of gene expression after log normal transformation. Each diagram contains the model profile ID at the top center. “Genes assigned” represents the number of genes actually assigned to the profile; “Genes expected” represents the expected number of genes based on a permutation test of time points; p-value represents the significance of the number of genes assigned versus expected.

Significant gene functions and genes of significant profiles

Results of the analysis of GO categories for profile 66 implied that the main enriched GO categories were focused on control checkpoints for cell size, cell wall biogenesis, deoxyribonucleotide catabolic processes, fungal-type cell wall biogenesis, fungal-type cell wall organization and biogenesis, lymphocyte chemotaxis, M phase-specific microtubule processes, membrane depolarization, microvillus biogenesis, mitochondrial depolarization, mitotic metaphase, etc. during growth from 0 to 21 d (Table 3). To further investigate the most significant functions, the annotation of the GO was analyzed using Amigo (a website tool (http://amigo.geneontology.org/cgi-bin/amigo/go.cgi). It was concluded that the process of the cell cycle initiated at a checkpoint and the immune response to fungi stimulated by fungal infection were the most significant phenomena that occurred from d 3 to d 21.

Table 3

Significant GO of profile 66

Go ontology	p-value	Enrichment
Cell size control checkpoint	0	5.94
Cell wall biogenesis	0	5.94
Deoxyribonucleotide catabolic process	0	5.94
Fungal-type cell wall biogenesis	0	5.94
Fungal-type cell wall organization and biogenesis	0	5.94
Lymphocyte chemotaxis	0	5.94
M phase specific microtubule process	0	5.94
Membrane depolarization	0	5.94
Microvillus biogenesis	0	5.94
Mitochondrial depolarization	0	5.94
Mitotic metaphase	0	5.94
Negative regulation of calcium-dependent cell-cell adhesion	0	5.94
Negative regulation of megakaryocyte differentiation	0	5.94
Negative regulation of mitochondrial depolarization	0	5.94
Phosphatidic acid biosynthetic process	0	5.94
Phosphatidic acid metabolic process	0	5.94
Pyrimidine deoxyribonucleotide catabolic process	0	5.94
Regulation of calcium-dependent cell-cell adhesion	0	5.94
Regulation of mitochondrial depolarization	0	5.94
Regulation of mitochondrial membrane potential	0	5.94
Release of cytochrome c from mitochondria	0	5.94
Deoxyribonucleoside monophosphate biosynthetic process	0	3.96
Deoxyribonucleoside monophosphate metabolic process	0	3.96
Leukocyte mediated cytotoxicity	0	3.96
Negative regulation of leukocyte mediated cytotoxicity	0	3.96
Negative regulation of T cell mediated cytotoxicity	0	3.96
Regulation of cell killing	0	3.96

In contrast, the main enriched GO categories for profile 13 were focused on folic acid transport, negative regulation of endocytosis, the phosphatidylcholine biosynthetic process, phosphatidylcholine metabolic processes, the response to cobalt, manganese and zinc ions, catabolic processes for aminoglycans, glycosaminoglycans, and leucine, leucine metabolic processes, a lipopolysaccharide-mediated signaling pathway, Mo-molybdopterin cofactor biosynthetic processes, and Mo-molybdopterin cofactor metabolic processes, etc., as presented in Table 4.

Table 4

Significant GO of profile 13

Go ontology	p-value	Enrichment
Phosphatidylcholine biosynthetic process	0	5.94
Folic acid transport	0	5.94
Response to manganese ion	0	5.94
Response to cobalt ion	0	5.94
Negative regulation of endocytosis	0	5.94
Response to zinc ion	0	5.94
Phosphatidylcholine metabolic process	0	5.94
Negative regulation of chemokine biosynthetic process	0	2.97
Aminoglycan catabolic process	0	2.97
Lipopolysaccharide-mediated signaling pathway	0	2.97
Glycosaminoglycan catabolic process	0	2.97
Regulation of lipopolysaccharide-mediated signaling pathway	0	2.97
Vitamin transport	0	2.97
Molybdopterin cofactor biosynthetic process	0	2.97
Leucine catabolic process	0	2.97
Mo-molybdopterin cofactor metabolic process	0	2.97
Negative regulation of lipopolysaccharide-mediated signaling pathway	0	2.97
Negative regulation of JNK activity	0	2.97
Molybdopterin cofactor metabolic process	0	2.97
Leucine metabolic process	0	2.97
Mo-molybdopterin cofactor biosynthetic process	0	2.97
Cellular zinc ion homeostasis	0	2.38

Distribution of significant pathways of differentially expressed genes

The pathway was analyzed using the KEGG and Genmapp databases. The significant pathways for differentially expressed genes were applied to the Fisher exact test and tests of multiple comparisons to obtain the significant functions. From the results of Table 6, it can be seen that the metabolism of xenobiotics by cytochrome P450, Alzheimer’s disease, antigen processing and presentation, the PPAR signaling pathway, complement and coagulation cascades, cell adhesion molecules (CAMs), autoimmune thyroid disease, allograft rejection, graft-versus-host disease, natural killer cell mediated cytotoxicity, and cytokine-cytokine receptor interaction were the significant pathway categories. It is well known that owing to the changes in diet and living environment that occur after the birth of piglets, adaptive changes occur in the internal environment. In this process various immune reactions are activated, such as antigen processing, complement cascades, allograft rejection, graft-versus-host disease, natural killer cell mediated cytotoxicity and cytokine-cytokine receptor interaction. It can be concluded. Therefore, that the immune response may be coupled with the growth of the organism (Table 5).

Table 6

Network structure parameter lists for each gene

Probe Set ID	NCBI-geneID	Clustering coefficient	Degree	k-core
Ssc.14544.1.S1_at	U89949.1	0.58	28	12
Ssc.21579.1.S1_at	CF789025	0.58	27	12
Ssc.1735.1.S1_at	BI118904	0.62	26	12
Ssc.5250.1.S1_at	BQ604261	0.58	26	12
Ssc.15069.1.S1_at	BI345060	0.63	25	12
Ssc.22910.1.A1_at	CO986224	0.67	24	12
Ssc.10748.1.A1_at	BQ597862	0.61	24	12
Ssc.16817.1.S1_at	BE233328	0.64	23	12
Ssc.25002.1.S1_at	BG833819	0.69	22	12
Ssc.6080.1.S1_at	NM_213870.1	0.68	22	12
Ssc.15430.1.A1_at	CO937301	0.67	22	12
Ssc.11246.1.A1_at	BI181165	0.55	22	12
Ssc.18215.1.S1_at	NM_214365.1	0.71	21	12
Ssc.13713.1.A1_at	BQ603271	0.62	20	12
Ssc.3714.1.S1_a_at	CK456589	0.71	19	12
Ssc.21141.1.S1_at	BX666500	0.6	19	12
Ssc.26028.1.S1_at	BX926563	0.43	19	9
Ssc.1209.1.S1_at	CK465976	0.76	18	12
Ssc.18033.1.A1_at	CF180565	0.75	18	12
Ssc.5052.1.S1_at	BI398918	0.69	18	12
Ssc.2795.1.S1_at	BI183754	0.68	18	12
Ssc.21952.1.S1_at	BX666921	0.54	18	10
Ssc.9503.1A1_at	CF789586	0.46	18	9
Ssc.24263.2.S1_at	AW480085	0.55	16	9
Ssc.12622.2.A1_a_at	CO948012	0.52	16	9
Ssc.27572.1.S1_at	CN161469	0.67	15	11
Ssc.10988.1.S1_at	CO947173	0.72	15	10

Representative public ID can be queried on NCBI. Clustering coefficient and degree describe the network properties of genes. K-core represents the subnetwork in which genes are located.

Table 5

Significant pathway categories (p<0.01, FDR<0.05)

Pathway	p-value	Enrichment	FDR
Metabolism of xenobiotics by cytochrome P450	0	4.25	0
Alzheimer’s disease	0	4.12	0
Antigen processing and presentation	0	4.1	0
PPAR signaling pathway	0	3.89	0
Complement and coagulation cascades	0	3.57	0
Cell adhesion molecules (CAMs)	0	3.47	0
Autoimmune thyroid disease	0	3.36	0
Allograft rejection	0	3.36	0
Graft-versus-host disease	0	3.27	0
Natural killer cell mediated cytotoxicity	0	2.7	0
Cytokine-cytokine receptor interaction	0	2.19	0.01

The gene co-expression network

The gene co-expression network was constructed with respect to gene function associations. Given that the elements of a network represent a variety of gene regulatory abilities, a large-scale gene network can be divided into subgraphs, named k-core networks, in which all genes affect each other (Huber et al., 2007). Generally, the complexity of gene relationships increases with the rank of the k-core value. In the light of the definition of the k-core network, the core status within a large-scale gene network consists of subgraphs of higher k value. In this work, we used the k-cores of protein-protein interaction networks to define the main gene functions in the main subgraph. Moreover, we aimed to find the main GO assigned by the maximum number of genes in a separate k-core and then to define the key gene functions at each level of the network.

It has been reported that the core functions at the core status of a network have a top k-core level (Altaf-Ul-Amin et al., 2006), and a higher degree but lower clustering coefficient (Barabasi and Oltvai, 2004). A k-core subnetwork with a k-value of 12 was located at the core status within the large scale gene network (Figure 3). Using the hierarchical categories of gene ontology, it could be concluded that most of the genes in this subnetwork are attributed to cell communication, which may be related to cell division or proliferation through cell adhesion and transport of nutrients (for details of the genes, Tables 3 and 4). In addition, most genes in the secondary network (with a k-core value of 11) play an important role in metabolic processes of nitrogen compounds and amines (for detailed information, Tables 3 and 4), except for those genes without detailed functional annotation. Among these genes, CN161469 is similar to methylcrotonoyl-Coenzyme A carboxylase 2 (beta) when analyzed by BLAST, and it takes part in the regulation of other genes within the secondary level subnetwork because of the higher value of its degree. Most notably, the gene U89949.1, which encodes a folate binding protein, possessed the highest k-core and degree, and the lowest clustering coefficient (Table 6). These features indicate that this is the core gene that is located at the center of both the large scale network and the 12 k-core subnetwork, and it regulates directly 28 neighboring genes.

Figure 3

Gene association network. Nodes denote genes; undirected links represent gene-gene interrelation; the size of the nodes represents the power of the interrelationship among the nodes. Different colors of nodes represent distinct k-core gene classes.

Validation of the gene chip data by RT-PCR

To confirm whether the genes identified were differentially expressed, quantitative RT-PCR analysis for specific transcripts was carried out as described previously. Based on our microarray analysis, 14 genes, including three antimicrobial peptides (LEAP2, LOC396871, and MS4A2), two proteins involved in the regulation of proliferation and apoptosis of epithelial cells (PIAP and AREG), and genes associated with a number of physiological and metabolic processes (ATRN, CD59, CLU, DPEP1, EPHX1, RETN and UBP) were selected. The core gene FBP was also involved in the RT-PCR analysis. For all of these genes, the expression ratios measured by RT-qPCR and by microarray were highly correlated (r = 0.74, 0.73, 0.95, 0.93 at d 3, 8, 14, 21 respectively) (Table 7). Genes involved in apoptosis, signal transduction, and energy and protein metabolism, such as CLU, MS4A2 and UBP, were up-regulated during suckling, as indicated by both the Genechip and RT-qPCR experiments. However, HCRTR2 displayed significant differential expression when detected by RT-qPCR but no significance on Genechip testing. The correlation of the results from the RT-qPCR and Genechip testing helped us to analyze the reliability of the assays so as to reveal true differences in the gene expression profiles in the intestine.

Table 7

Genes differentially expressed in the Genechip data

Symbol	Method	0-d	3-d old	8-d old	14-d old	21-d old
Clusterin	RT-PCR	1	0.584	0.822	8.830	10.255
CLU	Genechip	1	1.546	2.923	8.558	13.860
Amphiregulin	RT-PCR	1	2.263	1.004	1.546	1.969
AREG	Genechip	1	10.317	1.950	5.347	7.130
Folate binding	RT-PCR	1	1.005	0.380	0.850	3.235
FBP	Genechip	1	0.295	0.304	0.264	0.2606
Dipeptidase	RT-PCR	1	0.115	0.111	0.179	0.251
DPEP1	Genechip	1	0.255	0.248	0.375	0.468
Membrane-spanning 4-domains	RT-PCR	1	0.475	2.670	31.138	7.704
MS4A2	Genechip	1	0.831	1.798	8.877	4.304
Utative inhibitor of apoptosis	RT-PCR	1	0.905	1.449	6.358	2.027
PIAP	Genechip	1	2.366	2.106	4.854	3.761
Ubiquitin-specific protease	RT-PCR	1	5.199	1.770	5.554	28.016
UBP	Genechip	1	4.774	3.007	3.499	16.768
Attractin	RT-PCR	1	0.211	0.477	0.758	0.197
ATRN	Genechip	1	0.239	0.679	0.186	0.361
Liver expressed antimicrobial peptide 2	RT-PCT	1	1.775	4.712	2.652	2.841
LEAP2	Genechip	1	3.882	3.552	2.484	2.758
Resistin	RT-PCT	1	0.995	2.634	1.1760	0.309
RETN	Genechip	1	2.874	3.181	3.077	2.113
CD59 molecule	RT-PCT	1	0.080	0.099	0.346	0.197
CD59	Genechip	1	1.485	0.848	1.048	0.992
Epoxide hydrolase EPHX1	RT-PCT	1	6.169	22.177	64.977	33.188
	Genechip	1	12.195	26.257	25.856	35.314
Arginine rich antibacterial peptides	RT-PCT	1	4.342	23.007	10.273	1.649
LOC396871	Genechip	1	3.211	3.927	4.822	1.118
Cystatin C	RT-PCR	1	1.588	0.718	2.674	8.490
CST3	Genechip	1	1.474	2.631	2.469	1.996
r		1	0.74	0.72	0.95	0.93

DISCUSSION

Owing to the continuous improvement of microarray technology, the rapid development of bioinformatics and the depth of research on genomics and proteomics, global genome expression profile analysis provides a more comprehensive and dynamic analysis of intestinal development at the genome-wide transcriptional level (Ravasz et al., 2002; Schlitt et al., 2003; Stears et al., 2003; Barabasi and Oltvai, 2004; Han, 2008; Lee et al., 2008). The understanding and description of cells, organs and living individuals on a system level is a new area of biological research (Davidson et al., 2002; Kitano, 2002; Oltvai and Barabasi, 2002). Various kinds of model organism have been proved to be useful in this area (Covington et al., 2008; Sarkar et al., 2008; Schweikl et al., 2008; Wechter et al., 2008). To our knowledge, this is the first study involving whole genome analysis of the genes related to intestinal development in Landrace piglets at different stages of growth (0, 3, 8, 14, and 21 d after birth).

In the present study, more than 8,000 differentially expressed transcripts were identified by microarray technology. This technique was validated by RT-qPCR for utilization in a comprehensive study of gene expression profiles. However, many of these transcripts had no detailed information in the cDNA resources deposited for swine, which limited the analysis of their biological importance. This indicates that many porcine genes have yet to be analyzed. Complex molecular events pertinent to intestinal development occur during the suckling stage in piglets (Hall and Byrne, 1989; Buddington, 1994; Thaler and Cummings, 2008). We analyzed the development process of the gene expression profile, and 80 model profiles were obtained in all. These profiles do not cover all the possible profiles that could be generated in the intestine of Landrace piglets during suckling. Among these profiles, 30 were identified as significant and two of them (profiles 66 and 13) were highly significant. Significant profiles may indicate that common functions were mainly attributed to the co-expressed genes, and these functions influence the biological characteristics of the organism (Gracey et al., 2004). We can therefore conclude, according to the significance of the profiles, that the genes assigned to these two profiles may be highly specific for the intestinal development of piglets in early life. In order to obtain a better understanding of the genes assigned to these two profiles and to elucidate the molecular mechanism of development of the intestine, analysis of GO categories, pathway analysis, and topological analysis of the gene co-expression network were carried out.

Annotation of GO categories has proved to be remarkably useful for the mining of functional and biological significance from the results of microarray studies. In our study, significant GO categories were determined by the hypergeometric distribution and multiple comparison tests. The results obtained for the GO categories implied that the genes in profiles 66 and 13 were involved in many biological processes, among which the distinguished phenomena were initiation of the cell cycle at a checkpoint and the immune response to fungi for profile 66, while the folic acid transport process was most notable for profile 13. These results suggest that the genes in profile 66 may regulate cell proliferation and differentiation, and the genes in profile 13 may modulate the metabolism of nutrients.

Analysis of the pathway for the lists of differentially expressed genes may allow us to find the target pathway that regulates phenotypic differences. Pathway analysis in the present study indicated that the genes in profiles 66 and 13 focused mainly on the immune system, metabolism, cell adhesion molecules, etc. It may be concluded that these genes play an important role in the above pathways, which agrees with the results obtained from the analysis of the GO categories.

During biological processes, a macromolecular network can be constructed according to the yeast two-hybrid method (Y-2H) (Smidtas et al., 2006) or an algorithmic prediction based on the gene function correlation and expression profiles (Nikiforova and Willmitzer, 2007). The complexity of the network model based on the high throughput gene expression test using the algorithmic prediction made it feasible to determine snapshots of protein-protein interactions, gene expression regulatory networks, and metabolic networks among different gene groups (Carlson et al., 2006). To find the most effective genes in profiles 66 and 13, the gene co-expression network was constructed using the algorithmic prediction, and topological analysis was carried out, the results of which demonstrated that more than 20 effective genes were identified. Most of these genes were attributed to cell communication, which may be related to cell division or proliferation through cell adhesion and the transport of nutritients. Among these genes, CN161469 was found to be similar to methylcrotonoyl-Coenzyme A carboxylase 2 (beta), which suggests that it possesses the similar biological functions. In addition, the gene U89949.1, which encodes a folate binding protein that is involved in the transport of foliate (Vallet et al., 2001; Kim and Vallet, 2004), played an major role in the network. This finding was in accordance with the results of the assessment of the GO categories and pathway analysis. It has been demonstrated that down-regulated genes reflect a state of malnutrition, which may be induced by a lack of the essential molecules for growth, such as folic acid, trace elements, and aminoglycan, in early life. In profile 13, the gene U89949.1 was down-regulated, which may imply a decrease in folic acid transport that may subsequently induce malnutrition. This gene therefore provides a channel that allows regulation of the metabolism on a molecular level to balance malnutrition and over–nutrition.

IMPLICATIONS

In order to increase our understanding of the molecular mechanisms involved in intestinal development, it is essential to obtain more information on the functions of genes in the mammalian genome and their corresponding products. In the present study, we have provided a rich new information resource, which describes the significant pathways and effective transcripts that are related intimately with porcine intestinal development. Although there was a certain degree of ambiguity when we established the regulatory network based on the expression profiling data, this represents a valuable attempt to guide further in-depth research into the candidate genes and signal transduction pathways involved in intestinal development.