ORIGINAL RESEARCH article

Front. Physiol., 20 October 2017

Sec. Avian Physiology

Volume 8 - 2017 | https://doi.org/10.3389/fphys.2017.00816

Identification and Differential Abundance of Mitochondrial Genome Encoding Small RNAs (mitosRNA) in Breast Muscles of Modern Broilers and Unselected Chicken Breed

  • 1. Department of Poultry Science, Center of Excellence for Poultry Science, University of Arkansas, Fayetteville, AR, United States

  • 2. School of Human Environmental Sciences, University of Arkansas, Fayetteville, AR, United States

Abstract

Background: Although small non-coding RNAs are mostly encoded by the nuclear genome, thousands of small non-coding RNAs encoded by the mitochondrial genome, termed as mitosRNAs were recently reported in human, mouse and trout. In this study, we first identified chicken mitosRNAs in breast muscle using small RNA sequencing method and the differential abundance was analyzed between modern pedigree male (PeM) broilers (characterized by rapid growth and large muscle mass) and the foundational Barred Plymouth Rock (BPR) chickens (characterized by slow growth and small muscle mass).

Methods: Small RNA sequencing was performed with total RNAs extracted from breast muscles of PeM and BPR (n = 6 per group) using the 1 × 50 bp single end read method of Illumina sequencing. Raw reads were processed by quality assessment, adapter trimming, and alignment to the chicken mitochondrial genome (GenBank Accession: X52392.1) using the NGen program. Further statistical analyses were performed using the JMP Genomics 8. Differentially expressed (DE) mitosRNAs between PeM and BPR were confirmed by quantitative PCR.

Results: Totals of 183,416 unique small RNA sequences were identified as potential chicken mitosRNAs. After stringent filtering processes, 117 mitosRNAs showing >100 raw read counts were abundantly produced from all 37 mitochondrial genes (except D-loop region) and the length of mitosRNAs ranged from 22 to 46 nucleotides. Of those, abundance of 44 mitosRNAs were significantly altered in breast muscles of PeM compared to those of BPR: all mitosRNAs were higher in PeM breast except those produced from 16S-rRNA gene. Possibly, the higher mitosRNAs abundance in PeM breast may be due to a higher mitochondrial content compared to BPR. Our data demonstrate that in addition to 37 known mitochondrial genes, the mitochondrial genome also encodes abundant mitosRNAs, that may play an important regulatory role in muscle growth via mitochondrial gene expression control.

Introduction

Production efficiency in animal agriculture is critically important in order to meet the increasing needs of high quality meat protein of a growing global human population. Modern broiler chickens have been highly selected for rapid growth and feed efficiency compared to unselected (progenitor or heritage) chicken breeds. Consequently, the understanding of key regulatory mechanisms for rapid muscle growth and enhanced feed efficiency is imperative with regard to more efficient, and therefore sustainable, food animal production (Niemann et al., 2011).

Global gene expression studies using cDNA microarray assay (Kong et al., 2011; Bottje et al., 2012; Bottje and Kong, 2013), RNA sequencing (Bottje et al., 2017), and shotgun proteomics (Kong et al., 2016) have been carried out on breast muscle obtained from pedigree males (PeM) exhibiting either a high or low feed efficiency (FE) phenotype. Additionally, global transcriptomics has been conducted on breast obtained from PeM broilers (characterized by rapid growth, large muscle mass and higher feed efficiency) compared with the foundational Barred Plymouth Rock (BPR) chicken breed (exhibiting slow growth, small muscle mass and lower feed efficiency) (Kong et al., 2017). The BPR chicken exhibits a characteristic pattern of alternating white and black bars of feather pigmentation and was developed in the United States during the mid-nineteenth century (Dorshorst and Ashwell, 2009). The BPR breed was developed for dual purposes of both meat and egg production and has a much slower growth rate compared to commercial broilers (Lopez et al., 2007). Several global expression studies showed that production efficiency may be associated with various cellular mechanisms including mitochondrial oxidative stress, myogenic growth and differentiation, inflammatory response, protein degradation, stress responses, growth hormone signaling, cell cycle, apoptosis, and fatty acid transportation. Particularly, the association between production efficiency and mitochondrial functions on energy expenditure and oxidative phosphorylation system (OXPHOS) in chicken breast muscle has been well-documented (Bottje et al., 2002, 2006, 2009; Bottje and Carstens, 2009; Bottje and Kong, 2013).

A major function of mitochondria is in producing energy (ATP) through OXPHOS in most eukaryotic cells (Taylor and Turnbull, 2007). Mitochondria contain their own mitochondrial DNA (mtDNA) that exhibits characteristics that differ from nuclear DNA, including divergent genetic codes, transmission by maternal inheritance of mitochondrial type, higher rates of mutation due to proximity of mtDNA to the electron transport chain, polyploidy status, and a more compact organization (Fernandez-Silva et al., 2003). Similar to mammalian mitochondrial genome, chicken mtDNA harbors one control region (D-loop) and contains 37 genes encoding 2 rRNA subunits (12S and 16S rRNAs), 22 transfer RNAs (tRNAs) and 13 proteins that are all involved in OXPHOS (Zhao et al., 2016).

