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

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.

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.
Keywords: mitosRNA, mitochondria, chicken breast muscle, small RNA sequencing, mitochondrial contents 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) . 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 welldocumented (Bottje et al., 2002(Bottje et al., , 2006Bottje 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.

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 log 2 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 log 2 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 TM SYBR R 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.

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 proteinand 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 .

Differentially Expressed (DE) mitosRNA in PeM Compared to BPR
Since our previous global gene expression study  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.

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 . In the RNAseq analysis conducted in the previous study , 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 FIGURE 2 | Relative expression of total mitosRNAs in each mitochondrial gene. X axis indicates mitochondrial coding regions, while Y axis represents depth of coverage (raw sequence read counts as results of small RNA sequencing). differential mitochondrial protein abundance and their functions influenced by mitosRNAs.

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 . 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. According to an earlier transcriptomic investigation with the same tissue samples , 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 fastto-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.

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).