A Systematic Analysis on mRNA and MicroRNA Expression in Runting and Stunting Chickens

Runting and stunting syndrome (RSS), which is characterized by lower body weight, widely occurs in broilers. Some RSS chickens simply exhibit slow growth without pathological changes. An increasing number of studies indicate that broiler strains differ in susceptibility to infectious diseases, most likely due to their genetic differences. The objective of this study was to detect the differentially expressed miRNAs and mRNAs in RSS and normal chickens. By integrating miRNA with mRNA expression profiling, potential molecular mechanisms involved in RSS could be further explored. Twenty-two known miRNAs and 1,159 genes were differentially expressed in RSS chickens compared with normal chickens (P < 0.05). qPCR validation results displayed similar patterns. The differentially expressed genes were primarily involved in energy metabolism pathways. The antisense transcripts were extensively expressed in chicken liver albeit with reduced abundance. Dual-luciferase reporter assay indicated that gga-miR-30b/c directly target CARS through binding to its 3′UTR. The miR-30b/c: CARS regulation mainly occurred in liver. In thigh muscle and the hypothalamus, miR-30b/c are expressed at higher levels in RSS chickens compared with normal chickens from 2 to 6 w of age, and notably significant differences are observed at 4 w of age.


Introduction
Runting and stunting syndrome (RSS), which is characterized by low body weight, mainly occurs in broilers, and the condition is incredibly complicated [1][2][3]. Previous reports indicate that both genetic and environmental factors, such as diseases and feeding, are implicated in arrested development [4,5]. The body weights of RSS chickens are significantly reduced compared with normal chickens. RSS generally occurs early in life and always results in considerable economic losses, particularly in the commercial broiler industry. Previous reports indicated that RSS chickens are frequently associated with viral infections, such as avian leukosis virus [6], reticulo endotheliosis virus [7], digestive tract disease [8], and astroviruses [3]. Currently, no effective commercial vaccine available for the control of this disease exists, primarily owing to the absence of known etiologic agents. In addition, some RSS chickens grow slowly without pathological changes. These chickens typically survive until sold. An increasing number of studies indicate that broiler strains differ in susceptibility to infectious diseases, most likely due to their genetic differences [9][10][11].
RNA sequencing (RNA-seq) technology holds great potential for disease diagnosis and genetic differences analyses [12]. So far, large-scale mRNA expression analysis in RSS chickens is lacking. Additionally, an increasing number of reports indicate that microRNAs (miRNAs), such as let-7b, miR-1, miR-133 and miR-206 [13,14], function in animal growth by regulating their target genes. Therefore, an integrative analysis of miRNA and mRNA interaction in RSS chickens was conducted to investigate the genetic differences between RSS and normal chickens and identify potential markers for growth.
In the present study, we identified differentially expressed miRNAs and mRNAs. These genes were mainly enriched in energy metabolism pathways, and these genes could serve as candidate genes for RSS chickens. The data from this study provide a combined analysis of miRNAs and mRNAs and shed new light on the fine-tuning process of growth.

RSS chicken classification
In the experiment, we dissected a total of 400 RSS chickens and 24 (6%) of them exhibited no pathological changes except for lower body weight. The remaining chickens showed various symptoms, such as intestinal bleeding, yolk malabsorption, perihepatitis and proventriculitis (Figs 1 and 2). All the 24 RSS individuals were ALV and REV virus negative. Five negative RSS chicken and 5 normal healthy individuals were used for Solexa and Digital Gene Expression (DGE) sequencing.

Differentially expressed miRNAs between the livers of RSS and normal chickens
From Solexa miRNA sequencing, we obtained a total of 9,246,256 and 8,714,768 small RNA high quality reads in RSS and normal chicken livers, respectively. Length distribution analysis showed that most reads ranged from 21-23 nt with the percentage of the three types of reads being 77.18% and 79.10% for the two libraries, respectively (Fig 3). The high quality reads were subsequently annotated to different classes of RNA categories using different databases such as miRBase (V19.0), from Rfam (10.1) and Genbank to identify miRNAs, repeats-associated RNA, rRNA, tRNA, snRNA, snoRNA, etc. In total, 6,246,916 and 6,002,911 reads were annotated and the most abundant RNA species in the two libraries was classified as miRNAs (Table 1).
In this study, 202 and 194 known miRNAs were detected in RSS and normal chickens, respectively. Among these known miRNAs, 22 were differentially expressed (P < 0.05 and Fold-change !1) ( Table 2). Four up-regulated miRNAs (gga-miR-221, gga-miR-30b, gga-miR-30c and gga-miR-215) and 1 down-regulated miRNA (gga-miR-375) in RSS chicken were chosen randomly to validate the Solexa sequencing results using RT-qPCR. Our results indicated that the miRNA expressions displayed patterns similar to those observed using Solexa analysis (Table 3).

Differentially expressed genes between RSS and normal chicken livers
DGE-sequence obtained 3,768,716 and 3,632,492 high quality clean tags in RSS and normal chicken livers, respectively. In the obtained clean tags, 2,134,676 (56.64%) and 1,900,218 (52.31%) tags were mapped to 11,627 and 11,381chicken genes. Compared with normal chickens, a total of 1,159 genes (605 up-regulated and 554 down-regulated) were differentially expressed (P < 0.05 and Fold-change ! 1) in RSS chickens.
To verify the DGE sequencing results, two up-regulated genes (CKAP2L and TBPL1) and six down-regulated genes (MAP2K1, KDR, CARS, MAP3K7IP2, MAP3K5 and LARP7) were selected for qPCR analysis. The qPCR results displayed patterns similar to DGE sequencing (Table 4). Of the differentially expressed sequences, six genes (LOC415662, C1orf78, ZZEF1,   which has important functions in animal immunity [15]. Both genes exhibited significantly different expression in RSS versus normal chickens (P < 0.001).

Antisense transcription
In the present study, 7,237 genes in RSS and 7,363 genes in normal chickens were found to possess antisense transcripts, thereby indicating that the antisense transcripts were extensively expressed in chicken liver. The majority of the antisense tags were expressed in lower abundance than their corresponding sense tags. KEGG analysis of the genes with antisense tags implied that these genes were notably enriched in 15 pathways in RSS chickens and 17 pathways in  normal chickens (P < 0.05, Table 5). The amino acyl-tRNA biosynthesis pathway was particularly enriched in RSS chickens, whereas DNA replication, RNA degradation and ubiquitin-mediated proteolysis pathways were enriched in normal chickens (Table 5).
2.6. gga-miR-30b/c directly target CARS through binding to its 3 0 UTR Previous studies reported that gga-miR-30b/c were involved with growth and development [18][19][20][21][22]. The inhibitory relationship between gga-miR-30b/c and their predicted target gene CARS was verified using luciferase reporter assays. The wild-type and mutated target sequences of CARS 3 0 UTR were introduced into pmirGLO vector (Promega Biotech Co., Ltd., China), We co-transfected 200 ng pmirGLO-CARS or pmirGLO-CARS mutant, with 600 ng pCDNA3.1-miR-30b/c plasmid, into DF-1 cells. Thirty-six hours later, the cells were lysed, and the luciferase activity was measured using the Dual Luciferase Assay according to the manufacturer's protocol (Promega Biotech Co., Ltd., China). Compared with the cells that were transfected with pmirGLO-CARS mutant, the luciferase activity of CARS 3 0 UTR wild-type was significantly decreased. Also, the luciferase activity of the mock transfection was higher than that of the wild-type group (Fig 5). These results indicated that CARS 3 0 UTRwas directly targeted by miR-30b/c.

Temporal and spatial expression of miR-30b/c and CARS gene in RSS chickens
We firstly compared the miR-30b/c and CARS expression between normal and RSS chicken livers at three time points (Fig 6). qPCR results indicated that miR-30b/c were expressed in a higher level in RSS chicken than in the normal chickens, while CARS showed in an inversed trend. The differences of miR-30c and CARS between the two types of chickens reached a significant level at 5 or 6 w of age (Fig 6-B and 6-C) (P < 0.05). P value of miR-30b at 6 w of age was 0.08, very close to the statistically significant level (Fig 6-A). Secondly, we detected the miR-30b/c and CARS expression in pectoral muscle, thigh muscle, liver, hypothalamus and pituitary of chicken at 7 w of age and found that miR-30b in thigh muscles and CARS in hypothalami were differentially expressed between RSS and normal chickens (Fig 7). Hypothalamic-pituitary system can regulate hormone secretion and further affect animal growth. In addition, thigh muscle weight was an important economic trait in livestock production. Subsequently, qPCR was conducted to compare the expression of miR-30b/c and CARS in thigh muscle and hypothalami between RSS and normal chickens at various developmental stages (1 d, 2w, 4wand 6w of age). From 2 to 6 w, miR-30b/c expression was increased in the thigh muscles and hypothalami of RSS chickens compared with normal chickens, and it is notable that significant differences (P < 0.05) appeared at 4 w of age (Fig 8).
CARS expression varied inversely with miR-30b/c expression at 4 w and 6 w of age, although the differences did not attain to statistical significance (P > 0.05) (Fig 9).

Variation analysis of 3'UTR of CARS gene
We hypothesized that CARS gene variation most likely evened the expression differences between RSS and normal chickens. The 3'UTR of the CARS gene was amplified to detect variations in RSS or normal chicken DNA samples. However, the nucleotide sequence alignments indicated no nucleotide mutations in the amplified fragment.

Western blot analysis CARS protein in thigh muscle
Both in normal and in RSS chicken, CARS protein showed an increasing trend along with chicken development. But no obvious differences were observed in thigh muscles as compared normal with RSS chickens at the same stage because of individual differences (Fig 10).    Comparison of miR-30b/30c expression in thigh muscles and hypothalami of RSS and normal chickens at day 1 (D1), 2, 4 and 6 w of age. Quantitative RT-PCR assays of miR-30b and miR-30c expression using total RNA isolated from the chicken thigh muscle and hypothalamus at day1 (D1), 2, 4 and 6 w of age. "*" represents miR-30b or c was differentially expressed between RSS and normal chickens (P < 0.5).

Discussion
An increasing number of studies indicate that broiler strains differ in susceptibility to infectious diseases, most likely due to their genetic differences [9][10][11]. In this study, we found 6% RSS chickens exhibiting no pathological changes except for lower body weight. And further we observed significant differences in mRNA and miRNA profiles between RSS and normal chickens. qPCR results of the randomly selected miRNAs and mRNAs displayed good consistency with deep sequencing. In total, 1,159 genes were differentially expressed according to the DGE sequence. Of particular note, six genes (LOC415662, C1orf78, ZZEF1, CEBPA, LOC418836 and FOXP2) were specifically expressed in RSS chickens, whereas 5 genes (CINP, MYH3, EEF1A2, MAP3K13 and IMMP1L) were exclusively detected in the livers of normal chickens. MAP3K7IP2, also known as TAB2, is a member of the IL-1 signal transduction pathway, which is a central regulator of immune and inflammatory responses. Moreover, TAB2 is also involved in NF-kappaB activation in the TCR signaling pathway, which plays a key role in the adaptive immune response [23,24]. Previous cDNA microarray results indicated that immune-related genes were differentially expressed in the intestines of malabsorption syndrome broilers compared with control [10]. Malabsorption syndrome broilers exhibited growth retardation and weight gain depression, and these symptoms are similar to those observed in RSS chickens. Additionally, intestinal cell-mediated immunity also plays a crucial role in broiler development [25]. In short, immunity plays important roles in chicken growth and could be involved in RSS in chickens.
Pathway analysis revealed that the aforementioned genes were involved in 7 pathways at a significant level. It was notable that the enriched pathways were mainly related to energy metabolism, including oxidative phosphorylation, pyruvate metabolism, fatty acid metabolism and primary bile acid biosynthesis pathways. Previous comparative analyses indicated that meat-type chickens, which display a very rapid growth rate, have a higher oxidative phosphorylation efficiency compared with laying-type chickens [26]. In the current study, most genes in the oxidative phosphorylation pathway, such as NADH dehydrogenase (ubiquinone) 1 alpha subcomplex family genes, are expressed at lower levels in RSS chickens compared with normal chickens. Moreover, chickens displaying pale bird syndrome symptoms exhibit a substantial reduction in utilized dietary energy and absorbed lipids compared with normal birds from the commercial farm [27], thereby suggesting that energy metabolism is important in chicken growth.
Recently, antisense transcripts, especially long noncoding RNAs, have attracted increasing attention. Natural antisense transcripts have been identified in diverse organisms, such as human, mouse, chicken andzebrafish [28]. In the current study, the antisense transcripts were widely expressed in chicken liver. The previous known antisense transcripts of bFGF, IGF-II and programmed cell death 2 (Pdcd2) genes were also detected in our analysis, indicating that our sequencing results were reliable. Four predicted pathways of the genes with antisense transcripts were differentially compared in RSS and normal chickens. The amino acyl-tRNA biosynthesis pathway, which is specifically enriched in RSS chickens, regulates gene expression and amino acid biosynthesis [29]. The other three pathways (DNA replication, RNA degradation and ubiquitin-mediated proteolysis) specifically observed in normal chickens participate in numerous functions. The four pathways are important in animal development, suggesting that antisense transcripts could function in RSS chicken development. Notably, antisense tags of suppressor of cytokine signaling 3 (SOCS3); myosin, heavy chain 9, non-muscle (MYH9); and inhibitor of growth family member 1 (ING1) genes were also detected with lower abundance. These genes play key roles in the animal growth axis [13]. The role of these antisense tags in regulating chicken growth should be further experimentally verified and studied.
Identification of the complicated regulation mechanisms among natural antisense transcripts, miRNAs and mRNAs in animal growth and human health is a great challenge. The latest results on antisense transcription provide an important clue; through an RNA-RNA pairing interaction, the PTENpg antisense transcript affects PTEN protein output by altering PTENpg stability and miRNA sponge activity [30].
Twenty-two differentially expressed miRNAs from liver were identified from Solexa sequencing results. The integration of global profiling of miRNAs and mRNA expression potentially provides a unique opportunity in understanding animal growth [31,32]. Gga-miR-30b and gga-miR-30c share the same seed regions and belong to the gga-miR-30 family. Members of the miR-30 family function in osteoblast differentiation, adipogenic differentiation, angiogenesis, myocardial matrix remodeling and cell apoptosis by targeting various pathways [18][19][20][21][22]33]. In this study, gga-miR-30b/c expression was increased in liver in RSS chickens compared with normal chickens. The CARS gene was proven as the target of miR-30b/c. gga-miR-30b/c and CARS in liver showed inversed expression trends, especially at age of 6w (P < 0.05). Also, gga-miR-30b/c expression was increased in the thigh muscle and hypothalamus in RSS chickens compared with normal chickens at the indicated time points. The differences achieved statistical significance at 4 w of age when the maximum gap of body weight was observed between RSS and control chickens [10]. But in the two tissues, the CARS mRNA and protein expression did not demonstrate marked differences. In short, abnormal expression of gga-miR-30b/c and CARS might lead to RSS chickens by affecting liver metabolism and function.
In addition, variation analysis indicates that no mutations were present in the 3'UTR of this gene. Recent studies showed that the post-transcription regulation mechanism was far more complicated than we expected, and increased dark matters, such as long noncoding RNA, have been confirmed to function in growth and health [34]. Hence, further studies are needed to clarify the role of miR-30b/c in regulating muscle development.