In addition to mitochondrial RNAs, mtDNA encoding small RNAs (mitosRNAs) were recently discovered in human, mouse and rainbow trout (Ro et al., 2013; Ma et al., 2016). These novel mitosRNAs are hypothesized to regulate mitochondrial gene expression and mRNA stability (Ro et al., 2013). While mitosRNAs have been detected in mammalian cells, these have not been characterized in avian species to our knowledge. Therefore, the present study was conducted to identify mitosRNAs in avian breast muscle and to determine if differences in abundance may occur in birds exhibiting large differences in growth and muscle development in a line of birds highly selected for growth and development (PeM broilers) and a heritage breed, BPR chickens.

Materials and methods

Ethics statement

The present study was performed in accordance with the recommended guidelines for the care and use of laboratory animals of the National Institutes of Health. All procedures for animal care were followed by animal use protocol that were reviewed and approved by the University of Arkansas Institutional Animal Care and Use Committee (IACUC Protocol #14012).

Breast muscle tissue and RNA extraction

Pedigree male broilers within a single genetic line and RNA extraction were described previously (Bottje et al., 2002; Kong et al., 2017). Briefly, PeM and BPR chickens (6–8 weeks old, n = 6 per breed) were humanely killed by cervical dislocation. Breast muscle tissue was rapidly excised and immediately flash frozen in liquid nitrogen. Total RNA was isolated from the tissue using TRIzol reagent (Invitrogen Life Technologies, Thermo-Fisher Scientific, Carlsbad, CA) according to the manufacturer's instructions. After initial extraction, the RNA samples were treated with DNase I (New England Biolabs Inc., Ipswich, MA) and extracted a second time with TRIzol reagent. Assessing RNA quality using an Agilent 2200 TapeStation instrument (Santa Clara, CA) revealed that all RNA samples showed >8.0 values of RNA integrity number (RIN). These RNA samples were then used for small RNA sequencing (RNAseq) analysis outlined below.

RNAseq and data analysis

Small RNAseq library preparation for individual samples with barcoding and Illumina sequencing was carried out at the Research Technology Support Facility (Michigan State University, East Lansing, MI). The Illumina HiSeq system at this facility used a 1 × 50 bp single end read technology for sequencing the RNA samples. The RNA sequence FASTQ files that were obtained were mapped to the chicken mitochondrial genome (GenBank Accession: X52392.1) using the NGen program in Lasergene software package (DNAStar, Madison, WI). The total mapped counts were transformed to log2 values based on the number of reads per million (RPM) in order to stabilize the variance. The resulting normalized values were then subjected to further statistical analyses using the JMP Genomics 8 (SAS Institute Inc., Cary, NC). The small RNAs aligned to mitochondrial genome showing >100 average read counts derived by six individual samples, adjusted p-value (false discovery rate; FDR) <0.05 after t-test and log2 fold change >0.5 were considered as differentially expressed (DE) between PeM and BPR. The p-value correction (FDR calculation) was performed by multiple tests of Benjamini and Hochberg method (1995) using JMP Genomics 8.

Nomenclature of identified mitosRNAs

The mitosRNA sequences were identified based on previous reports (Ro et al., 2013; Ma et al., 2016), and their naming was simplified somewhat to contain five parts: [species]-mitosRNA-[strand]-[location]-[variants/fragments]. As an example, gga-mitosRNA-H-2344-1 refers to chicken mitochondrial small RNA derived from the 2344 locus on the heavy strand of the first fragment, whereas gga-mitosRNA-L-16707 refers to chicken mitochondrial small RNA derived from the 16707 locus on the light strand. If multiple fragments determined to have been derived from the same gene locus, the fragments were sequentially numbered.

Quantitative PCR of mitosRNA

Quantitative PCR (qPCR) was conducted following the method reported previously (Ro et al., 2013; Ma et al., 2016) with modifications. Briefly, 60 μg of total RNA obtained from 6 muscle samples from PeM and BPR breast muscle was used for general validation of the RNAseq results and for specific confirmation of the most abundant mitosRNAs for each mitochondrial gene region. Small RNAs were purified from total RNA samples using a miRNA isolation kit (mirVana, Thermo Fisher Scientific, Waltham, MA), polyadenylated with a E. coli Poly-A Polymerase (Thermo Fisher Scientific) and reverse transcribed to cDNA using an adapter primer containing poly T residues at 3′-end (Table 1) and SuperScript III reverse transcriptase (ThermoFisher Scientific) following the manufacturer's instructions. The resulting cDNA samples were diluted (1:10) and a portion (2 μL) was used for qPCR determination (total volume of 25 μL) with an ABI prism 7500HT system (ThermoFisher Scientific) with PowerUp™ SYBR® Green Master Mix (ThermoFisher Scientific). Specific oligonucleotide primers for individual mitosRNA and the universal reverse primer are shown in Table 1. The conditions of real-time qPCR amplification were as follows: 1 cycle at 95°C (10 min), 40 cycles at 95°C (15 s each), followed by 60°C for 1 min. The chicken small 5S rRNA gene was used as the internal control. All qPCR reactions were conducted in duplicate and the values of average cycle threshold (Ct) determined for each sample (n = 6 for each PeM and BPR). The 2−ΔΔCt values were calculated using the average ΔCt (PeM) − average ΔCt (BPR) indicating that the BPR group values were set to 1 and used for relative quantification by linear fold-change. The statistical significance (p < 0.05) between PeM and BPR values were determined by t-test.

Table 1

Primer namePrimer sequence
gga_mitosRNA_H_1781_FACATGTATCCGCCTGAGAACT
gga_mitosRNA_H_2346_6_FTCTAGCCCGACAAACTCGTA
gga_mitosRNA_H_4040_FACATGACCCTGCCCACCCTAA
gga_mitosRNA_H_4957_FATGCCTCTGACATACCAGCAT
gga_mitosRNA_H_8083_FAAGTACTCCAACCCGAATTA
gga_mitosRNA_H_8272_FACCAATTACATAGACCTGTC
gga_mitosRNA_H_8487_FAGATGCCCAAGAAGTTGAAC
gga_mitosRNA_H_9163_FTCACTCTAACAAACAACCCT
gga_mitosRNA_H_10617_FTCGGATTTGAAGCAGCAGCCT
gga_mitosRNA_H_11350_FACCCCATCATTCGCCCTTGT
gga_mitosRNA_H_11664_FTCAACTCCCCTCTTAGTACT
gga_mitosRNA_H_15821_FAACGAACAATAACCTTCCGA
5S_rRNA-F1AAGCCTACAGCACCCGGTAT
RTQ-UNI_RCGAATTCTAGAGCTCGAGGCAGG
RT-Primer with adapterCGAATTCTAGAGCTCGAGGCAGGCGACATGG
CTGGCTAGTTAAGCTTGGTACCGAGCTCG
GATCCACTAGTCCT TTTTTTTTTTTTTTTTTTTTTTTTVN
ND1_FACCCTAGCCATCATCCTGTT
ND1_RTCCTGAGACTAGCTCTGACT
ND2_FAGCATAACCAACGCCTGATC
ND2_RGATGTTAGGAGGAGGAGTGT
COX1_FTCCTTCTCCTACTAGCCTCA
COX1_RAGGAGTAGTAGGATGGCAGT
COX2_FAGGCTTTCAAGACGCCTCAT
COX2_RGTGAGATCAGGTTCGTCGAT
ATP8_FATGCCCCAATTAAACCCTTTCCCA
ATP8_RTTAGGTTCATGGTCAGGTTCA
ATP6_FAATTCTCAAGCCCCTGCCTA
ATP6_RAGGAGGCCTAGGAGGTTAAT
COX3_FTAGTTGACCCAAGCCCATGA
COX3_RGTAGGCCCTTTTGGACAGTT
ND3_FTCTACTAAGCGCTGCACTAA
ND3_RAGGGCGATTTCTAGGTCGAA
ND4L_FTCCCCTACACTTCAGCTTCT
ND4L_RTTCGCATGCTGAGAAGGCTA
ND4_FTCGATCAGCCTACACTGACT
ND4_RTGGGATTAGGGTTGCTTCGA
ND5_FAGCCTCAATGGAGAACAAGACA
ND5_RTGTGGCAAGTAGTGTAAGTGA
CYTB_FTGCCTCATGACCCAAATCCT
CYTB_RAGTGTGAGGAGGAGGATTACT
ND6_FAGACAACCCACGCACAAGCT
ND6_RCTAGGTTTTGTCTTGGTGGT
GAPDH_FACAGCAACCGTGTTGTGGAC
GAPDH_RCAACAAAGGGTCCTGCTTCC

Primers used for qPCR.

The first column indicates primer names and the second column shows sequences for the forward and reverse primers.

Determination of mitochondrial contents using qPCR

Total DNAs of muscle samples were purified using Zymo Tissue and Insect DNA mini Prep kit (Zymo Research, Irvine, CA). Ten nanogram of total DNA were subjected to qPCR reaction with primers for mitochondrial genes that included; (a) NADH-ubiquinone oxidoreductase chain (ND): ND1, ND2, ND3, ND4, ND5, ND6, (b) cytochrome oxidase (COX): COX1, COX2, COX3, (c) cytochrome b (CYTB), and (d) ATPase (ATP) 6, and ATP8, as well as nuclear encoded genes (glyceraldehyde-3-phosphate dehydrogenase; GAPDH) using an ABI prism 7500HT system (ThermoFisher Scientific) with PowerUp™ SYBR® Green Master Mix (ThermoFisher Scientific). The specific oligonucleotide primers are listed in Table 1 and quantifying methods described above.

Results and discussion

Identification of breast muscle MitosRNAs in chicken

Twelve small RNAseq libraries were constructed using RNA isolated from breast muscles obtained from PeM and BPR chickens (n = 6 per group) and sequenced with the 1 × 50 bp single end read method. A total of ~41 million 50 bp sequence reads were produced with an average of ~3.3 million reads per sample (data not shown). After quality assessment, adapter trimming of raw reads, and alignment to chicken mitochondrial genome (GenBank Accession: X52392.1) using the NGen program, ~1.4 million reads with an average of ~120,000 reads per sample (~350 × average coverage) were aligned to the entire sequence of the mitochondrial genome (data not shown). About 97% of reads used in alignment were in the sense orientation of the mitochondrial genome (data not shown). Totally 183,416 unique small RNA sequences were identified as potential chicken mitochondrial genome encoding small RNAs (mitosRNAs) (data not shown). After stringent filtering, 117 mitosRNAs showing >100 raw read counts were abundantly produced from all 37 mitochondrial genes (except D-loop region), all of them were in sense orientation, and the length of mitosRNAs ranged from 22 to 46 nucleotides (Supplementary Table 1).

Only four mitosRNAs (gga_mitosRNA_L_6542, gga_mitosRNA_L_16177, gga_mitosRNA_L_16707, and gga_mitosRNA_L_16745) were derived from the light strand while 113 mitosRNAs were derived from heavy strand mitochondrial DNA (Supplementary Table 1). Similar to findings by Ma et al. (2016), the mitosRNAs were not evenly distributed across the mitochondrial genome (Figure 1). Within the D-loop, rRNA, tRNA, and protein coding regions, the protein- and rRNA coding regions produced much higher number of mitosRNAs (Figure 2), which differs from reports in other species that showed higher mitosRNA frequencies in tRNA coding regions (Lee et al., 2009; Mercer et al., 2011; Kumar et al., 2014; Goodarzi et al., 2015; Hirose et al., 2015; Ma et al., 2016). The unique condition of higher mitosRNAs in mitochondrial protein- and rRNA genes suggests that mitochondrial functions of electron transport chain associated with energy expenditure and OXPHOS in avian muscle may be regulated by unique mechanisms different from other species. Interestingly, our earlier results of global gene expression, in which same RNA samples were used, indicated that muscle growth and fiber composition may be regulated by distinct mitochondrial functions in chicken that might be different from other species (Kong et al., 2017).

Figure 1

Figure 2

Differentially expressed (DE) mitosRNA in PeM compared to BPR

Since our previous global gene expression study (Kong et al., 2017) indicated differential functions of mitochondria in muscles between PeM (rapid growth and large muscle mass) and foundational BPR chickens (slow growth and small muscle mass), differentially expressed (DE) mitosRNA between PeM and BPR were analyzed in the present study. Of the 117 potential mitosRNAs, 44 DE mitosRNAs with adjusted p-values (FDRs) <0.05 were identified. All of the DE mitosRNAs were upregulated in PeM breast, except those that were encoded by the 16S-rRNA gene region (Table 2). Similar to the abundance of total mitosRNAs, only 3 mitosRNAs (gga_mitosRNA_H_4040, gga_mitosRNA_H_8272, and gga_mitosRNA_H_8277) were derived from tRNA coding region while 41 mitosRNAs were derived from mitochondrial protein- and rRNA coding regions (Table 2). The mitosRNAs may be expected to regulate expression and stability of mRNA that is the origin of the designated mitosRNA (Ro et al., 2013; Ma et al., 2016). To validate small RNAseq results, a subset of 12 mitosRNAs that were the most abundant in each gene were subjected to qPCR (Table 3). Results indicated that expression levels for 11 of 12 mitosRNAs were well-matched between small RNAseq and qPCR analyses in terms of the direction and magnitude of change (Table 3). One possible reason why the DE value of a mitosRNA was not matched between RNAseq and qPCR may be related to the use of different normalization approaches between the two methods. The mitosRNAs showing similar DE predictions can be considered the most important of our reported data.

Table 2

NameRNA sequenceLengthRegionA*M**FDR
gga_mitosRNA_H_1426GTAGCCCAAGACGCCTTGCTTAAGCC2612S rRNA10.171.790.0376
gga_mitosRNA_H_1781ATCACACATGTATCCGCCTGAGAACTACGAGCACA3512S rRNA13.331.660.0395
gga_mitosRNA_H_1816AACGCTTAAAACTCTAAGGACTTGGCGGTGCCCCA3512S rRNA9.982.230.0143
gga_mitosRNA_H_2346_1CTTGCCCCCCCTCTAGCCCGACAAACTAGTA3116S rRNA10.72−1.470.0165
gga_mitosRNA_H_2346_2CTTGCCCCCCCTCTAGCCCGACAAACTCATAC3216S rRNA10.11−1.500.0218
gga_mitosRNA_H_2346_3CTTGCCCCCCCTCTAGCCCGACAAACTCCTAC3216S rRNA12.02−1.460.0451
gga_mitosRNA_H_2346_4CTTGCCCCCCCTCTAGCCCGACAAACTCGCACC3316S rRNA11.59−1.440.0218
gga_mitosRNA_H_2346_6CTTGCCCCCCCTCTAGCCCGACAAACTCGTACCCTTAACATAAAAA4616S rRNA17.95−1.190.0387
gga_mitosRNA_H_2346_8CTTGCCCCCCCTCTAGCCCGACAAACTCGTCCC3316S rRNA12.23−1.620.0285
gga_mitosRNA_H_2346_10CTTGCCCCCCCTCTAGCCCGACAAACTGGTA3116S rRNA10.49−1.540.0165
gga_mitosRNA_H_2346_11CTTGCCCCCCCTCTAGCCCGACAAACTTGTA3116S rRNA11.21−1.910.0156
gga_mitosRNA_H_3067TTATTAACAGAACTCAACTTATACCCCCA2916S rRNA11.482.360.0080
gga_mitosRNA_H_3421ATAAGACGAGAAGACCCTGTGGAACTTTAA3016S rRNA11.551.750.0218
gga_mitosRNA_H_4040TACCCCGGACATGACCCTGCCCACCCTAACA31tRNA-Leu13.551.160.0156
gga_mitosRNA_H_4554TTAAGCACCCTGGCCATCACCCAAGAACCC30ND111.533.910.0156
gga_mitosRNA_H_4957TAACCCTAGCCTTATGCCTCTGACATACCAGCATACC37ND111.592.390.0092
gga_mitosRNA_H_5727CTAATCGGAGGCTGAATGGGCCTAAACC28ND211.332.200.0218
gga_mitosRNA_H_5954AATACTAAATGCAACTGTAATACTAACCC29ND211.272.820.0048
gga_mitosRNA_H_6276_1CTACAGAAACTTAGGATTAACTGTCACC28ND213.711.620.0156
gga_mitosRNA_H_6678AACCACAAAGACATTGGCACTCTTTACCT29COX110.372.840.0048
gga_mitosRNA_H_7156TAAAACCCCCCGCACTGTCACAATACC27COX111.301.630.0048
gga_mitosRNA_H_8083GAAAAGTACTCCAACCCGAATTAACT26COX111.861.760.0156
gga_mitosRNA_H_8272AAACCAATTACATAGACCTGTCAAGACTA29tRNA-Asp12.711.660.0349
gga_mitosRNA_H_8277AATTACATAGACCTGTCAAGACTAA25tRNA-Asp11.471.760.0080
gga_mitosRNA_H_8375CATCATAGAAGAGCTCGTTGAATTCCAC28COX210.302.600.0218
gga_mitosRNA_H_8487AACACCGTAGATGCCCAAGAAGTTGAACTAATCTGAACC39COX212.412.250.0221
gga_mitosRNA_H_9163AAACTTCTTTCATTCACTCTAACAAACAACCCTGCA36ATP812.921.650.0080
gga_mitosRNA_H_9235TGACCATGAACCTAAGCTTCTTCGACC27ATP88.312.040.0354
gga_mitosRNA_H_9601AACCCTCCGCCTCCTTAGGACACCTACTCCCTGAAGGCACCCC43ATP612.182.040.0266
gga_mitosRNA_H_9757AACTTATCTCTACAGCCACAATCGCCCTACTACC34ATP613.361.440.0080
gga_mitosRNA_H_9800ATCAATCTCCGCCCTAACGGCACTCA26ATP610.672.890.0080
gga_mitosRNA_H_9871AAGCCTACGTCTTCGTCCTCCTCCTAAGCCTCTACTTACA40ATP612.851.500.0156
gga_mitosRNA_H_10617ACTTCGGATTTGAAGCAGCAGCCTGATACTG31COX310.572.680.0218
gga_mitosRNA_H_11350AAACCCCATCATTCGCCCTTGTACC25ND4_L9.382.030.0218
gga_mitosRNA_H_11633AAAACCCTAACCCTCTGAACAGGCATAGACCA32ND411.422.360.0049
gga_mitosRNA_H_11664AAATCTCAACTCCCCTCTTAGTACTCTCCTGC32ND412.932.050.0049
gga_mitosRNA_H_11665AATCTCCACTCCCCTCTTAGTACTCTCCTGC31ND410.702.650.0005
gga_mitosRNA_H_11905AGAACGACTTAGCGCAGGCATTTACC26ND49.481.550.0218
gga_mitosRNA_H_12707CTCTACATACTACTCTCAACCCAACGAGGCACTC34ND410.091.560.0218
gga_mitosRNA_H_12766CAACTCAAACACTCGAGAACATCTTCTCATAACCC35ND410.112.410.0218
gga_mitosRNA_H_12824CTCATCCTCAAACCAGAACTAATCTCAGGAACCC34ND411.782.580.0080
gga_mitosRNA_H_15821AATCTAAACAACGAACAATAACCTTCCGACC31CYTB12.751.980.0048
gga_mitosRNA_H_15860AAACCCTATTCTGACTTCTAGTAGCCAACC30CYTB11.931.870.0218
gga_mitosRNA_H_15904TGAATCGGAAGCCAACCAGTAGAACACCCC30CYTB9.533.170.0218

The differentially expressed mitosRNAs in PeM (N = 6) compared to BPR (N = 6).

*

A denotes Log2 value of average expressions.

**

M means Log2 fold change of differential expressions.

Table 3

Gene symbolPeM (N = 6) vs. BPR(N = 6)
RNAseqqPCRGene
FC*FDRFC*p-value
gga_mitosRNA_H_17813.170.03955.670.000785012 S rRNA
gga_mitosRNA_H_2346_6−2.280.0387−1.360.047653416 S rRNA
gga_mitosRNA_H_40402.230.01562.020.0007420tRNA-Leu
gga_mitosRNA_H_49575.230.00923.550.0000156ND1
gga_mitosRNA_H_80833.400.01563.240.0014579COX1
gga_mitosRNA_H_82723.170.0349−2.400.0255489tRNA-Asp
gga_mitosRNA_H_84874.770.02214.220.0001418COX2
gga_mitosRNA_H_91633.140.00803.030.0000015ATP8
gga_mitosRNA_H_106176.410.02182.480.0007107COX3
gga_mitosRNA_H_113504.090.02184.390.0000091ND4_L
gga_mitosRNA_H_116644.130.00494.620.0000016ND4
gga_mitosRNA_H_158213.940.00485.810.0000010CYTB

Comparison of fold changes (FC) between RNAseq and qPCR.

*

Values denote linear fold changes.

Comparison between DE MitosRNA and DE mRNA

In order to determine the abundance of total mitosRNAs encoded by each mitochondrial gene, the read counts of all mitosRNAs in a gene region were added to represent a gene and then compared between PeM and BPR. In these results, we noted the same differential pattern as that shown in the DE individual mitosRNA (data not shown). Since DE mitosRNAs were mostly found in protein coding genes (mRNA), the DE mitosRNAs were compared to gene expression derived from the RNAseq analysis conducted on the same muscle samples (Kong et al., 2017). In the RNAseq analysis conducted in the previous study (Kong et al., 2017), differential expression of mitochondrial mRNA was not observed between PeM and BPR (Table 4). However, when data from both DE mitosRNAs and non-DE mitochondrial mRNA were considered, distinct phenotypes (rapid growth, greater muscle mass, fiber composition change) associated with mitochondrial functions shown in modern PeM may be regulated by higher mitosRNA production than the mitochondrial mRNAs (Table 4). Thus, potential regulatory activities of mitosRNAs may influence the stability of the source mRNA and, in turn, ultimately affect the protein abundance and protein functionalities in electron transport chain for OXPHOS and energy expenditure. Future studies on mitosRNAs in chicken muscles will be important to characterize the regulation of differential mitochondrial protein abundance and their functions influenced by mitosRNAs.

Table 4

Gene RegionmitosRNAmRNA
Log2FCp-valueLog2FCp-value
ATP61.330.000130.260.195682
CYTB1.540.00030.310.129304
ND21.510.000320.360.198144
ND41.300.000540.050.753202
ND11.340.00116−0.020.928934
ATP81.090.001340.580.078447
ND50.940.004420.140.306487
COX11.090.004810.090.68997
COX30.830.00697−0.070.800796
ND4L1.090.008890.010.978282
COX21.020.01297nana
ND30.900.053640.170.552786
ND6−0.110.57830.830.003621

Comparison of differential expression of mitosRNA and mRNA in gene regions of mitochondrial genome in PeM breast muscle (N = 6) compared to BPR (N = 6).

Differential abundance of mitosRNAs correlated with mitochondrial contents

While higher abundance of mitosRNAs in PeM could be a function or result from higher mitochondrial content in the muscle tissue, a global transcriptomic study conducted on the same tissue samples reported a significant skew of the mitochondrial proteome in the BPR muscle compared to the PeM which suggests that the mitochondrial content was greater in the BPR tissue than in the PeM tissue (Kong et al., 2017). This finding conflicts with the results of the present study. For this reason, mitochondrial contents were determined by qPCR with muscle DNA and primers for various gene regions (e.g., ND1, ND2, ND3, ATP6, ATP8, COX2, COX3 etc.) of mitochondrial genome, compared with stationary contents of nuclear DNA sequences (GAPDH). When the same amount of DNA was used, the average contents of nuclear genes were very similar among samples in a group and even between groups of PeM and BPR (data not shown), suggesting nuclear DNA copy numbers are the same between samples. In contrast, the average content of mitochondrial genes are slightly higher in PeM compared to BPR (Figure 3). Mitochondrial DNA sequences were ~1.58-fold higher in average throughout the mitochondrial genes in PeM samples compared to BPR muscles, indicating a modest increase of mitochondrial content in modern breeds selected for growth and FE (Figure 3). This result supports the argument that the higher abundance of mitosRNAs in PeM muscles may be derived from the higher mitochondrial contents. Further studies are warranted to fully characterize functions of the lower abundance of mitosRNAs derived from 16S ribosomal RNA in PeM.

Figure 3

According to an earlier transcriptomic investigation with the same tissue samples (Kong et al., 2017), upregulated muscle-related genes suggested a possible higher slow muscle fiber composition in PeM compared to BPR. In addition, an independent report by Mutryn et al. (2015) suggested that wooden breast myopathy (which in commercial broilers showed a dramatic upregulation of troponin I1, myoglobin, myosin binding protein C1, and myozenin 2 genes compared to unaffected muscle) may undergo fast-to-slow fiber type switching. Slow-twitch or oxidative fibers usually possess a higher mitochondrial content, slower contraction rates, increased reliance on OXPHOS, higher fatigue resistance, and higher representation of postural muscles. In contrast, fast-twitch or glycolytic fibers have lower mitochondrial content, rapid contractions, decreased reliance on OXPHOS, low resistance to fatigue, and high representation in muscle groups used for directional movement (Mishra et al., 2015). Chickens have white breast muscle composed primarily of fast twitch glycolytic fibers bearing low mitochondrial content (e.g., Kiessling, 1977). Thus, the fact that higher mitosRNAs with higher mitochondrial contents retained in PeM breast may be indicative evidence for fiber type switching to increase slow-twitch fibers in modern broiler breasts.

In mammals, mitosRNAs are known to be generated by unidentified ribonucleases retained in mitochondria since mammalian mitosRNAs contain 5′ phosphate and 3′ hydroxyl termini that are found in microRNA and small nucleolar RNAs (snoRNAs) (Ro et al., 2013). Depending on the orientation, sense mitosRNAs enhance the expression of host mitochondrial genes by targeting anti-sense transcripts. In contrast, anti-sense mitosRNAs inhibit mitochondrial gene expression by targeting sense transcripts (Ro et al., 2013). Chicken muscle mitosRNAs, which were significantly expressed in this study, are all sense orientations: suggesting that the mitosRNAs may enhance the expression of mitochondrial genes by suppressing antisense RNAs in modern PeM broilers compared to unselected BPR. Although the mRNA expression of mitochondrial genes were not different between PeM and BPR, mRNA stability may be influenced by mitosRNA and, in turn protein production may be affected. Thus, further studies on mitochondrial protein encoding and synthetic machinery may need to be conducted in the future in order to fully understand the role that mitosRNA plays in growth and development of muscle in modern broilers.

Statements

Ethics statement

The present study was performed in accordance with the recommended guidelines for the care and use of laboratory animals of the National Institutes of Health. All procedures for animal care were followed by animal use protocol that were reviewed and approved by the University of Arkansas Institutional Animal Care and Use Committee (IACUC Protocol #14012).

Author contributions

WB and BCK conceived and designed this study. BK, SS, DS, and SK analyzed small RNAseq data and performed qPCR. BM, SO, CO, and NA conduct the animal care, sample collection, and data analyses. JP and JK contribute to data analyses associated to mitochondria. The paper was written through contributions and critical review of the manuscript by all authors (WB, BK, SS, DS, BM, SO, JP, SK, CO, NA, JK, and BCK).

Funding

USDA-NIFA (#2013–01953), Arkansas Biosciences Institute (Little Rock, AR) and the Agricultural Experiment Station (University of Arkansas, Fayetteville).

Acknowledgments

We thank University of Arkansas Research and Sponsored Program and Cell and Molecular Biology Graduate Program for providing site licenses for Lasergene Genomics Suites (DNAStar, Madison, WI) and JMP Genomics 8 (SAS Institute, Cary, NC) for data analyses.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer DC and handling Editor declared their shared affiliation.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2017.00816/full#supplementary-material

Supplementary Table 1

The 117 chicken muscle mitosRNA.

References

  • 1

    BenjaminiY.HochbergY. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Statist. Soc. Ser.57, 289–300. 10.2307/2346101

  • 2

    BottjeW.BrandM. D.Ojano-DirainC.LassiterK.ToyomizuM.WingT. (2009). Mitochondrial proton leak kinetics and relationship with feed efficiency within a single genetic line of male broilers. Poult. Sci.88, 1683–1693. 10.3382/ps.2009-00100

  • 3

    BottjeW. G.CarstensG. E. (2009). Association of mitochondrial function and feed efficiency in poultry and livestock species. J. Anim. Sci.87, E48–E63. 10.2527/jas.2008-1379

  • 4

    BottjeW. G.KongB. W.SongJ. J.LeeJ. Y.HargisB. M.LassiterK.et al. (2012). Gene expression in breast muscle associated with feed efficiency in a single male broiler line using a chicken 44K microarray. II. Differentially expressed focus genes. Poult. Sci.91, 2576–2587. 10.3382/ps.2012-02204

  • 5

    BottjeW.IqbalM.TangZ. X.CawthonD.OkimotoR.WingT.et al. (2002). Association of mitochondrial function with feed efficiency within a single genetic line of male broilers. Poult. Sci.81, 546–555. 10.1093/ps/81.4.546

  • 6

    BottjeW.KongB. W. (2013). Cell Biology Symposium: feed efficiency: mitochondrial function to global gene expression. J. Anim. Sci.91, 1582–1593. 10.2527/jas.2012-5787

  • 7

    BottjeW.KongB. W.ReverterA.WaardenbergA. J.LassiterK.HudsonN. J. (2017). Progesterone signalling in broiler skeletal muscle is associated with divergent feed efficiency. BMC Syst. Biol.11:29. 10.1186/s12918-017-0396-2

  • 8

    BottjeW.PumfordN. R.Ojano-DirainC.IqbalM.LassiterK. (2006). Feed efficiency and mitochondrial function. Poult. Sci.85, 8–14. 10.1093/ps/85.1.8

  • 9

    DorshorstB. J.AshwellC. M. (2009). Genetic mapping of the sex-linked barring gene in the chicken. Poult. Sci.88, 1811–1817. 10.3382/ps.2009-00134

  • 10

    Fernandez-SilvaP.EnriquezJ. A.MontoyaJ. (2003). Replication and transcription of mammalian mitochondrial DNA. Exp. Physiol.88, 41–56. 10.1113/eph8802514

  • 11

    GoodarziH.LiuX.NguyenH. C.ZhangS.FishL.TavazoieS. F. (2015). Endogenous tRNA-derived fragments suppress breast cancer progression via YBX1 displacement. Cell161, 790–802. 10.1016/j.cell.2015.02.053

  • 12

    HiroseY.IkedaK. T.NoroE.HiraokaK.TomitaM.KanaiA. (2015). Precise mapping and dynamics of tRNA-derived fragments (tRFs) in the development of Triops cancriformis (tadpole shrimp). BMC Genet. 16:83. 10.1186/s12863-015-0245-5

  • 13

    KiesslingK. H. (1977). Muscle structure and function in the goose, quail, pheasant, guinea hen, and chicken. Comp. Biochem. Physiol. B.57, 287–292. 10.1016/0305-0491(77)90055-4

  • 14

    KongB. W.HudsonN.SeoD.LeeS.KhatriB.LassiterK.et al. (2017). RNA sequencing for global gene expression associated with muscle growth in a single male modern broiler line compared to a foundational barred plymouth rock chicken line. BMC Genomics18:82. 10.1186/s12864-016-3471-y

  • 15

    KongB. W.LassiterK.Piekarski-WelsherA.DridiS.Reverter-GomezA.HudsonN. J.et al. (2016). Proteomics of breast muscle tissue associated with the phenotypic expression of feed efficiency within a pedigree male broiler line: i. highlight on mitochondria. PLoS ONE11:e0155679. 10.1371/journal.pone.0155679

  • 16

    KongB. W.SongJ. J.LeeJ. Y.HargisB. M.WingT.LassiterK.et al. (2011). Gene expression in breast muscle associated with feed efficiency in a single male broiler line using a chicken 44K oligo microarray. I. Top differentially expressed genes. Poult. Sci.90, 2535–2547. 10.3382/ps.2011-01435

  • 17

    KumarP.AnayaJ.MudunuriS. B.DuttaA. (2014). Meta-analysis of tRNA derived RNA fragments reveals that they are evolutionarily conserved and associate with AGO proteins to recognize specific RNA targets. BMC Biol.12:78. 10.1186/s12915-014-0078-0

  • 18

    LeeY. S.ShibataY.MalhotraA.DuttaA. (2009). A novel class of small RNAs: tRNA-derived RNA fragments (tRFs). Genes Dev.23, 2639–2649. 10.1101/gad.1837609

  • 19

    LopezG.de LangeK.LeesonS. (2007). Partitioning of retained energy in broilers and birds with intermediate growth rate. Poult. Sci.86, 2162–2171. 10.1093/ps/86.10.2162

  • 20

    MaH.WeberG. M.WeiH.YaoJ. (2016). Identification of mitochondrial genome-encoded small RNAs related to egg deterioration caused by postovulatory aging in rainbow trout. Mar. Biotechnol.18, 584–597. 10.1007/s10126-016-9719-3

  • 21

    MercerT. R.NephS.DingerM. E.CrawfordJ.SmithM. A.ShearwoodA. M.et al. (2011). The human mitochondrial transcriptome. Cell146, 645–658. 10.1016/j.cell.2011.06.051

  • 22

    MishraP.VaruzhanyanG.PhamA. H.ChanD. C. (2015). Mitochondrial dynamics is a distinguishing feature of skeletal muscle fiber types and regulates organellar compartmentalization. Cell Metab.22, 1033–1044. 10.1016/j.cmet.2015.09.027

  • 23

    MutrynM. F.BrannickE. M.FuW.LeeW. R.AbashtB. (2015). Characterization of a novel chicken muscle disorder through differential gene expression and pathway analysis using RNA-sequencing. BMC Genomics16:399. 10.1186/s12864-015-1623-0

  • 24

    NiemannH.KuhlaB.FlachowskyG. (2011). Perspectives for feed-efficient animal production. J. Anim. Sci.89, 4344–4363. 10.2527/jas.2011-4235

  • 25

    RoS.MaH. Y.ParkC.OrtogeroN.SongR.HennigG. W.et al. (2013). The mitochondrial genome encodes abundant small noncoding RNAs. Cell Res.23, 759–774. 10.1038/cr.2013.37

  • 26

    TaylorR. W.TurnbullD. M. (2007). Mitochondrial DNA transcription: regulating the power supply. Cell130, 211–213. 10.1016/j.cell.2007.07.002

  • 27

    ZhaoF. P.FanH. Y.LiG. H.ZhangB. K. (2016). Complete mitochondrial genome sequence and gene organization of Chinese indigenous chickens with phylogenetic considerations. Genet. Mol. Res. 15:gmr.15028200 10.4238/gmr.15028200

Summary

Keywords

mitosRNA, mitochondria, chicken breast muscle, small RNA sequencing, mitochondrial contents

Citation

Bottje WG, Khatri B, Shouse SA, Seo D, Mallmann B, Orlowski SK, Pan J, Kong S, Owens CM, Anthony NB, Kim JK and Kong BC (2017) Identification and Differential Abundance of Mitochondrial Genome Encoding Small RNAs (mitosRNA) in Breast Muscles of Modern Broilers and Unselected Chicken Breed. Front. Physiol. 8:816. doi: 10.3389/fphys.2017.00816

Received

15 August 2017

Accepted

04 October 2017

Published

20 October 2017

Volume

8 - 2017

Edited by

Sandra G. Velleman, The Ohio State University Columbus, United States

Reviewed by

Woo Kyun Kim, University of Georgia, United States; Daniel Lee Clark, The Ohio State University Columbus, United States

Updates

Copyright

*Correspondence: Byungwhi C. Kong

This article was submitted to Avian Physiology, a section of the journal Frontiers in Physiology

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics