Deep Sequencing of Guinea Pig (Cavia Porcellus) Lung Small RNAs Reveals Differential Expression of microRNAs between Non-infected and BCG-Infected Lungs

Background: Pulmonary tuberculosis (TB) caused by Mycobacterium tuberculosis (Mtb) infection remains a major public health burden worldwide. It has been well documented that a group of small noncoding RNAs, microRNAs (miRNAs) are involved in the development and pathogenesis many diseases, including the TB. Guinea pigs are considered as one of the best animal models for biomedical research in TB, the potential roles of miRNAs in the innate immune regulation of guinea pig lung against Mtb infection are not well understood. Methods: In this study, we investigated the differential expression of miRNA proles between the uninfected lungs and Mycobacterium bovis bacillus Calmette-Guérin (BCG)-infected lungs of guinea pigs via deep sequencing and bioinformatics analysis. Results: A total of 2508 miRNAs were identied, among them 1187 were conserved miRNAs and 56 were novel miRNAs in the uninfected lungs, and 1202 were identied as conserved miRNAs and 63 were novel miRNAs in the BCG-infected lungs. Interestingly, comparison analysis further identied 902 co-expressed miRNAs and 585 distinct miRNAs between these two groups. Of the 15 most abundantly conserved miRNAs in guinea pig lungs, which belong to 7 miRNA families, including miR23, miR29, miR145, miR320, miR378, miR451, and miR423. 13 of these 15 most abundant miRNAs were signicantly downregulated and 2 of them were signicantly upregulated in the BCG-infected lungs. Individually, miRNA Let-7f-5p, let-7f, let-7-5p and let-7b-5p were the most abundant in both proles of the non-infected and BCG-infected guinea pig lungs. The predicted target genes of specic miRNAs found in guinea pig lungs were involved in regulation signaling pathways related to immune responses, including Toll-like receptors (TLRs), nuclear factor (NF)-kappa B, Wnt, mitogen-activated protein kinase (MAPK), and transforming growth factor (TGF)-beta signaling, as well as related to autophagy signaling mTOR and apoptotic molecule p53. Conclusions: These data of comprehensive analysis of miRNA transcriptome demonstrated the differential expression proles of miRNAs during M. tuberculosis infection of guinea pig lungs. These results could facilitate the future exploitation of the roles of miRNAs in regulation of immune responses to M. tuberculosis infection using the guinea pig model. guinea pig lungs infected by BCG may have induced the downregulation of miR-145 and inhibited the apoptosis of the host cell. When miR-145-5p was overexpressed in H9c2 cells, the inammatory factor of IL-1beta, TNF-alpha, and IL-6 production triggered by hypoxia were suppressed effectively, and apoptosis was suppressed via the negative regulation of the target gene CD40. This proved that miR-145-5p provides a strategy for the modulation of the inammatory response and apoptosis [46]. Furthermore, These results suggest that miR-944 may be used as a novel biomarker for improving prognosis and as a potential therapeutic target [47] and tumorigenesis [48] by targeting CDK9. These results suggest that miR-145, miR-145-5p, miR-145a-5p, miR-423 and miR-423-5 may play important roles in guinea pig lungs against M. tuberculosis infection.


Introduction
MicroRNA (miRNAs) are a class of short (21 ~ 25 nt), single-stranded, noncoding RNAs that bind the 3' untranslated region of messenger RNA (mRNA) in a sequence-speci c manner, forming an RNA-induced silencing complex (RISC) that either degrades target mRNAs speci cally or suppresses the translation of target mRNAs, and previous studies reveal that miRNAs have a direct inhibitory effect on at least 30% of the protein-encoding genes within the genome [1]. In recent years, studies on how miRNAs function in gene regulation have gained an increased interest in the eld of life science. miRNAs are extensively expressed in various systems, organs, tissues and cells, and play important physiological and pathological roles in cells by regulating signaling transduction, including signaling pathways involved in immune regulation.
Indeed, a compelling body of studies has demonstrated the potential roles of miRNAs in the regulation of immune systems by targeting key signal transduction molecules, since the rst report of miRNA-based regulation of the immune response was published [1]. In this context, miRNAs are able to control the growth, differentiation and function of immune cells, bridge the innate and adaptive immune responses upon an infection of microorganism pathogens. Intriguingly, pathogenic microorganisms in turn can alter the host immune gene signal regulatory networks by directly or indirectly encoding their own miRNAs or changing the target cell miRNAs, which allow them to escape host immunological surveillance [2]. Tuberculosis (TB) is an old disease with new challenges, despite its pathogen Mycobacterium tuberculosis (Mtb) has been discovered for almost one and half century, Particularly, the pathogenic mechanism of immune regulation between host cells and pathogen remains largely unknown [3,4] . Mechanistically, experimental studies have shown that signaling pathways regulating macrophage proliferation or apoptosis may be affected by the mycobacteria in the immunomodulatory process, which can also be regulated by miRNAs. For examples, the mRNA expression of caspase 3 and caspase 7 was downregulated in human macrophages after 48 hours of Mycobacterium avium infection, accompanying with an increased expression of let-7e and miR-29a [5]. Similarly, miR-144 is one of 28 overexpressed miRNAs that is mainly expressed in the T cells of patients with active tuberculosis, and it may play a role in anti-TB immunity by inhibiting the cytokine producing of TNF-α and IFN-γ and proliferation of T cells [6] .
The in vivo study in the transgenic mice expressing a 'sponge' target to compete with endogenous miR-29 exhibited higher serum concentrations of IFN-γ, enhanced T helper type 1 (T(H)1) responses and greater resistance to infection with Mycobacterium bovis bacillus Calmette-Guérin (BCG) and Mycobacterium tuberculosis, as miR-29 was able to inhibit the proliferation of IFN-γ-producing natural killer cells, CD4(+) T cells and CD8(+) T cells [7]. miR-155 and miR-155* are another miRNAs able to be induced by TBspeci c antigens, suggesting that they can be potential diagnostic markers under the challenge of speci c MTB antigens [8].
Animal models are important tools in biomedical researches on pathogenic mechanisms of TB. Among them, mice are overtly the most widely used experimental species for Mtb infection study. Although the reagents used in murine model studies are more diverse because of their relatively resistant properties, the mice are not susceptible to infection with Mtb, and the typical symptoms of TB patients and delayed type hypersensitivity (DTH) are not observed in mouse TB model. In addition, the lungs of mice do not appear to show the typical clinical symptoms of TB, such as coughing, lung cavity alterations, or other pathological changes even in the presence of high numbers of Mtb bacilli. Therefore, mice are not the ideal model for understanding the mechanism of TB pathogenesis [4]. On the other hand, however guinea pigs are very susceptible to Mtb and usually develop lung injury associated with a wide range of chronic diseases seen in human TB lungs, when they are infected with a low-dose aerosol of Mtb. For this reason, guinea pigs are used as the sentinel animals for airborne Mtb contamination [9]. Their study showed that when guinea pig lungs were infected with Mtb, a striking resemblance to the incubation period of TB emerged. After Mtb infection, the guinea pig lungs develop pathological features, such as cheese-like, calci ed granulomas, basic and hemorrhagic lung injury, pulmonary brosis and dissemination. In this context, the guinea pig lung cavity undergoes structural alterations after low-dose and chronic infection with Mtb, and tuberculosis in guinea pigs is a blood-borne disease process. Therefore, the guinea pig is one of the most reliable animal models for the study of tuberculosis immune mechanisms [10].
In the present work, we generated two small RNA libraries by deep sequencing of small RNA transcriptomes of lungs in non-infected and BCG-infected guinea pig. Based on the analysis of small RNA libraries created in this study, we preformed the miRNA expression pro les and differential analysis in conserved and novel miRNAs between these two groups Furthermore, the target genes and signaling pathways of identi ed miRNAs were also predicted by pathway analysis. This work may provide useful information for future study in the innate immunity roles of miRNAs against M. tuberculosis infection in the guinea pig model, particularly it may help us to identify miRNAs candidates involving the innate immunity against M. tuberculosis infection.

