Abstract Background Co-sequencing of messenger ribonucleic acid (mRNA) and micro ribonucleic acid (miRNA) across a time series (1, 3, 6, 24, and 48 h post injury) was used to identify potential miRNA-gene interactions during pancreatic injury, associate serum and tissue levels of candidate miRNA biomarkers of pancreatic injury, and functionally link these candidate miRNA biomarkers to observed histopathology. RNAs were derived from pancreatic tissues obtained in experiments characterizing the serum levels of candidate miRNA biomarkers in response to acute pancreatic injury in rats. Results No correlation was discovered between tissue and serum levels of the miRNAs. A combination of differential gene expression, novel delayed anti-correlation analysis and experimental database interrogation was used to identify messenger RNAs and miRNAs that experienced significant expression change across the time series, that were negatively correlated, that were complementary in sequence, and that had experimentally supported relationships. This approach yielded a complex signaling network for future investigation and a link for the specific candidate miRNA biomarkers, miR-216a-5p and miR-217-5p, to cellular processes that were in fact the prominent histopathology observations in the same experimental samples. RNA quality bias by treatment was observed in the study samples and a statistical correction was applied. The relevance and impact of that correction on significant results is discussed. Conclusion The described approach allowed extraction of miRNA function from genomic data and defined a mechanistic anchor for these miRNAs as biomarkers. Functional and mechanistic conclusions are supported by histopathology findings. Electronic supplementary material The online version of this article (10.1186/s12864-018-4657-2) contains supplementary material, which is available to authorized users. Keywords: mRNA, miRNA, Co-sequencing, Pancreatic injury, Biomarkers, Delayed anti-correlation Background MicroRNAs (miRNAs) are short ribonucleic acid (RNA) sequences that are hypothesized to primarily function as negative regulators of gene expression [[27]1]. They accomplish this regulation largely through suppressing translation or catalyzing degradation of messenger RNAs (mRNAs) that serve as the template for protein synthesis in the cell [[28]2]. Observed increases in serum miRNAs during multiple disease, dysfunction, and toxicity scenarios have stimulated interest in their use as non-invasive biomarkers of injury. The highly conserved nature of miRNAs across species and, therefore, their large translational potential have further heightened the scrutiny on these small RNAs especially those that are highly enriched in specific tissues [[29]3]. Recent studies have demonstrated the potential of several pancreas-enriched miRNAs as non-invasive biomarkers of acute tissue injury with promise of high sensitivity and specificity [[30]4–[31]9]. Tools are constantly evolving to identify and provide validated interaction information for miRNAs and the mRNAs that they regulate [[32]10]. Within these relationships lay the foundations for the hypothesis and design of future mechanistic investigations. Next Generation Sequencing (NGS) provides an extremely powerful tool for generating comprehensive datasets for discovery and characterization of these interactive relationships [[33]11, [34]12]. As the capability evolves to better analyze and interpret these NGS datasets, so does the ability to more completely capture all of the signaling, process, and regulatory elements in a given biological location at a specific point in time. Faced with these vast and complex NGS datasets, the scientific community, much as it did with microarray gene expression data, must refine quality parameters for, as well as the utility and shortcomings of, these datasets [[35]13–[36]16]. Previously, the caerulein model of acute pancreatic injury was used in rats and the magnitude and temporal responses of miR-216a-5p and miR-217-5p were characterized and their relationship to histopathology defined changes was reported [[37]6]. These two miRNAs were evaluated because they are highly enriched in the pancreas. Therefore, any significant changes in their circulating levels would likely be due to pancreas specific injury and would represent a non-invasive and highly tissue specific marker of cellular injury. The present study was designed to interrogate paired NGS-generated miRNA and mRNA data across the experimental time points in this previous study with three objectives. The first objective was to identify the data defined signaling network, relevant canonical pathways, and potential miRNA-mRNA interactions during early, mild pancreatic injury. The second objective was to characterize the relationship of serum levels of specific candidate miRNA biomarkers of pancreatic injury, miR-216a-5p and miR-217-5p, to their tissue levels during pancreatic injury and recovery. The final objective was to interrogate the refined dataset and determine the miRNA-mRNA relationships of these two candidate miRNA biomarkers and compare the theoretical outcome of regulation of gene expression by these miRNAs to actual histopathology observations in the pancreas. To achieve the present objectives, total RNA was exacted from frozen pancreatic samples generated in that previous study [[38]6]. Extracted RNA was sequenced in separate runs for large and small RNA species to capture mRNA and miRNA results, respectively. This paired NGS data was analyzed for temporal change in expression patterns in both miRNA and mRNA over the first 48 h following caerulein induced pancreatic injury. The subset of miRNA that demonstrated significant expression change following treatment was examined for negative correlations in expression change with their predicted mRNA targets using a novel delayed anti-correlation approach. To satisfy the first study objective, pathway analysis was completed for the gene targets of all negatively correlated miRNA-mRNA pairs to identify potential components of a large miRNA-mRNA regulatory network in early pancreatic injury. The second objective was addressed by correlating quantitative real time polymerase chain reaction (qRT-PCR) measured changes of miR-216a-5p and miR-217-5p in serum with their NGS measured expression change in tissue across the same time points. The final objective was accomplished by demonstrating negative correlation between expression changes in miR-216a and miR-217 and expression changes in their predicted gene (mRNA) targets and then supporting these negative associations through the experimental literature. Subsequent investigation identified the functional pathways and processes associated with these literature-supported genes implicating these miRNAs in the regulation of those pathways and processes. Results Details of all the animal study methodology that yielded the samples from which the present molecular study were derived and the other findings from the animal study including all histopathology findings, characterization of miRNA levels in serum following acute pancreatic injury, and the relationship of the miRNAs to histopathology, have been previously published [[39]6]. Figure [40]1 diagrams the analytical workflow and provides summary data for the present study on the scope of miRNA, mRNA, and potential miRNA-mRNA interactions defined in the data set. NGS identified 126 miRNAs and 12,586 mRNAs that were differentially expressed between treated and control samples in at least one time point (p ≤ 0.05). These data are supplied as additional materials (Additional files [41]1 and [42]2). Time series analysis of these differentially expressed RNAs revealed 4 distinct miRNA profiles (Fig. [43]2) and 6 distinct mRNA profiles (Fig. [44]3) during pancreatic injury and recovery. An integrated analysis combining computational miRNA target prediction, experimental miRNA target database interrogation, and miRNA/mRNA delayed anti-correlation of time profiles identified 619 miRNA-mRNA pairs (Additional file [45]3) that formed a regulatory network of miRNA-mRNA relationships with most miRNAs impacting multiple mRNAs. Pathway analysis suggested this network would affect cell survival/death-related pathways (e.g., Fibroblast Growth Factor and Epidermal Growth Factor mediated pathways) and immune response pathways (e.g., IL-2 and IL-6 signaling pathways). The top pathways associated with the network are listed in Table [46]1. Fig. 1. Fig. 1 [47]Open in a new tab Research Scheme and Results of mRNA and miRNA Co-Sequencing during Pancreatic Injury. Co-sequencing revealed 126 miRNAs and 12,586 mRNA that were differentially regulated at some time in the experimental time course. Differentially expressed miRNAs were then linked to their predicted targets. Those predicted miRNA-mRNA pairs that had an anti-correlation formed a 619-pair signaling network active during pancreatic injury Fig. 2. Fig. 2 [48]Open in a new tab Temporal patterns of differential miRNA expression during pancreatic injury. Differentially expressed miRNA during pancreatic injury sorted into these 4 temporal patterns Fig. 3. Fig. 3 [49]Open in a new tab Temporal patterns of differential mRNA expression during pancreatic injury. Differentially expressed mRNA during pancreatic injury sorted into these 6 temporal patterns Table 1. Top pathways identified in large miRNA-mRNA defined signaling network Pathway Source % of genes False discovery rate ErbB1 signaling NCI 1.0 1.2e-6 Signaling by NGF Reactome 3.8 7.2e-6 PIP3/AKT signaling Reactome 0.9 7.2e-6 Signaling by EGFR KEGG 0.7 6.2e-7 Signaling by SCF-KIT Reactome 2.9 2.7e-5 IL-2 signaling Reactome 2.6 2.7e-5 Signaling by PDGF NCI 0.5 2.7e-5 Signaling by ERBB4 Reactome 3.0 2.7e-5 Signaling by FGFR3 Reactome 2.6 2.7e-5 IL-6 signaling NCI 0.5 2.9e-5 [50]Open in a new tab Serum and tissue levels of two candidate miRNA biomarkers of acute pancreatic injury, miR-216a-5p and miR-217-5p, were compared across the time course (Additional file [51]4). Although serum levels of miR-216a-5p and miR-217-5p increased dramatically, no significant corresponding changes in tissue levels were detected at early time points. At 24 to 48 h post-treatment, as serum levels returned to or neared baseline levels, significant decreases in tissue levels were observed. No significant correlation was demonstrated between an animal’s tissue and serum levels of these miRNAs (Table [52]2). The relationships of these miRNAs to their presumptive mRNA targets were interrogated in this data set using a novel method that takes into consideration the possible delayed target gene expression change in response to its miRNA regulator (see [53]Methods for details). This analysis reveals 9 mRNA negatively correlated with miR-216a-5p and 9 mRNA negatively correlated with miR-217-5p. These negative correlations are temporally depicted for miR-216a-5p and miR-217-5p in Figs. [54]4 and [55]5, respectively. The delayed anti-correlation scores (see [56]Methods) for the individual miRNA-mRNA relationships are presented in Table [57]3. In Fig. [58]2, miR-216a-5p and miR-217-5p would be found in the group depicted in the upper left in which miRNA values initially increase followed by a gradual decline to or below control levels. Their associated mRNAs are immediately depressed relative to controls or decline after the miRNAs increase and then increase again as the miRNAs decrease. Relative to Fig. [59]3, these mRNAs would fall in clusters 1, 3, or 6. A focused pathway analysis of the down regulated genes associated with miR-216a-5p and miR-217-5p revealed that one mRNA (Pten) was a shared target of the two miRNAs and that all but two (Twistnb and Tmem178b) of the identified genes had some previously documented or hypothesized participation in pathways associated with cell autophagy and/or cell death stimulated by cellular stress (Fig. [60]6). Table [61]4 provides summaries of and references