# Whole Blood Transcriptome Sequencing Reveals Gene Expression Differences between Dapulian and Landrace Piglets

Authors: Jiaqing Hu; Dandan Yang; Wei Chen; Chuanhao Li; Yandong Wang; Yongqing Zeng; Hui Wang
Journal: BioMed Research International
 (2016)

---

## Abstract

There is little genomic information regarding gene expression differences at the whole blood transcriptome level of different pig breeds at the neonatal stage. To solve this, we characterized differentially expressed genes (DEGs) in the whole blood of Dapulian (DPL) and Landrace piglets using RNA-seq (RNA-sequencing) technology. In this study, 83 DEGs were identified between the two breeds. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses identified immune response and metabolism as the most commonly enriched terms and pathways in the DEGs. Genes related to immunity and lipid metabolism were more highly expressed in the DPL piglets, while genes related to body growth were more highly expressed in the Landrace piglets. Additionally, the DPL piglets had twofold more single nucleotide polymorphisms (SNPs) and alternative splicing (AS) than the Landrace piglets. These results expand our knowledge of the genes transcribed in the piglet whole blood of two breeds and provide a basis for future research of the molecular mechanisms underlying the piglet differences.


---

## Body

## 1. Introduction

The domestic pig is a major livestock animal in China, where the long-term breeding and geographical separation have generated over 80 indigenous pig breeds, many of which have particular features [1]. Among them, Dapulian (DPL) pig is an indigenous breed, which is mainly distributed in Jining, Shandong Province, China. It has many desirable idioplasmatic traits such as high prolificacy, strong disease resistance, good meat flavor, and good adaptability to crude feed [1, 2]. However, its growth rate is slow and fat content is high. Landrace pigs have a better growth rate and higher lean content than DPL pigs but are weaker in disease resistance compared to DPL [1, 2]. Pork producers have great economic interest to improve breeding stocks with better growth performance and immune capacity, thus to meet the increasing consumer demand for higher quality pork and lower production cost. However, the genetic basis of such differences is still not well understood. Therefore, it will be very valuable to unravel the genetic mechanisms underlying the differences between these two breeds.Growth performance and immune capacity play essential roles in the economics of pig production. However, understanding the gene expression variation during body growth and immune capacity development at young ages is in its primitive stage. Blood cells constitute the first line of the immune defense system and pervade the whole body [3]. Moreover, peripheral blood cells share more than 80% of the transcriptome with nine tissues including brain, colon, heart, kidney, liver, lung, prostate, spleen, and stomach [4]. Thus, immune status and metabolic changes in other tissues may manifest through whole blood gene expression. Therefore, blood is now widely used as a surrogate tissue to identify the extent and kinetics for various economically important traits [4].The advent of next-generation sequencing (NGS) technology and progresses in bioinformatics have provided a useful platform to easily explore DEGs [5]. Moreover, NGS technology can simultaneously assess mRNA transcription patterns for all the genes in various species. RNA-seq using NGS technology has become widely used for high-throughput researches of gene expression. RNA-seq has been used to identify genes associated with economically important quantitative traits, such as body growth, immune, disease resistance, and meat quality [6–9]. Although RNA-seq technology has been applied in many species, including pigs [10, 11], limited information is available regarding its application in the transcriptome of newborn piglets. SNPs can be found inside candidate genes for artificial or natural selection and might provide a cost-effective way to obtain more information than other molecular markers. AS is an important mechanism to increase the genomic and proteomic diversity, which occurs widely in the process of gene transcription in eukaryotes. Up to date, many new AS events have been discovered in humans [12],Arabidopsis [13], andDrosophila melanogaster [14] and other organisms. There are six common types of AS events: Skipped Exon (SE), Intron Retention (IR), Alternative 3′ Splice Sites (A3SS), Alternative 5′ Splice Sites (A5SS), Alternative First Exon (AFE), and Alternative Last Exon (ALE) [15].To the best of our knowledge, few studies have been performed to investigate blood transcriptome for divergence in gene expression in newborn piglets [16, 17]. The growth performance and immune capacity are the main factors that affect pig herd productivity [18]. Newborn piglets of DPL and Landrace were used to discover DEGs, SNP, and ASs between two breeds, which may help to explain the growth/production and immune differences of the two breeds. It will advance our understanding of the mechanism of piglet breed differences and may be helpful to improve stock breeding with better growth performance and immune capacity.


## 2. Materials and Methods

### 2.1. Experimental Design and Ethical Statement

Three female piglets from each purebred DPL and Landrace were used in this study. After birth, the piglets were weaned at once and sent to a common nursery with the ambient temperature set at25±1°C. For one day, cow milk was supplied for these piglets to meet the nutrient needs. 1 mL whole blood samples were collected from the external jugular vein of one-day-old piglets and then immediately mixed with 3 mL RNALock Reagent (Tiangen Biotech Co., Beijing, China). The mixture was kept at room temperature for 2 h and then stored at −80°C until use. All animal experiments in this study were approved by the Animal Care and Use Committee of Shandong Agricultural University (Approval number: 2004006).


### 2.2. Total RNA Extraction and RNA Quality Analysis

Total RNA was extracted using a RNAsimple Total RNA kit (Tiangen Biotech Co., Beijing, China) according to the manufacturer’s instructions. RNA concentration and purity were assessed using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). RNA integrity was assessed using a RNA Nano 6000 Assay Kit on an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).