Deep Sequencing of Guinea Pig Lung Small RNAs
To identify the miRNAs involved in immune regulation of Mtb infection in guinea pig lungs, small RNA libraries from the non-infected and BCG-infected lungs were sequenced side-by-side using the Illumina HiSeq™ 2000 sequencing system. The statistics of the small RNA sequences from the two libraries were summarized in Table 1 and Table 2. A total of 9,289,133 and 9,600,000 raw reads were obtained from non-infected lungs and BCG-infected lungs, respectively. After removing the low-quality sequences, adapter sequences, and sequences smaller than 18 nt, the remaining clean reads from the two libraries were aligned to the guinea pig genome. A total of 9,237,270 and 9,551,301 clean reads were perfectly matched to the genome ( Table 3). The two libraries generated 302,562 unique small RNA reads, indicating a highly complex small RNA population. For the total small RNAs, 98.43% were found in both libraries. For the unique small RNAs, only 16.09% were shared in both libraries, 42.86% were non-infectedlung-speci c and 41.05% were BCG-infected-lung-speci c, indicating a dynamically similar small RNA population expression in the non-infected and BCG-infected lungs of guinea pigs (Fig. 1). The size of small RNAs was not evenly distributed in both libraries; however, the overall size distribution of the sequenced reads from the two sequencing efforts were very similar. The majority (97%) of small RNAs from both libraries were 20-24 nt, with 22 nt being the most abundant, followed by 23 and 21 nt classes ( Fig. 2). This result was consistent with a recent report in the mouse model [11,12] and with two results in rat [13], human [14] and porcine [12,15] models.

Conserved miRNAs
Conserved miRNA families have been found in many mammalianspecies and have unique expression pro les in the cells of the innate and adaptive immune systems, in which they play pivotal roles in the regulation of both cell development and function [16]. To identify the conserved miRNAs in our dataset, all small RNA sequences were BLASTN searched against the known mature miRNAs and their precursors in the miRNA sequence database, miRBase. There were currently 255 families containing 1187 and 1202 known, unique miRNAs in non-infected and BCG-infected guinea pig lungs, respectively, which match the mouse genome in miRBase 20 (Additional File 1 Table S1). The statistics of conserved miRNA families in non-infected and BCG-infected guinea pig lungs were listed in Additional File 2 Table S2 [17].
In the non-infected and BCG-infected guinea pig lung groups of the miRNA expression pro le, we observed that the majority of miRNAs with abundant counts are represented by a few miRNAs. As shown in Fig. 3, the miRNA expression pro le consisted of unevenly distributed counts of sequences in which the top 13 unique miRNAs with the highest expression levels account for 82.47% and 82.98% (by count) of the total counts of all 13,144,294 and 13,606,027 unique miRNAs of the mappable sequences in noninfected and BCG-infected guinea pig lung libraries, respectively. Interestingly, the histogram showed similar results for both the non-infected and BCG-infected miRNA libraries, with same unique miRNAs making up the top 13 positions. The miR-140 family (miR-140, miR-140-3p) takes the top ranking in both non-infected (257,733) and BCG-infected (277,646) guinea pig lungs, and these account for 1.96% and 2.04% (by count) of mappable sequences, respectively. Furthermore, 7 of top 13 miRNAs are from the let-7 family (Fig. 3). This nding is consistent with previous reports indicating that the let-7 family miRNAs are widely expressed in various cell and tissue types at high expression levels and are involved in the basic process of cell life activities as major regulators of cell proliferation pathways [18,19].
However, the relative abundance of these miRNAs in non-infected and BCG-infected lungs of guinea pig may also suggest their basic cellular roles and that they may be the most important regulatory miRNA groups during the resistance to foreign pathogenic microorganism invasion in lung [20]. In addition, recent evidences also support our hypothesis. Indeed, the anomalous expression of these top 13 ranked miRNAs may be related to the body's immune response [21,22]. These ndings may support the prediction that the 13 miRNAs are related to guinea pig lung physiology, indicating that these may be the most important regulatory miRNA groups and a signature of infectious lung diseases.
It is well known that miRNAs often have implications in signal transduction regulating the binding of the 3' untranslated region (UTR) of its target mRNA, which would be translated into the key signal protein molecules at a later time in immune-related cells. The miRNAs that are differentially expressed between non-infected and BCG-infected guinea pig lungs may represent a major regulationory mechanism for immune responses that could contribute to BCG infection in the lung. To determine the signi cance of differences in known miRNA expression between non-infected and BCG-infected lungs, we showed the expression of miRNA in two samples by plotting the Log2-ratio gure and Scatter Plot ( Fig. 4A and Additional le 3 Table S3). As shown in Fig. 4B, of 1487 unique miRNAs, 902 (60.66%) unique miRNAs were co-expressed in the non-infected and BCG-infected guinea lung libraries, while 285 (19.17%) and 300 (20.30%) miRNAs appear to be preferentially expressed in the non-infected and BCG-infected libraries, respectively. The analysis of library sequencing data resulted in the identi cation of 381 unique miRNAs (out of 1187; 32.10%) with statistically signi cant differential expression between the non-infected and BCG-infected guinea pig lung libraries. Out these 381 differential expression unique miRNAs, 248 (Noninfected-speci c: 42, co-expressed: 206) and 133 (BCG-infected-speci c: 55, co-expressed: 78) unique miRNAs are down-and upregulated in non-infected versus BCG-infected lungs, respectively (Fig. 4B).
In the miRNA differential expression libraries of the non-infected and BCG-infected guinea pig lung, the 1487 unique miRNAs belong to 55 of the 255 identi ed conserved miRNA families. In addition, of the 15 most abundant conserved miRNAs in miRNA differential expression libraries of non-infected guinea pig lungs, 13 were signi cantly downregulated and 2 were signi cantly upregulated, compared to those of the BCG-infected lungs (Fig. 5). The 13 downregulated miRNAs belong to 6 miRNA families: miR23, miR29, miR145, miR320, miR378, and miR451. The 2 upregulated miRNAs belonged to the miR423 family (Fig. 6). These data was in line with previous ndings that miR-29 was expressed during active pulmonary tuberculosis [23], while miR-23, miR-29, miR-145, and miR -452 were either up-or downregulated in chronic obstructive pulmonary disease (COPD), non-small cell lung cancer (NSCLC) and in ammatory responses [24][25][26][27].

