Investigating circulating miRNA in transition dairy cows: What miRNAomics tells about metabolic adaptation

Abstract

In the current study, we investigated dairy cows' circulating microRNA (miRNA) expression signature during several key time points around calving, to get insights into different aspects of metabolic adaptation. In a trial with 32 dairy cows, plasma samples were collected on days -21, 1, 28, and 63 relative to calving. Individually extracted total RNA was subjected to RNA sequencing using NovaSeq 6,000 (Illumina, CA) on the respective platform of IGA Technology Services, Udine, Italy. MiRDeep2 was used to identify known and novel miRNA according to the miRbase collection. Differentially expressed miRNA (DEM) were assessed at a threshold of fold-change > 1.5 and false discovery rate < 0.05 using the edgeR package. The MiRWalk database was used to predict DEM targets and their associated KEGG pathways. Among a total of 1,692 identified miRNA, 445 known miRNA were included for statistical analysis, of which 84, 59, and 61 DEM were found between days -21 to 1, 1 to 28, and 28 to 63, respectively. These miRNA were annotated to KEGG pathways targeting the insulin, MAPK, Ras, Wnt, Hippo, sphingolipid, T cell receptor, and mTOR signaling pathways. MiRNA-mRNA network analysis identified miRNA as master regulators of the biological process including miR-138, miR-149-5p, miR-2466-3p, miR-214, miR-504, and miR-6523a. This study provided new insights into the miRNA signatures of transition to the lactation period. Calving emerged as a critical time point when miRNA were most affected, while the following period appeared to be recovering from massive parturition changes. The primarily affected pathways were key signaling pathways related to establishing metabolic and immune adaptations.

Keywords: immune-related pathways; miRNA; negative energy balance; post-calving; posttranscriptional regulation; systemic inflammation.

FIGURE 1

Schematic representation of (A) the animal experimental design and sampling time points, consisting of four groups CTRL: control (n = 8), EFA: essential fatty acids (n = 8), CLA: conjugated linoleic acid (n = 8), and EFA + CLA (n = 8), and four plasma sampling (days −21, +1, +28, and +63 relative to parturition). The focus of this study is primarily on time comparisons since the effect of treatments was negligible (B) Sequencing analysis and bioinformatics. CTRL: control, EFA: essential fatty acids, CLA: conjugated linoleic acid.

FIGURE 2

The top 25 miRNAs (by mean counts (million) per time point) in circulation during days −21, 1, 28, and 63 relative to parturition. Different colors correspond to different time points.

FIGURE 3

Multidimensional scaling (MDS) plot indicates samples (individual sequencing) separation during days -21, 1, 28, and 63 relative to parturition. Different colors correspond to different time points.

FIGURE 4

(A). Volcano plot representing differentially expressed miRNA (DEM) between day -21 and +1 relative to parturition (day +1/-21); increased (red dots) and decreased (green dots) miRNA in d+1/d-21 are highlighted (p < 0.05 and log2 fold change (FC) > 1.3). (B). KEGG pathways annotated to DEM; Bars indicate proportional to the false discovery rate (FDR) adjusted p-value, and the box on each bar represents the gene counts (GC).

FIGURE 5

(A). Volcano plot representing differentially expressed miRNA (DEM) between day +1 and +28 relative to parturition; increased (red dots in top right) and decreased (green dots in top left) miRNA in d+28/d+1 are highlighted (p < 0.05 and log2 fold change (FC) > 1.3). (B). KEGG pathways annotated to DEM; Bars indicate proportional to the false discovery rate (FDR) adjusted p-value, and the box on each bar represents the genes count (GC).

FIGURE 6

(A). Volcano plot representing differentially expressed miRNA (DEM) between day +28 and +63 relative to parturition; increased (red dots in top right) and decreased (green dots in top left) miRNA in d+63/d+28 are highlighted (p < 0.05 and log2 fold change (FC) > 1.3). (B). KEGG pathways annotated to DEM; Bars indicate proportional to the false discovery rate (FDR) adjusted p-value, and the box on each bar represents the genes count (GC).

FIGURE 7

Venn diagram representing the overlap between differentially expressed miRNA during time series comparisons, including Day +1/-21 (blue), Day +28/+1 (red), and Day +63/+28 (green).

FIGURE 8

(A) A bar chart shows differentially expressed miRNAs whose up- or downregulation occurred on day 1 postpartum and returned antepartum levels on day 28 postpartum (B) A bar chart shows differentially expressed miRNAs whose up- or downregulation occurred on day 28 postpartum and returned to day 1 levels on day 63 postpartum.

FIGURE 9

miRNA/mRNA network analysis. The regulation network of downregulated (A) and upregulated (B) miRNAs and predicted target mRNAs in d+1/-21 is illustrated by Cytoscape software. The nodes represent the mRNA (green), and the boxes represent the miRNA (blue). The box size corresponds to the centrality and betweenness of miRNA in constructing the network.

FIGURE 10

miRNA/mRNA network analysis. The regulation network of downregulated (A) and upregulated (B) miRNAs and predicted target mRNAs in d+28/+1 is illustrated by Cytoscape software. The nodes represent the mRNA (green), and the boxes represent the miRNA (blue). The box size corresponds to the centrality and betweenness of miRNA in constructing the network.

FIGURE 11

miRNA/mRNA network analysis. The regulation network of downregulated (A) and upregulated (B) miRNAs and predicted target mRNAs in d+63/+28 is illustrated by Cytoscape software. The nodes represent the mRNA (green), and the boxes represent the miRNA (blue). The box size corresponds to the centrality and betweenness of miRNA in constructing the network.