Ethics statement
The Animal Care Committee of South China Agricultural University (Guangzhou, People's Republic of China) approved this study (approval number SCAU#0017). Animals involved in this study were humanely sacrificed as necessary to ameliorate suffering.

Sample collection
Four hundred female RSS yellow-feathered broiler chickens were dissected at 52 d of age. The ones, exhibiting no obvious pathological changes except for lower body weight, were selected for ALV and REV viruses detection. Five negative RSS chickens and 5 normal healthy individuals were used for Solexa and DGE sequencing.

Liver DNA and RNA extraction
Genomic DNA samples were used for ALV and REV virus detection and isolated from livers by a phenolic extraction protocol. Total RNA was extracted from 5 RSS and 5 normal chicken livers using Trizol (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. Total RNA was tested via agarose electrophoresis, Qubit-fluorometer and Agilent 2100. The five RNA samples isolated from each type of chicken were pooled with equal amounts from individual chickens.

Solexa sequencing of miRNA and Q-RT-PCR validation
The small RNAs (sRNA) were enriched from chicken livers to construct libraries for Solexa sequencing according to previous report [35]. Briefly, 18-40 base pair long RNA fragments were isolated from total RNA using a Novex 15% TBE-Urea gel (Invitrogen). Then, a 5' and 3' adaptor (Illumina, San Diego, CA, USA) were ligated to purified small RNAs followed by purification of ligation products on a Novex 15% TBE-Urea gel. These ligation products were subsequently reverse transcribed and PCR amplified. The purified DNA fragments were used for sequencing using an Illumina Genome Analyzer at the Beijing Genomics Institute, Shenzhen, China. After removing the adaptor/acceptor sequences, filtering the low quality tags and cleaning up the contamination formed by the adaptor-adaptor ligation, the occurrence of each unique sequence reads was counted as sequence tags. We compared the miRNA expression between two kinds of sample and found out the differentially expressed miRNAs according to previous report [36,37]. Briefly, we normalized the expression of miRNA in two samples (RSS and normal chicken) to get the expression of transcript per million (TPM-Normal and TPM-RSS) at first. Secondly, we calculated fold-change and P-value from the normalized expression. According to previous research, denote the number of unambiguous clean tag from gene A as x, as every gene's expression occupies only a small part of the library, the p(x) is in the Poisson distribution [36,37].
(λ represents the population mean tag number of x; The capital X represents random variables for the read counts.) The probability of gene A expressed equally between two samples can be calculated with: (The total clean tag number of the sample 1 is N1, and total clean tag number of sample 2 is N2; gene A holds x tags in sample1 and y tags in sample2; The x and y are not normalized. The capital Y represents random variables for the read counts) To compute the confidence intervals, we made use of the cumulative distributions: Cðy y min jxÞ ¼ X y y min y¼0 pðyjxÞandD ðy ! y max jxÞ ¼ X 1 y!y max pðyjxÞ which allow the computation of an interval [y min , y max ] ε and serve as a significance test when comparing [36,37]. The intervals are given for the 95% (2ε = 0.05) and 99% (2ε = 0.01) confidence levels and the P-value can be calculated as follows when the transcript count y and x of one gene between the two samples were compared:Without loss of geneality, assume that y > x FDR (False Discovery Rate) is the method to determine the threshold of P value in multiple test and analysis through manipulating the FDR value [38]. Assume that we have picked out R differentially expressed genes in which S genes are really show differential expression and the other V genes are false positive. We preseted the FDR to a number no larger than 0.001 according to previous report [38] because we decide that the error ratio "Q = V/R" must stay below a cutoff (e.g. 0.1%). For sequencing data, fold change was calculated according to the following formula, Fold change ¼jlog2Ratioj¼ j log2ðTPM À Normal= TPM À RSSÞ j One miRNA will be considered as differentially expressed between two types of samples when P < 0.05 and fold change !1. The differentially expressed miRNAs were selected to verify the results of Solexa sequencing using a miScript Reverse Transcription Kit and miScript SYBR GreenPCR Kit (Qiagen).

Digital gene expression sequencing of mRNA and Q-RT-PCR validation
DGE sequencing of mRNA was documented previously [28]. Briefly, cDNA was digested with NIaIII to obtain a cDNA fragment from the most 3' CATG to the poly(A)-tail. Subsequently, GEX adapter 1 was ligated to the free 5' end of the RNA, and the construct was digested with MmeI, which cuts 17 bp downstream of the CATG site. PCR was performed with primers to enrich the desired fragments. A preprocessed database of all possible CATG+17nt tag sequences was formed using the chicken genome and transcriptome. One gene will be considered as differentially expressed between two types of samples when P < 0.05 and fold change !1. The false discovery rate (FDR) was analyzed according to the method reported by Benjamini Yekutieli [38]. P value and fold change were calculated with formulas listed above. Eight differentially expressed genes were selected to confirm the results of DGE sequencing by Q-RT-PCR (ReverTraAce qPCR RT Kit).

Dual-luciferase activity assay
The wild-type and mutated target sequences of CARS 3 0 UTR predicted as target of gga-miR-30b/c was amplified by PCR from chicken cDNA and inserted into the downstream of luciferase reporter gene in the pmirGLO vector (Promega Biotech Co., Ltd., China) (pmirGLO-CARS). The pri-miR-30b/c and the flanking sequences, obtained by PCR from chicken DNA, were inserted into pCDNA3.1+ expression vector (Invitrogen) and subsequently were described as pCDNA3.1-miR-30b/c. pmirGLO-CARS and pCDNA3.1-miR-30b/c plasmids were cotransfected into chick DF-1 cell lines with lipofectamine 3000 reagent (Invitrogen). The relative luciferase activity was detected 36 hours later according to the manufacturing protocol of dual luciferase reporter gene assay kit (Promega Biotech Co., Ltd., China).

Western blot
Leg muscle tissues were collected in lysis buffer (50 mM Tris-HCl, pH 7.5, 150 mMNaCl, 1% Triton X-100, and 0.25% sodium deoxycholate) and were sonicated to shear DNA and reduce sample viscosity. Protein concentration was measured by BCA Protein Assay Kit (Thermo Scientific Pierce, USA). Samples were run on a 10% PAGE gel and transferred onto polyvinylidene fluoride membranes. After blocked in 5% nonfat dry milk in TPBS (0.1% Tween 20 in PBS) for 1 h, membranes were incubated with monoclonal anti-CARS antibody (ab83256, Abcam, SA, USA) overnight at 4°C. After five washes in 5% milk/TPBS 5 min each, membranes were incubated in goat anti-rabbit IgG conjugated with horseradish perioxidase for 1 h followed by two washes in 5% nonfat milk in TPBS, TPBS and PBS 5 min each, respectively. The signals were developed in BeyoECL Plus (Beyotime Institude of Biotech, China). The protein levels were normalized to GAPDH.

Statistical analysis
T-test was used to make miR-30b/c and CARS temporal and spatial expression comparisons between RSS and normal chickens.
Supporting Information S1