Novel miRNAs
To identify more potential miRNAs in the BCG-infected guinea pig lungs, the unclassi ed reads were further processed by Mireap software (http://sourceforge.net/projects/mireap). Only those tags meeting the default parameters and with read counts greater than 56 and 63 in non-infected and BCG-infected guinea pig lungs, respectively, were de ned as candidate novel mature miRNAs. The target genes (433 of non-infected and 434 of BCG-infected guinea pig lung) of novel miRNAs were predicted from Mireap, and the number of miRNAs (only miRNAs with predicted target genes will be included here) and the number of corresponding targets is shown in Table 4 (for detailed information on the miRNA sequence and target gene sequence, see Additional File 4 Table S4 for novel miRNAs and Additional File 5 for target prediction of novel miRNA). The rules used for target prediction are based on those suggested by Allen [28] and Schwab et al. [29].
On the basis of this analysis, 79 novel miRNAs were identi ed in the guinea pig lung, and these showed different expression levels between the non-infected and BCG-infected lungs. The expression of miRNAs in these two groups was presented by plotting the Log2-ratio gure and a scatter plot (Fig. 7A). Of the 79 novel unique miRNAs, 21 (26.58%) novel miRNAs were coexpressed in non-infected and BCG-infected guinea lung libraries, while 28 (35.44%) and 30 (37.97%) novel miRNAs appear to be preferentially expressed in non-infected and BCG-infected libraries, respectively. The analysis of library sequencing data resulted in the identi cation of 29 novel miRNAs (out of 79, 36.71%) with statistically signi cant differential expression between the non-infected and infection guinea pig lung libraries. Out of these 29 differentially expressed novel miRNAs, 14 (non-infected lung-speci c) and 15 (BCG-infected-speci c: 13, coexpressed: 2) novel miRNAs were found to be down-and upregulated in non-infected versus BCG-infected lung samples, respectively ( Fig. 7B and Additional File 6 Table S5 show the novel miRNA differential expression).
The Venn diagram displays the distribution of 49 novel miRNAs between non-infected (left, blue circle) and BCG-infected guinea pig lung (right, red circle) libraries. The dashed circles indicate the miRNAs that were signi cantly differentially expressed (p < 0.0001) in the two lung samples.

Nucleotide bias of miRNAs in guinea pig lungs
We analyzed the nucleotide bias of the conserved and novel miRNAs in non-infected and BCG-infected guinea pig lungs. U was the most common rst nucleotide at the 5' end of the conserved and novel miRNAs. U was the dominant rst nucleotide of the conserved 20-24 nt miRNAs in accounting for 97% of 49977 in non-infected guinea pig lungs and for 81% of 43124 in BCG-infected lungs. Like the conserved miRNA, the novel 21-24 nt miRNA accounted for 59% of 1871 in non-infected group and 61% of 1697 in the BCG-infected group (Fig. 8). Some novel miRNAs predicted in our miRNA expression pro le may be just small RNAs rather than miRNAs, which could be why the percentage of novel miRNAs was lower than the percentage of conserved miRNAs. Like the other studies, these results also con rmed U to be the most common nucleotide at the 5' end of miRNAs [30].

Prediction of target genes and signaling pathways
In organisms, genes usually interact with each other to play different roles in certain biological functions. An analysis based on pathways could allow us to better understand the biological functions of genes. KEGG is the major public pathway-related database.
The target genes were predicted on the guinea pig (Cavia porcellus) genome sequence (http://asia.ensembl.org/Cavia_porcellus/Info/Index?db=core) to seek out miRNAs involved in the lungs of guinea pigs as a model against the M. tuberculosis infection. A total of 19606 potential target genes were predicted (Additional le 7 Table S6). We identi ed 308 signaling pathways on the Kyoto Encyclopedia of Genes and Genomes (KEGG), including at least 12 pathways (Table 6) involved in the innate immune response of guinea pigs against M. tuberculosis infection, such as BCG, Toll-like receptor (TLR) and NF-kappa B signaling pathways, which could regulate host immune cell recognition and phagocytosis to pathogenic bacteria [31]; NF-kappa B [32], Wnt, MAPK, and TGF-beta signaling pathways are involved in the host in ammatory response and innate immune response [33]. The mTOR and p53 signaling pathways regulate the autophagy and apoptosis processes in immune cells [34]. A number of targets genes were detected. 223 belonged to signaling pathways involving tuberculosis, 119 belonged to the TLR signaling pathway, 158 belonged to the NF-kappa B signaling pathway, and 114 belonged to the TGF-beta signaling pathway (Table 5).
In the innate immunity of hosts infected with M. tuberculosis, TLRs are important pattern-recognition receptors for recognizing microbial pathogen-associated molecule patterns (PAMPs) and ghting aM. tuberculosis attack by initiating an immune response via tactivating TLR, NF-kappa B and many other signaling pathways that regulate the immune and in ammatory gene expression [35]. To better understand the biological functions of the conserved miRNAs in differential expression pro les that target to TLRs signaling pathways, putative target genes of miRNAs were predicted using TargetScan (http://www.targetscan.org/vert_72/) (Fig. 9). The results showed that 52 miRNAs were downregulated and 6 were upregulated on activated TLR and NF-kappa B signaling pathways in BCG-infected guinea pig lung models. This reminds us that these miRNAs may be involved in the host innate immunity response by regulating target genes in many signaling pathways.

Validation of miRNAs by quantitative polymerase chain reaction
To con rm the expression of the conserved miRNAs, quantitative RT-PCR (qRT-PCR) was performed on 6 randomly selected miRNAs, namely, miR-1421 g-5p, miR-3689a-3p, miR-2227, miR-4958-5p, miR-29d-3p and miR-1591-3p in non-infected and BCG-infected guinea pig lungs. Aliquots of the RNA samples used for sequencing were subjected to the qRT-PCR assay, with similar expression levels detected for all the 6 conserved miRNAs candidates by deep sequencing (Fig. 10).

Discussion
Tuberculosis (TB) is caused by Mycobacterium tuberculosis (Mtb) infection and is one of the most fatal infectious diseases. It is a main leading cause of mortality, especially in co-infection with HIV, with 95% of these deaths occurring in developing countries [36]. The causative pathogen of M. tuberculosis has the capability of avoiding the host's immune system for its intracellular survival. Increasing studies have suggested that miRNAs exert its function at the posttranscriptional level and are able to shape immunity via regulating the gene pool expressed in immune-related cells, especially by regulating the host innate immune response against Mtb infection. Guinea pigs are commonly used as tuberculosis vaccine rh, and because the pathogenesis of guinea pigs infected with Mycobacterium tuberculosis shares many similarities with the clinical symptoms of human tuberculosis, they have become an increasingly useful model for tuberculosis research [37,38]. We established a guinea pig model infected with BCG and performed an analysis of miRNA expression pro les in guinea pig lungs by deep small RNA sequencing. The miRNA libraries were constructed from non-infected and BCG-infected guinea pig lungs, which displayed a similar read length distribution. More than 60% of miRNA were 22 nt, followed by 23 and 21 nt, representing the typical dicer-derived products [39].
The most abundant miRNAs in non-infected and BCG-infected guinea pig lungs were let-7f, let-7, let-7b, let-7e, let-7g, miR-140-3p and miR-107, with more than 10000 reads each. Many of them were present in the host immunity regulation against M. tuberculosis infection. M. tuberculosis macrophage infection causes ESAT-6-dependent miRNA let-7f downregulation, and the overexpression of Let-7f in macrophages could reduce M. tuberculosis survival via targets A20 (TNFAIP3) to modulate NF-κB activity during M. tuberculosis infection [40]. The levels of hsa-let-7b-5p were con rmed to be signi cantly upregulated in pulmonary TB patients, compared to both control groups, which suggests that circulating hsa-let-7b might be associated with TB development by regulating the target genes involved in the TLR-NF-kB mediated signal pathway [41]. The miRNA let-7e was downregulated in the blood serum of patients with active pulmonary tuberculosis [42]. The expression of miR-107 target cyclin-dependent kinase 6 (CDK6) is increased by TLR4 as a result of the decrease in miR-107 when macrophages respond to lipopolysaccharide (LPS) [43]. However, miR-140-3p was differentially expressed in macrophages infected in vitro in active TB patients. These miRNAs may be related to the secondary signal transduction pathway, which is likely involved in MTB infection. Four signaling pathways, similar to those predicted by KEGG, including the Wnt signaling pathway, insulin signaling pathway, TGF-β signaling pathway and glycosaminoglycan biosynthesis, are involved in active MTB & LTBI studies [44].
The most differential of the 6 miRNAs in non-infected and BCG-infected guinea pig lungs were miR-145, miR-145-5p, miR-145a-5p and miR-320c, with more than 5000 each. They showed downregulated expression after BCG infection. On the other hand, the expressions of miR-423 and miR-423-5p increased signi cantly. The role of miR-145 in the posttranscriptional regulation of c-Myc by p53 was also assessed. p53 transcriptionally induces the expression of miR-145 by interacting with a potential p53 response element (p53RE) in the miR-145 promoter. In this context, miR-145 provides a direct link between p53 and c-Myc in this gene regulatory network [44]. Additionally, previous results indicated that miR145 over-expression promotes the cleavage of caspase3 and PARP, suppresses proliferative ability and induces the apoptosis of A549 cells, possibly by decreasing the MMP2 and MMP9 expression, the Bax/Bcl2 ratio and the activity of the caspase-3 cascade [45]. These studies were inconsistent with our results in guinea pig lung. One possible reason may be that the guinea pig lungs infected by BCG may have induced the downregulation of miR-145 and inhibited the apoptosis of the host cell. When miR-145-5p was overexpressed in H9c2 cells, the in ammatory factor of IL-1beta, TNF-alpha, and IL-6 production triggered by hypoxia were suppressed effectively, and apoptosis was suppressed via the negative regulation of the target gene CD40. This proved that miR-145-5p provides a strategy for the modulation of the in ammatory response and apoptosis [46]. Furthermore, These results suggest that miR-944 may be used as a novel biomarker for improving prognosis and as a potential therapeutic target [47] and tumorigenesis [48] by targeting CDK9. These results suggest that miR-145, miR-145-5p, miR-145a-5p, miR-423 and miR-423-5 may play important roles in guinea pig lungs against M. tuberculosis infection.
It has been established that miRNAs play a critical role in the host innate immunity responses. An increased number fo evidence has revealed the roles of miRNAs in regulating the host in ammatory responses in reaction to innate immunity signal activation via the TLR signaling pathways. Among them, miR-146a, predicted in our analysis in Fig. 9, was one of the most important miRNA to negatively regulate the host in ammatory signals triggered by the TLR pathway in human hepatic stellate cells. In addition, miR-146a is another best characterized miRNA and is shown to negatively regulate host in ammatory signaling through its target genes, IRAK1 and TRAF6, in the TLR-4 signaling pathways [49], which was also upregulated in alveolar macrophages [50] and some other nonimmune cells by the activation of the TLR pathway, and regulation of in ammatory responses through mechanisms by inhibiting target gene expression [51]. Overall, accumulating miRNAs were shown to be involved in the innate immunity response to regulate in ammatory cytokines and aid in host defenses against the infection of various pathogens, particularly M. tuberculosis [52].

Conclusions
This study provides the rst report of the differential expression pro les of miRNAs in a guinea pig model infected with BCG. It provides informative data regarding the potential involvement of miRNAs in immune regulation of M. tuberculosis in guinea pig lungs. miRNAs regulate many aspects of the immune response, especially in innate immunity and autophagy during mycobacterial infection, which targets genes in the TLR, Wnt, and mTOR signaling pathways. These speci cally expressed miRNAs provide attractive candidates for further studies to reveal the role of miRNAs in M. tuberculosis infection using guinea pig models.

Establishment of a guinea pig model infected with BCG
The SPF Hartley Guinea Pig (Charles River Laboratories, China) were housed in the Experimental Animal Center of Ningxia Medical University (Yinchuan, Ningxia, China) in an standard independent ventilation cage (IVC). The oor space of IVC was 45cm×29cm, and height of 25cm. One guinea pig was housed in each IVC. The animals were bred in the SPF unit with temperature of 21 ± 2°C, relative humidity of 30-70% and a light/dark cycle of 12/12 hour. The animals were fed with ad libitum Guinea pig Breeding & Maintenance Diets (Irradiated) and water animals were fed with ad libitum Guinea pig Breeding & Maintenance Diets (Irradiated) and water. For infection with BCG, the animals were anesthetized with Ketamine (40 mg/kg) + Dexmedetomidine (0.15 mg/kg) by intramuscular injection (i.m). The reversal agent of anesthesia was atipamezole (1 mg/kg) by IM [53].
For experimental designing, six guinea pigs (male: female = 1:1, with weights of 300-400 g and negative results for 10 IU PPD tests) were randomly divided into two groups: 3 were in the control group and 3 were in the BCG-infected group. Infection method: BCG (Bacillus Calmette Guerin, 60 mg/ 6.0×10 7 CFU, S20123007, CDIBP, China) was diluted to 10 6 CFU/mL with 5 mL saline and added to a medical ultrasonic nebulizer (WH-2000 ultrasonic atomizer, Guangdong Yuehua, Guangzhou, China), which was connected to a spray chamber with aerosol inhalation for 90 minutes [54]. After 24 hours, the guinea pigs infected with BCG, and animals were euthanized with overdose pentobarbital (150 mg/kg) by intraperitoneal injection (i. p.) at 48 hours post infection, the lungs were then harvested and immediately frozen in liquid nitrogen and stored at − 80°C for further use [55,56].

RNA isolation
The total RNA was extracted from each guinea pig lung sample separately using TRIzol reagent (Thermo Fisher Scienti c, USA), following the manufacturer's protocol. The purity and quantity of RNA were assessed using a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scienti c) and denaturing gel electrophoresis, and the samples were then stored at − 80°C.

Small-RNA library construction and sequencing
We generated small-RNA libraries of two groups of samples from six lungs of guinea pigs using the mirVana™ miRNA Isolation Kit (Thermo Fisher, USA), according to manufacturer's protocol. Three biological replicates of each group were prepared for each sample, with each biological replicate comprising RNA extracted from 3 guinea pigs. The total RNA was ligated with 3' and 5' RNA adaptors, and fragments with adaptors on both ends were ampli ed by PCR after reverse transcription. The subsequent cDNAs were puri ed and enriched by 6% denaturing polyacrylamide gel electrophoresis to isolate the expected size fractions and eliminate unincorporated primers, primer dimer products, and dimerized adaptors. Finally, the three RNA libraries of each groups were mixed as one sample and were sequenced using an Illumina/Solexa Genome Analyzer at BGI (Huada Genomics Institute Co. Ltd., Shenzhen, China).

Analysis of the sequencing data
The 50 nt sequence reads from deep sequencing were purged of low-quality reads and adaptor sequences to obtain the nal clean reads. Read lengths between 18 nt and 30 nt were used for further bioinformatic analysis. The reads mapping to the guinea pig genome were prepared for the next steps of analysis. The reads mapped to the guinea pig genome were subsequently analyzed to annotate the small RNA reads with rRNA, scRNA, snoRNA, snRNA and tRNA, and degraded mRNA fragments, miRNAs, repeatassociated sRNA, and noncoding RNA sequences were removed by blasting against the Rfam database (http://useast.ensembl.org/Cavia_porcellus/Info/Index) and the GenBank database (http://www.blast.nvbi.nlm.nih.gov/). The remaining sequences were identi ed as the conserved miRNAs in guinea pigs by blasting against the miRBase 21.0, allowing no more than two mismatches. miRNA information of guinea pigs in miRbase 21 was obtained when sequences were aligned to the miRNA precursor/mature miRNA of corresponding species in miRbase. This produced detailed information on the alignment. If there was no miRNA information on guinea pigs in miRBase21, then they were aligned to the miRNA precursor/mature miRNA of all other animals in miRBase21 and shown alongside the sequence information of the miRNA families (no speci c species), which could be found in the samples. Existing guinea pig miRNA existed in the miRBase with no base mismatches. The sequences that did not match existing or conserved miRNAs were used to identify potentially novel miRNA candidates, which were characteristic of hairpin structures of the miRNA precursors [34]. These candidates were identi ed using Mireap (https://sourceforge.net/projects/mireap/) to predict novel miRNAs by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the unannotated small RNA tags, which could be mapped to the genome. The procedures to compare the known miRNA expression between two libraries to identify differentially expressed miRNA are shown below. To non-infectedize the expression of miRNA in the two samples (the control and treatment) to get the expression of transcript per million (TPM), the non-infectedization formula used was as follows: Non-infectedized expression = Actual miRNA count/Total count of clean reads× 1,000,000; Calculate fold-change and P-value from the non-infectedized expression. Then, the Log2-ratio plot and scatter plot were generated. The fold-change formula was as follows: Fold-change = log2(treatment/control) [57]. A positive value indicated miRNA upregulation while a negative value indicated downregulation.
Target mRNAs of the miRNAs were identi ed by using Targetscan software (http://www.targetscan.org/vert_72/). KEGG pathway analysis was used for the target gene candidates [58] , and pathway analysis of the predicted target genes was analyzed using the KEGG pathway database (http://www.genome.jp/kegg/pathway.html).

Quantitative reverse transcriptase PCR
The expression pro les of 6 randomly selected conserved miRNAs were investigated by quantitative PCR to validate their expression changes. The total RNA (800 ng) was converted to cDNA using the Reverse Transcription System (Promega, USA) following the instructions of the manufacturer. Real-time PCR was performed using the Bio-Rad IQ5 instrument according to the standard protocol. The 20-µL PCR reaction mixture included 10 µL SYBR Green Supermix (2×), 0.4 µL miRNA-speci c forward primers (10 µM), 0.4 µL universal primer (10 µM), and 1 µL PCR template (cDNA). The cycle parameters were as follows: 2 minutes at 95℃ 1C, followed by 40 cycles of 15 seconds at 95℃ 1C and 30 seconds at 60℃ 1C [59].
Melting curve analyses were performed to con rm the Speci city of the PCR products. All primers used in the real-time PCR experiments are shown in Additional le 8 Table S7. The relative miRNA expression levels were calculated using the 2-ΔΔCt method after the threshold cycle (Ct). Each sample was run in triplicate. snRNA U6 was used as an endogenous control for the real-time PCR of miRNAs. Declarations The Author con rms : that the work described has not been published before (except in the form of an abstract or as part of a published lecture,review,or thesis); that it is not under consideration for publication elsewhere; that its publication has been approved by all co-authors, if any; that its publication has been approved (tacitly or explicitly) by the responsible authorities at the institution where the work is carried out. Figure 1 Comparison of the common and speci c sequences between non-infected and BCG-infected Guinea pig lungs.

Figure 2
Size distribution of mappable reads related to miRNAs in the non-infected and BCG-infected guinea pig lungs of miRNA libraries Figure 3 Size distribution of mappable reads related to miRNAs in the non-infected and BCG-infected guinea pig lungs of miRNA libraries Comparison of differentially expressed miRNAs between the non-infected and BCG-infected guinea pig lungs.

Figure 5
The 15 miRNAs with the highest (-std) value in the non-infected lung differential expression pro les matched to BCG-infected guinea pig lungs. The 15 miRNAs count with the highest (-std) value in the non-infected lung differential expression pro le belong to 7 miRNA families. Comparison of differentially expressed novel miRNAs between the non-infected and BCG-infected guinea pig lung samples.

Figure 8
First nucleotide bias in different length tags in the conserved and novel miRNAs of non-infected and BCGinfected guinea pig lungs.