### 2.3. cDNA Library Construction and Next-Generation Sequencing (RNA-seq)

Poly(A) mRNAs were further purified from the extracted total RNA using magnetic oligo (dT) beads (Invitrogen, Carlsbad, CA, USA). Sequencing libraries were constructed using the NEBNext® Ultra RNA Library Prep Kit from Illumina (New England Biolabs, Ipswich, MA, USA) with multiplexing primers, following the manufacturer’s instructions. The average insert size for the paired-end libraries was 200 bp (150~250 bp). AMPure XP Beads (Beckman Coulter Inc., Brea, CA, USA) were used for cDNA purification. The purified double-stranded cDNAs were subjected to end repair and adapter ligation. Agencourt AMPure XP beads (Beckman Coulter Inc.) were then used for the selection of suitable fragments. Finally, the cDNA libraries were generated through PCR enrichment. A 125 bp paired-end run was performed on the Illumina HiSeq2500 platform at the Biomarker Technologies Corporation (Beijing, China). The RNA-seq data was deposited into NCBI database, and its accession is SAMN04423129.


### 2.4. Bioinformatics Analysis

#### 2.4.1. Reads Trimming and Reference Genome Alignment

After discarding low quality reads (reads with adaptors, unknown nucleotides larger than 5%, or Q20 < 20%) by perl script, clean reads that were filtered from the raw reads were mapped to the reference genome ofSus scrofa (Sscrofa 10.2; ftp://ftp.ensembl.org/pub/release-75/fasta/sus_scrofa/) using the TopHat2 software. Gene expression levels were measured using the FPKM (fragments per kilobase of exon per million fragments mapped) values by the Cufflinks software [22]. Genes (FPKM > 0.1) that were measured commonly in all the six samples were assigned as cogenes.


#### 2.4.2. Gene Annotation and DEGs Analysis

The FPKM value ratio was used to estimate gene expression level between those samples. The false discovery rate (FDR) control method was used to identify thep value in multiple tests. Fold changes (log 2 ratio) were estimated according to the normalized gene expression level in each sample. The DEGs were analyzed using the R package DESeq [23] between the two breeds. The input data for gene differential expression are the readcount data obtained from the gene expression level analysis. The DESeq software internally performs an FPKM conversion on the count value. The genes with fold changes (log 2 ratio) ≥ 2 and FDR significance score < 0.01 were considered as DEGs. The p value of the original hypothesis test was corrected by the accepted Benjamini-Hochberg correction method. Finally, FDR was used as the key index for DEGs screening.The identified DEGs were used for subsequent GO and KEGG pathway analysis. To annotate the genes with GO terms, the Nr BLAST results were imported into the blast2 GO program. GO annotations for the genes were obtained by blast2GO. This analysis mapped all of the annotated genes to GO terms in the database and counted the number of genes associated with each term. Perl script was then used to plot GO functional classification for the unigenes with a GO term hit to view the distribution of gene functions. The obtained annotation was enriched and refined using the TopGo (R package) with the “elim” method and Kolmogorov–Smirnov test. The genes were classified into three categories including biological process, cellular component, and molecular function. Moreover, the enriched pathways of DEGs were analyzed using the KEGG database (http://www.genome.ad.jp/kegg/) using the right sided Fisher’s exact test.


#### 2.4.3. Identification of SNPs and AS

RNA-seq technology is an effective way to quickly and reliably discover SNPs and AS. Comparisons were made between the two pig breed transcriptomes and the reference genome ofSus scrofa. Single base mismatches between sequencing samples and the reference genomes were identified using the SAMtools software [24] to find potential SNPs in the gene region. The criteria for SNPs identification included nucleotide with reads mapping quality score ≥ 50; consensus quality score ≥ 20; and sequence depth ≥ 5. For the identification of AS, the Cufflinks method was used to predict AS events [22]. To detect AS events, the blasted results using TopHat were assembled using Cufflinks and compared with initial annotated results using Cuffcompare. The variable splicing types and corresponding expression levels were obtained using the ASprofile software [25].


### 2.5. Quantitative Real-Time PCR (qRT-PCR)

To validate the RNA-seq data, qRT-PCR analysis was conducted using the MX3000p qPCR system (Stratagene, La Jolla, CA, USA) with SYBR green as the fluorescent dye, according to the manufacturer’s instructions (Takara Biotechnology, Dalian, China). Six functionally important genes were selected from the identified DEGs for qRT-PCR validation. The primers for qRT-PCR were designed either using the online Prime3 program (http://primer3.ut.ee/) or according to the published references [19–21, 26]. All the primers used in this study were listed in Table 1. All qRT-PCR reactions were conducted in triplicate with corresponding negative controls. Using the B2M and PPIA gene for normalization [27], the expression fold differences of the selected DEGs were calculated using the 2-ΔΔCT method [28].Table 1
Primer sequences for qRT-PCR.
Gene Sequence 5′-3′ Size (bp) Tm Reference or gene accession number G0S2 F: gtcgccttacgtttggacttgc 156 58 [19] R: caggtaactccgctcaggtgc ATP4B F: cccctgcaggtggaatactt 160 58 NM_001001258.1 R: ggatcttgcacacgatgacc GSTO1 F: attacctcatctggccctgg 150 58 NM_214050.1 R: ctcgaaggtctctcggttca CST3 F: cgagtacaacaaagcgagca 193 60 NM_001044602.1 R: cagcgttttcttctgcaggt FECH F: gatcgcgtttaccagtgacc 180 60 NM_001170523.1 R:cgctcattggactggatgtg MX2 F: aggcaaccaagagggaaatc 131 60 NM_001097416.1 R: cacactgatatgcccgatga B2M F: ttcacaccgctccagtag 166 60 [20] R: ccagatacatagcagttcagg PPIA F: cacaaacggttcccagtttt 171 60 [21] R: tgtccacagtcagcaatggt


## 2.1. Experimental Design and Ethical Statement

Three female piglets from each purebred DPL and Landrace were used in this study. After birth, the piglets were weaned at once and sent to a common nursery with the ambient temperature set at25±1°C. For one day, cow milk was supplied for these piglets to meet the nutrient needs. 1 mL whole blood samples were collected from the external jugular vein of one-day-old piglets and then immediately mixed with 3 mL RNALock Reagent (Tiangen Biotech Co., Beijing, China). The mixture was kept at room temperature for 2 h and then stored at −80°C until use. All animal experiments in this study were approved by the Animal Care and Use Committee of Shandong Agricultural University (Approval number: 2004006).


## 2.2. Total RNA Extraction and RNA Quality Analysis

Total RNA was extracted using a RNAsimple Total RNA kit (Tiangen Biotech Co., Beijing, China) according to the manufacturer’s instructions. RNA concentration and purity were assessed using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). RNA integrity was assessed using a RNA Nano 6000 Assay Kit on an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).


## 2.3. cDNA Library Construction and Next-Generation Sequencing (RNA-seq)

Poly(A) mRNAs were further purified from the extracted total RNA using magnetic oligo (dT) beads (Invitrogen, Carlsbad, CA, USA). Sequencing libraries were constructed using the NEBNext® Ultra RNA Library Prep Kit from Illumina (New England Biolabs, Ipswich, MA, USA) with multiplexing primers, following the manufacturer’s instructions. The average insert size for the paired-end libraries was 200 bp (150~250 bp). AMPure XP Beads (Beckman Coulter Inc., Brea, CA, USA) were used for cDNA purification. The purified double-stranded cDNAs were subjected to end repair and adapter ligation. Agencourt AMPure XP beads (Beckman Coulter Inc.) were then used for the selection of suitable fragments. Finally, the cDNA libraries were generated through PCR enrichment. A 125 bp paired-end run was performed on the Illumina HiSeq2500 platform at the Biomarker Technologies Corporation (Beijing, China). The RNA-seq data was deposited into NCBI database, and its accession is SAMN04423129.


## 2.4. Bioinformatics Analysis

### 2.4.1. Reads Trimming and Reference Genome Alignment

After discarding low quality reads (reads with adaptors, unknown nucleotides larger than 5%, or Q20 < 20%) by perl script, clean reads that were filtered from the raw reads were mapped to the reference genome ofSus scrofa (Sscrofa 10.2; ftp://ftp.ensembl.org/pub/release-75/fasta/sus_scrofa/) using the TopHat2 software. Gene expression levels were measured using the FPKM (fragments per kilobase of exon per million fragments mapped) values by the Cufflinks software [22]. Genes (FPKM > 0.1) that were measured commonly in all the six samples were assigned as cogenes.


### 2.4.2. Gene Annotation and DEGs Analysis

The FPKM value ratio was used to estimate gene expression level between those samples. The false discovery rate (FDR) control method was used to identify thep value in multiple tests. Fold changes (log 2 ratio) were estimated according to the normalized gene expression level in each sample. The DEGs were analyzed using the R package DESeq [23] between the two breeds. The input data for gene differential expression are the readcount data obtained from the gene expression level analysis. The DESeq software internally performs an FPKM conversion on the count value. The genes with fold changes (log 2 ratio) ≥ 2 and FDR significance score < 0.01 were considered as DEGs. The p value of the original hypothesis test was corrected by the accepted Benjamini-Hochberg correction method. Finally, FDR was used as the key index for DEGs screening.The identified DEGs were used for subsequent GO and KEGG pathway analysis. To annotate the genes with GO terms, the Nr BLAST results were imported into the blast2 GO program. GO annotations for the genes were obtained by blast2GO. This analysis mapped all of the annotated genes to GO terms in the database and counted the number of genes associated with each term. Perl script was then used to plot GO functional classification for the unigenes with a GO term hit to view the distribution of gene functions. The obtained annotation was enriched and refined using the TopGo (R package) with the “elim” method and Kolmogorov–Smirnov test. The genes were classified into three categories including biological process, cellular component, and molecular function. Moreover, the enriched pathways of DEGs were analyzed using the KEGG database (http://www.genome.ad.jp/kegg/) using the right sided Fisher’s exact test.


### 2.4.3. Identification of SNPs and AS

RNA-seq technology is an effective way to quickly and reliably discover SNPs and AS. Comparisons were made between the two pig breed transcriptomes and the reference genome ofSus scrofa. Single base mismatches between sequencing samples and the reference genomes were identified using the SAMtools software [24] to find potential SNPs in the gene region. The criteria for SNPs identification included nucleotide with reads mapping quality score ≥ 50; consensus quality score ≥ 20; and sequence depth ≥ 5. For the identification of AS, the Cufflinks method was used to predict AS events [22]. To detect AS events, the blasted results using TopHat were assembled using Cufflinks and compared with initial annotated results using Cuffcompare. The variable splicing types and corresponding expression levels were obtained using the ASprofile software [25].


## 2.4.1. Reads Trimming and Reference Genome Alignment

After discarding low quality reads (reads with adaptors, unknown nucleotides larger than 5%, or Q20 < 20%) by perl script, clean reads that were filtered from the raw reads were mapped to the reference genome ofSus scrofa (Sscrofa 10.2; ftp://ftp.ensembl.org/pub/release-75/fasta/sus_scrofa/) using the TopHat2 software. Gene expression levels were measured using the FPKM (fragments per kilobase of exon per million fragments mapped) values by the Cufflinks software [22]. Genes (FPKM > 0.1) that were measured commonly in all the six samples were assigned as cogenes.


## 2.4.2. Gene Annotation and DEGs Analysis

The FPKM value ratio was used to estimate gene expression level between those samples. The false discovery rate (FDR) control method was used to identify thep value in multiple tests. Fold changes (log 2 ratio) were estimated according to the normalized gene expression level in each sample. The DEGs were analyzed using the R package DESeq [23] between the two breeds. The input data for gene differential expression are the readcount data obtained from the gene expression level analysis. The DESeq software internally performs an FPKM conversion on the count value. The genes with fold changes (log 2 ratio) ≥ 2 and FDR significance score < 0.01 were considered as DEGs. The p value of the original hypothesis test was corrected by the accepted Benjamini-Hochberg correction method. Finally, FDR was used as the key index for DEGs screening.The identified DEGs were used for subsequent GO and KEGG pathway analysis. To annotate the genes with GO terms, the Nr BLAST results were imported into the blast2 GO program. GO annotations for the genes were obtained by blast2GO. This analysis mapped all of the annotated genes to GO terms in the database and counted the number of genes associated with each term. Perl script was then used to plot GO functional classification for the unigenes with a GO term hit to view the distribution of gene functions. The obtained annotation was enriched and refined using the TopGo (R package) with the “elim” method and Kolmogorov–Smirnov test. The genes were classified into three categories including biological process, cellular component, and molecular function. Moreover, the enriched pathways of DEGs were analyzed using the KEGG database (http://www.genome.ad.jp/kegg/) using the right sided Fisher’s exact test.


## 2.4.3. Identification of SNPs and AS

RNA-seq technology is an effective way to quickly and reliably discover SNPs and AS. Comparisons were made between the two pig breed transcriptomes and the reference genome ofSus scrofa. Single base mismatches between sequencing samples and the reference genomes were identified using the SAMtools software [24] to find potential SNPs in the gene region. The criteria for SNPs identification included nucleotide with reads mapping quality score ≥ 50; consensus quality score ≥ 20; and sequence depth ≥ 5. For the identification of AS, the Cufflinks method was used to predict AS events [22]. To detect AS events, the blasted results using TopHat were assembled using Cufflinks and compared with initial annotated results using Cuffcompare. The variable splicing types and corresponding expression levels were obtained using the ASprofile software [25].


## 2.5. Quantitative Real-Time PCR (qRT-PCR)

To validate the RNA-seq data, qRT-PCR analysis was conducted using the MX3000p qPCR system (Stratagene, La Jolla, CA, USA) with SYBR green as the fluorescent dye, according to the manufacturer’s instructions (Takara Biotechnology, Dalian, China). Six functionally important genes were selected from the identified DEGs for qRT-PCR validation. The primers for qRT-PCR were designed either using the online Prime3 program (http://primer3.ut.ee/) or according to the published references [19–21, 26]. All the primers used in this study were listed in Table 1. All qRT-PCR reactions were conducted in triplicate with corresponding negative controls. Using the B2M and PPIA gene for normalization [27], the expression fold differences of the selected DEGs were calculated using the 2-ΔΔCT method [28].Table 1
Primer sequences for qRT-PCR.
Gene Sequence 5′-3′ Size (bp) Tm Reference or gene accession number G0S2 F: gtcgccttacgtttggacttgc 156 58 [19] R: caggtaactccgctcaggtgc ATP4B F: cccctgcaggtggaatactt 160 58 NM_001001258.1 R: ggatcttgcacacgatgacc GSTO1 F: attacctcatctggccctgg 150 58 NM_214050.1 R: ctcgaaggtctctcggttca CST3 F: cgagtacaacaaagcgagca 193 60 NM_001044602.1 R: cagcgttttcttctgcaggt FECH F: gatcgcgtttaccagtgacc 180 60 NM_001170523.1 R:cgctcattggactggatgtg MX2 F: aggcaaccaagagggaaatc 131 60 NM_001097416.1 R: cacactgatatgcccgatga B2M F: ttcacaccgctccagtag 166 60 [20] R: ccagatacatagcagttcagg PPIA F: cacaaacggttcccagtttt 171 60 [21] R: tgtccacagtcagcaatggt


## 3. Results

### 3.1. Overview of the Sequencing Data

Across all 6 samples, the total read length was 30.2 gigabases (Gb). The fold change of reading per animal compared to the pig genome (2.8 Gb) is 2.04-, 1.93-, 1.82-, 1.63-, 1.80-, and 1.57-fold, respectively. The main statistical data of the six samples were described in Table2. Sequencing read length of all the samples was 125 bp paired-end reads, with Q30 > 88%. More than two-thirds of the reads could be mapped to the reference genome, with mapped ratio from 67.7% to 79.41%. In addition, more than two-thirds of the mapped reads had unique genomic locations (except DPL 3). The variation between individual blood samples was evaluated through pairwise correlation. Except one sample (DPL 3), a strong correlation among other five samples was found (Spearman correlation coefficient across all genes ranging from 0.850 to 0.985); therefore, the DPL 3 sample was excluded from further analysis.Table 2
Summary of the RNA-sequencing and mapping.
Sample Landrace 1 Landrace 2 Landrace 3 DPL 1 DPL 2 DPL 3 Total reads 45768956 43112026 40759848 36558806 40206856 35205066 Total base pairs 5766888456 5432115276 5135740848 4606409556 5066063856 4435838316 Mapped reads 35504697 (77.57%) 34236396 (79.41%) 31893141 (78.25%) 26752.117 (73.18%) 31099885 (77.35%) 23833962 (67.70%) Unique mapped reads 34286547 (74.91%) 33172939 (76.95%) 30790776 (75.54%) 25714679 (70.34%) 29884275 (74.33%) 22,473,200 (63.84%) % ≥ Q30 89.68 89.65 89.12 88.97 89.67 89.00 GC content% 58.35 59.19 58.53 57.77 57.22 56.25


### 3.2. Identification and Analysis of DEGs

By comparing the transcriptome data of the two breeds, a total of 83 DEGs were identified, including 57 genes with significantly higher expression and 26 genes with significantly lower expression in the Landrace piglets compared to the DPL piglets. The major DEGs are listed in Supplementary Table S1 in Supplementary Material available online athttp://dx.doi.org/10.1155/2016/7907980 including their log fold difference and FDR values. A heatmap that depicted the majority of DEGs was presented in Figure 1.Figure 1
Heatmap of the DEGs between the DPL and Landrace piglets. Columns represent individual samples. Rows indicate genes with significant expression differences between the two breeds (DPL, left columns 1 and 2; Landrace, right columns 1, 2, and 3).


### 3.3. GO and KEGG Pathway Enrichment Analysis of DEGs

The GO and KEGG enrichment analysis were performed to gain insights into the biological implications of the identified DEGs. GO analysis was performed using the blast2GO software [29]. The GO analysis showed that the DEGs clustered in molecular functions, biological processes, and cellular components (Figure 2). The GO annotation indicated that the DEGs were involved in many biological processes, such as organ morphogenesis, skeletal system development, immune response, protein binding, and integral component of membrane.Figure 2
GO analysis of DEGs between the DPL and Landrace piglets. The DEGs are classified into three categories: cellular component, molecular function, and biological process. The percentage of genes in each category and the number of genes are shown above.The DEGs were also mapped into the KEGG pathway database to predict the significantly enriched pathways and the detailed information was shown in Supplementary S2. The major enriched pathways included genetic information processing, metabolism, human disease, cellular processes, environmental information processing, and organismal systems.


### 3.4. DEGs Related to Immunity, Growth, and Metabolism

Among the DEGs, we found many genes related to immunity, which may be possible reasons for the different immune capability between the two breeds. For instance, myxovirus resistance 2 (MX2), a novel innate immunity restriction factor with important roles in defensing virus infection, had higher expression in the DPL piglets than the Landrace piglets (1.56-fold). CEBPE, a member of the basic-leucine zipper transcription factor family, mediating innate immune response and systemic lipid metabolism, had higher expression in the DPL piglets compared to the Landrace piglets (2.27-fold). Other immune-related genes including GNLY, GZMA, CSF1R, ADAM8, and CD59 were all differentially expressed between the two breeds.We also detected some DEGs relevant with lipid metabolism and growth (such as GOS2 and RGS2) had higher expression in the DPL piglets compared to the Landrace piglets (6.0-fold and 1.58-fold, resp.). However, body growth relevant genes PHOSPHO1, KIT, and CST3 had higher expression in the Landrace piglets compared to the DPL piglets (2.32-fold, 3.92-fold, and 1.66-fold, resp.).


### 3.5. Statistical Analysis of SNPs

The average SNPs in the DPL piglets were 37344, twofold more than the Landrace piglets. The genic SNPs were higher in DPL piglets; however, the intergenic SNPs were higher in the Landrace piglets. The frequencies of transition, transversion, and heterozygosity were similar in the two breeds. The details of different SNPs in the two breeds were listed in Table3 and Supplementary S3.Table 3
The distribution of the SNPs in the two breeds.
Animal SNP number Genic SNP Intergenic SNP Transition Transversion Heterozygosity L1 16,943 13,104 3,839 73.77% 26.23% 54.84% L2 10,135 8,148 1,987 73.34% 26.66% 58.59% L3 13,754 11,011 2,743 73.68% 26.32% 57.93% D1 38,153 27,022 11,131 74.53% 25.47% 44.02% D2 36,535 27,461 9,074 75.19% 24.81% 58.32%


### 3.6. Identification of AS

AS is a mechanism that brings remarkable diversity to proteins and allows a gene to generate different mRNA transcripts that are translated into distinct proteins. Skipped Exon, Intron Retention, Alt3splice, and Alt First Exon were the major AS in the two breeds. There were more Skipped Exon and Alt First Exon in the DPL than the Landrace piglets. However, other AS were similar in the two breeds. The distribution of the AS events was shown in Table4.Table 4
The distribution of AS in the two breeds.
Animal Skipped Exon Intron Retention Alternative 5′ Splice Sites Alt 3′ Splice Sites Alt Last Exon Alt First Exon L1 1,782 477 270 575 208 678 L2 1,049 440 199 373 143 456 L3 1,476 483 221 508 166 560 D1 2,048 623 286 653 215 788 D2 2,519 598 359 780 266 1,068


### 3.7. qRT-PCR Confirmation of the Gene Expression Data from RNA-seq

To verify the gene expression data by RNA-seq, qRT-PCR was performed for six genes according to their biological function, including two genes with lower expression levels and four genes with higher expression levels in the Landrace piglets compared to the DPL piglets. The expression patterns of the six genes were generally consistent with the RNA-seq results (Figure3), suggesting that the results of the RNA-seq experiments were accurate and reliable.Figure 3
Validation of the 6 DEGs from RNA-seq analysis by qRT-PCR. The direction and magnitude of the fold changes obtained using the qRT-PCR technique were similar to those of the RNA-seq data (p∗<0.05; p∗∗<0.01).


## 3.1. Overview of the Sequencing Data

Across all 6 samples, the total read length was 30.2 gigabases (Gb). The fold change of reading per animal compared to the pig genome (2.8 Gb) is 2.04-, 1.93-, 1.82-, 1.63-, 1.80-, and 1.57-fold, respectively. The main statistical data of the six samples were described in Table2. Sequencing read length of all the samples was 125 bp paired-end reads, with Q30 > 88%. More than two-thirds of the reads could be mapped to the reference genome, with mapped ratio from 67.7% to 79.41%. In addition, more than two-thirds of the mapped reads had unique genomic locations (except DPL 3). The variation between individual blood samples was evaluated through pairwise correlation. Except one sample (DPL 3), a strong correlation among other five samples was found (Spearman correlation coefficient across all genes ranging from 0.850 to 0.985); therefore, the DPL 3 sample was excluded from further analysis.Table 2
Summary of the RNA-sequencing and mapping.
Sample Landrace 1 Landrace 2 Landrace 3 DPL 1 DPL 2 DPL 3 Total reads 45768956 43112026 40759848 36558806 40206856 35205066 Total base pairs 5766888456 5432115276 5135740848 4606409556 5066063856 4435838316 Mapped reads 35504697 (77.57%) 34236396 (79.41%) 31893141 (78.25%) 26752.117 (73.18%) 31099885 (77.35%) 23833962 (67.70%) Unique mapped reads 34286547 (74.91%) 33172939 (76.95%) 30790776 (75.54%) 25714679 (70.34%) 29884275 (74.33%) 22,473,200 (63.84%) % ≥ Q30 89.68 89.65 89.12 88.97 89.67 89.00 GC content% 58.35 59.19 58.53 57.77 57.22 56.25


## 3.2. Identification and Analysis of DEGs

By comparing the transcriptome data of the two breeds, a total of 83 DEGs were identified, including 57 genes with significantly higher expression and 26 genes with significantly lower expression in the Landrace piglets compared to the DPL piglets. The major DEGs are listed in Supplementary Table S1 in Supplementary Material available online athttp://dx.doi.org/10.1155/2016/7907980 including their log fold difference and FDR values. A heatmap that depicted the majority of DEGs was presented in Figure 1.Figure 1
Heatmap of the DEGs between the DPL and Landrace piglets. Columns represent individual samples. Rows indicate genes with significant expression differences between the two breeds (DPL, left columns 1 and 2; Landrace, right columns 1, 2, and 3).


## 3.3. GO and KEGG Pathway Enrichment Analysis of DEGs

The GO and KEGG enrichment analysis were performed to gain insights into the biological implications of the identified DEGs. GO analysis was performed using the blast2GO software [29]. The GO analysis showed that the DEGs clustered in molecular functions, biological processes, and cellular components (Figure 2). The GO annotation indicated that the DEGs were involved in many biological processes, such as organ morphogenesis, skeletal system development, immune response, protein binding, and integral component of membrane.Figure 2
GO analysis of DEGs between the DPL and Landrace piglets. The DEGs are classified into three categories: cellular component, molecular function, and biological process. The percentage of genes in each category and the number of genes are shown above.The DEGs were also mapped into the KEGG pathway database to predict the significantly enriched pathways and the detailed information was shown in Supplementary S2. The major enriched pathways included genetic information processing, metabolism, human disease, cellular processes, environmental information processing, and organismal systems.


## 3.4. DEGs Related to Immunity, Growth, and Metabolism

Among the DEGs, we found many genes related to immunity, which may be possible reasons for the different immune capability between the two breeds. For instance, myxovirus resistance 2 (MX2), a novel innate immunity restriction factor with important roles in defensing virus infection, had higher expression in the DPL piglets than the Landrace piglets (1.56-fold). CEBPE, a member of the basic-leucine zipper transcription factor family, mediating innate immune response and systemic lipid metabolism, had higher expression in the DPL piglets compared to the Landrace piglets (2.27-fold). Other immune-related genes including GNLY, GZMA, CSF1R, ADAM8, and CD59 were all differentially expressed between the two breeds.We also detected some DEGs relevant with lipid metabolism and growth (such as GOS2 and RGS2) had higher expression in the DPL piglets compared to the Landrace piglets (6.0-fold and 1.58-fold, resp.). However, body growth relevant genes PHOSPHO1, KIT, and CST3 had higher expression in the Landrace piglets compared to the DPL piglets (2.32-fold, 3.92-fold, and 1.66-fold, resp.).


## 3.5. Statistical Analysis of SNPs

The average SNPs in the DPL piglets were 37344, twofold more than the Landrace piglets. The genic SNPs were higher in DPL piglets; however, the intergenic SNPs were higher in the Landrace piglets. The frequencies of transition, transversion, and heterozygosity were similar in the two breeds. The details of different SNPs in the two breeds were listed in Table3 and Supplementary S3.Table 3
The distribution of the SNPs in the two breeds.
Animal SNP number Genic SNP Intergenic SNP Transition Transversion Heterozygosity L1 16,943 13,104 3,839 73.77% 26.23% 54.84% L2 10,135 8,148 1,987 73.34% 26.66% 58.59% L3 13,754 11,011 2,743 73.68% 26.32% 57.93% D1 38,153 27,022 11,131 74.53% 25.47% 44.02% D2 36,535 27,461 9,074 75.19% 24.81% 58.32%


## 3.6. Identification of AS

AS is a mechanism that brings remarkable diversity to proteins and allows a gene to generate different mRNA transcripts that are translated into distinct proteins. Skipped Exon, Intron Retention, Alt3splice, and Alt First Exon were the major AS in the two breeds. There were more Skipped Exon and Alt First Exon in the DPL than the Landrace piglets. However, other AS were similar in the two breeds. The distribution of the AS events was shown in Table4.Table 4
The distribution of AS in the two breeds.
Animal Skipped Exon Intron Retention Alternative 5′ Splice Sites Alt 3′ Splice Sites Alt Last Exon Alt First Exon L1 1,782 477 270 575 208 678 L2 1,049 440 199 373 143 456 L3 1,476 483 221 508 166 560 D1 2,048 623 286 653 215 788 D2 2,519 598 359 780 266 1,068


## 3.7. qRT-PCR Confirmation of the Gene Expression Data from RNA-seq

To verify the gene expression data by RNA-seq, qRT-PCR was performed for six genes according to their biological function, including two genes with lower expression levels and four genes with higher expression levels in the Landrace piglets compared to the DPL piglets. The expression patterns of the six genes were generally consistent with the RNA-seq results (Figure3), suggesting that the results of the RNA-seq experiments were accurate and reliable.Figure 3
Validation of the 6 DEGs from RNA-seq analysis by qRT-PCR. The direction and magnitude of the fold changes obtained using the qRT-PCR technique were similar to those of the RNA-seq data (p∗<0.05; p∗∗<0.01).


## 4. Discussion

In this study, we compared the gene expression differences of the two genetically different pig breeds, the DPL, and the Landrace pigs. To the best of our knowledge, our study was the first report on transcriptome gene expression profiling of whole blood in DPL and Landrace piglets using RNA-seq technology. The current study generated 30.2 Gb clean data, more than 10-fold the pig genome size. Moreover, the map ration to reference genome ranged from 77.35% to 79.4%, meeting the requirements for the subsequent analysis. In this study, we usedSus scrofa as reference sequence; the low mapping rates could be related to the difference between the sample genome and the reference genome. Moreover, compared with human and mouse genomic sequence, pig genomic sequence is incomplete. Therefore, the mapping rate of pigs is generally less than 80% [30]. The identified DEGs, SNPs, and AS may help to illustrate the underlying differences in the whole blood transcriptome between these two breeds and will be valuable for future studies on the identification of major genes that may affect pig growth performance, immune capacity, and meat quality.In this study, totally 83 genes were shown to be significantly differentially expressed between the two breeds. The GO annotation and KEGG pathways analysis showed that the DEGs were mainly involved in single-organism process, cellular process and metabolic process of biological process, the cell part, cell and organelle of a cellular component, and the binding and catalytic activity of molecular function. Our finding was in accordance with the transcriptome analysis of other studies [10, 31]. It is well known that the molecular regulation of animal traits is very complicated and the relationship between genes and traits can be one-to-many or many-to-one. Therefore, it was not surprising that our identified DEGs were enriched not only in immune-related pathways but also in those involved in metabolic process. The immune-related and metabolic pathways are important for improving pig traits such as disease resistance and growth. Based on our data, further functional studies with these DEGs are warranted to identify key genes influencing growth performance and immune capacity of swine.Our current study has identified several differentially expressed genes between DPL and Landrace that are important in disease resistance, such as MX2 and CEBPE. The MX2 gene was identified as a novel IFN-induced gene, which inhibited a variety of virus infection [32–34]. Studies have shown that the expression levels of MX2 was upregulated in piglets after coinfection of porcine reproductive and respiratory syndrome virus (PRRSV) andMycoplasma hyopneumoniae [35]; it was also increased with the sudden change of feeding regime after weaning in piglet gut [36]. Moreover, MX2 also showed antiviral ability against influenza virus [37], vesicular stomatitis virus [33], and PRRSV in swine [38]. Additionally, previous studies have shown that DPL piglets exhibited strong resistance to PRRSV [39, 40]. In the current study, the expression level of MX2 gene was higher in the DPL piglets than the Landrace piglets. Therefore, higher transcription level of MX2 in DPL could be an indicator for resistance to PRRSV. CEBPE, which expressed only in myeloid cells including monocytes and macrophages, has been shown to be required not only for hematopoietic cell development, but also for cytokine expression [41–43]. Studies have shown that CEBPE-deficient mice developed normally but failed to generate functional neutrophils [44]. Furthermore, the phagocytic function of CEBPE-deficient macrophages was also impaired [45]. The higher expression level of CEBPE in the DPL piglets may also suggest CEBPE as an indicator for disease resistance. In addition, other immune-related genes including GNLY, GZMA, CSF1R, ADAM8, and CD59 were also differentially expressed in the two breeds. The expression of the immune relevant DEGs in DPL and Landrace pigs could be related to the differences of the breed and the immune system. Our data suggest that these genes might play important roles in disease resistance after birth.Our study also identified DEGs between DPL and Landrace that are important for body growth, including KIT, PHOSPHO1, and CST3. Compared with other industrial pig breeds, Chinese indigenous pigs have the disadvantage of slow growth rate and high fat content, but with the advantage of high fertility, special meat flavor, and strong adaptability [1]. Our data showed that KIT gene had higher expression in the Landrace piglets compared to DPL. It was reported that the KIT genotypes affected piglet birth weight and body weight gain until weaning [46]. Two other important candidate genes for body growth, namely, PHOSPHO1 and CST3, had higher expression in the Landrace piglets compared to the DPL. PHOSPHO1 is a member of the haloacid dehalogenase superfamily and it plays an essential for the initiation of skeletal mineralization [47]. Recently, PHOSPHO1 was found to have an important function during growth and skeletal development in mouse, which could avoid low bone mineral density, spontaneous greenstick fractures, osteomalacia, and prominent scoliosis [47, 48]. Moreover, PHOSPHO1 was shown to play an important role in activating inside chondrocyte- and osteoblast-derived matrix vesicles [48]. CST3 serves an important physiological role as the controlling inhibitor of extracellular cysteine proteinase activity [49]. Previous studies also suggested it as an excellent candidate gene to influence fetal development and growth [50]. The current study has revealed highly enriched DEGs in piglets that are involved in organ morphogenesis, skeletal system development, and immune. These data provide a basis for future functional research that may aid in the discovery of genetic mechanism of body growth in pigs.Genes that are important in adiposity were also found to be differently expressed between the DPL and Landrace piglets, including G0S2 and CEBPE. Adiposity can increase the cost for pig farming because it reduces feed efficiency and carcass yield. Lipids and their metabolites can regulate gene expression and therefore affect many physiological processes. Our data showed that G0S2 and CEBPE were highly expressed in the DPL piglets compared to Landrace. G0/G1 switch gene 2 (G0S2) was initially found to be differentially expressed in lymphocytes during lectin-induced switch from the G0 to G1 phase of the cell cycle [51, 52]. Moreover, G0S2 was abundantly expressed in porcine adipose tissue and liver among various tissues [53]. Other data also showed that G0S2 was associated with fat content in bovine species [54]. Fatty acid composition is a critical aspect of meat quality, and its variation affects flavor, color, firmness, and softness of the fat in meat [55]. Since DPL pigs have higher fat content than Landrace, G0S2 may be a promising candidate gene for meat quality. CEBPE not only associates with immunity but also is involved in systemic lipid metabolism [56]. Macrophages from CEBPE-deficient mice showed less lipid accumulation than control mice [57]. At present, limited information is available about the roles of G0S2 and CEBPE in pigs; therefore, further investigations are needed to explain the role of these genes in porcine lipid metabolic. In addition, other genes unrelated to growth performance, immune capacity, or meat quality were found to be differently expressed in the two breeds. These genes included HRH4, REXO2, NTAN1, CA2, REXO2, and FAM213A. However, there are very few studies that have been performed to study these genes’ functions. Further studies are needed to gain more insights into the function of these genes.NGS has been widely used to analyze SNPs in eukaryotes [31, 58]. The DPL piglets have twofold more SNPs than the Landrace piglets. These different SNPs in the two breeds will be valuable molecular markers for future study. SNPs identification may be affected due to RNA editing, sequencing errors (in spite of sequencing quality filtration), mapping error, or the reference sequencing errors [59]. Therefore, further studies are needed to validate the identified SNPs.Another important advantage of RNA-seq is its ability to detect AS events [22]. AS is an important model of gene expression regulation and is not generally assessable using microarray method. We found that Skipped Exon was the most common type of AS in pigs, which is similar to human [24] and yeast [60]. Our data showed that DPL had more AS than Landrace, which is consistent with the SNPs result. The agreement of SNPs and AS result may suggest that many genes undergo the same transcriptional and translational regulation in both breeds. Further research will be able to determine the detailed regulation in each breed, and theses observed differences will be highly useful in future studies.In summary, this study was conducted to screen DEGs from the whole blood of DPL and Landrace piglets using RNA-seq technology. We found that 83 genes were significantly differentially expressed between the two breeds. Bioinformatics analyses identified significantly enriched pathways that are important for disease resistance and growth performance. The identification of specialized biological functions and key genes will be very valuable for further studies for breed improvement programs. Additionally, more SNPs and AS were found in the DPL than Landrace pigs. These findings provide new insights into the postnatal status of the piglets and will facilitate future genomic and gene function research on the different pig breeds.



---

*Source: 7907980-2016-12-26.xml*