All animal studies were conducted according to the guidelines of the Institutional Animal Care and Use Committee of Shandong Agricultural University (Approval Number: #SDAUA-2018-028) and the “Guidelines for Experimental Animals” of the Ministry of Science and Technology (Beijing, China). New Zealand rabbits, 38-days-old, were obtained from Qingdao Kangda Rabbit Industry Development Co., Ltd. All nasopharyngeal cultures of these animals tested negative for P. multocida. Moreover, rabbit sera were also negative for P. multocida based on ELISA using a specific antibody against the bacterium. Rabbits were adapted to their environment for one week before use in experiments.

Twenty-five rabbits were randomized into test (n = 20) and control (n = 5) groups. For the challenge assay, rabbits were injected subcutaneously with 10⁷ CFU of strain type A3 (LD 50 = 2.3 × 10^7.8 CFU/mL) in 1mL PBS (24). One day post-infection (P1, dpi) and 3 dpi (P3), five rabbits with symptoms of depression, anorexia, snuffles, serous nasal exudate, and dyspnea were euthanized by pentobarbital overdose (100 mg/kg, intravenous). Unchallenged rabbits (P0) were used as a negative control. All the lung samples were collected under aseptic conditions and immediately stored at −80°C. The remining alive rabbits, post-experiment, were euthanized by the same procedure.

Every three samples we used for the RNA-seq were taken at P0, P1, and P3 to examine the lung tissues at different stages of infection. Total RNA was extracted from lung tissues using a TRIzol Kit (Invitrogen, CA, USA), as recommended by standard procedure. RNA purity and integrity were checked through NanoDrop 2000 Spectrophotometer and Agilent 2100 Bioanalyzer. MiRNA enrichment was done with the PureLink miRNA isolation Kit (Invitrogen, USA) for library construction. Briefly, 15% urea denaturing polyacrylamide gel electrophoresis (PAGE) was employed to purify the miRNAs molecules (<30 nucleotides). Next, miRNAs were ligated with the specified adaptors at the 3′ and 5′ ends. Subsequently, ligated miRNAs were converted into cDNA, and products were successively amplified with PCR. A second process of size selection was performed, and target products were excised from gels. Finally, products were purified and submitted for miRNA transcriptome sequencing (Huada Gene, Shenzhen, Guangdong, China).

Adaptor sequences, redundant contaminants, low-quality reads were eliminated to ensure reliable clean data. Clean reads after data filtering were mapped with known small-RNA data resources—miRbase, snoRNA, Rfam, piRNA, and siRNA—to remove additional non-coding RNAs (rRNA, tRNA, snoRNA, and snRNA) and repeat sequences (25, 26). miRDeep2 (27) was used to predict novel miRNA by exploring secondary structure, and Piano (28) was used to predict piRNAs. miRNA expression level was calculated by Fragments per Kilobase Million (FPKM) method (29). Differential expression (DE) analysis used DEGseq (30) with a threshold of Q value ≤ 0.001 and |Log2 (fold changes)| ≥ 1.

To predict genes targeted by DE-miRNAs, putative 3′-UTR of rabbit mRNA were used with two computational algorithms RNAhybrid (31) and miRanda (32). Moreover, explicit roles of miRNA-targeted genes in distinct biological processes were explored with the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Notable terms and pathways were adjusted using a Q value with a rigorous Bonferroni cutoff (33).

DE-miRNAs identified by sequencing were selected for the validation of expression results. mRNA abundance for target genes was also examined. An miRcute miRNA isolation kit (Tiangen, Beijing, China) was employed to isolate miRNAs, and Trizol reagent (Invitrogen, CA, USA) was performed to extract total RNA from rabbit lung samples following manufacturer instructions. RNA was processed with the miRcute Plus miRNA First-Strand cDNA Synthesis Kit (Tiangen, Beijing, China) and PrimeScript™ RT reagent kit with gDNA Eraser (Takara, Dalian, China) on the basis of manufacturer recommended. DE-miRNA expression level was validated with an miRcute Plus miRNA qPCR Kit (SYBR Green, Tiangen), and U6 was quantified as the internal control gene. mRNA expression profiles were quantified with RT-qPCR, normalized to the level of GAPDH mRNA (34). Each sample was run in triplicate, and relative miRNA target gene expression was calculated using the 2^−ΔΔCT approach (35). The primer sequences were listed in Table 1.

Rabbit kidney (RK13) cells were grown in DMEM medium (Gibco, USA) containing 10% heat-inactivated fetal bovine serum (FBS, Gibco, USA) and 1% antibiotics. Transient transfections of RK13 cells with the ocu-miR-29-5p mimic or a negative control mimic or inhibitor or negative control inhibitor used lipofectamine 2000 (Invitrogen), and all procedures strictly followed standard instructions. After optimization, a final oligonucleotide concentration of 30 nM per well was adopted. Knockdown of ocu-miR-29-5p expression used an inhibitor or control inhibitor at a final oligonucleotide concentration of 80 nM. The transfection medium was incubated for 6 h, then replaced with complete culture medium and incubated for another 24 h.

RK13 cells were challenged with P. multocida at the multiplicity of infection (MOI) of 10 and harvested at different times (0, 6, 12, and 24 h) for RNA extraction. Uninfected RK13 cells served as negative control, and each trial contained three biological replicates. To determine the effect of ocu-miR-29-5p and inhibitor on P. multocida replication, the RK13 cells (2 × 10⁵ cells/well) were pre-transfected with 30 nM ocu-miR-29-5p mimic or negative control mimic or inhibitor or negative control inhibitor. Three hours after challenge, cells were harvested and washed three times with PBS containing gentamicin (100 ug/mL) for 3 h to kill extracellular P. multocida. Cells were then cultured for 3 h in DMEM containing gentamicin and the cells were lysed using 1% Triton X-100 in PBS for 20 min. Intracellular bacterial loads (CFU/mL) were assessed via plating of lysates onto nutrient agar.

To further confirm the interactions between miR-122 and target genes, the 3′-untranslated region (3′ UTR) sequence of epithelial membrane protein-2 (EMP2) and T-box gene 4 (TBX4) was cloned into the XhoI/NotI restriction sites of the psiCheck2 vector. Recombinant plasmids were verified by DNA sequencing and re-named EMP2 3′ UTR WT and TBX4 3′ UTR WT, respectively. The EMP2 and TBX4 3′ UTR sequence complementary to the ocu-miR-29-5p seed sequence (UUUGGUC) was then mutated to GGGAAGA, and the reconstructed mutated clones were re-named EMP2 3′ UTR MT and TBX4 3′ UTR MT. Luciferase reporter vectors were transfected with either negative control or synthetic ocu-miR-29-5p mimics.

Luciferase activity was measured 48 h after transfection using the Dual-Luciferase Reporter Assay System following the manufacturer's instructions. Luciferase activity were detected using a dual-luciferase reporter assay system (Promega) and a Modulus Single Tube Multimode Reader (Turner Biosystems, Sunnyvale, CA).

All data were analyzed in the GraphPad Prism (version 8.0) software. The results are exhibited as the mean ± SD (standard deviation). Differences between groups were statistically analyzed using Student's t-tests. Thresholds for statistical significance between groups were indicated by ^*P < 0.05; ^**P < 0.01.

To obtain miRNAs related to P. multocida insults, we performed miRNA sequencing by the BGISEQ-500 system. Unreliable sequences were removed from raw reads, and high-quality clean reads were compared with the rabbit reference genome to map miRNA reads (Table 2). Q 30 base percentage was greater than 94%, and clean data exceeded 22.28 M for each sample. Further, the average similarity of samples to the genome was 80%. Pearson's correlation coefficient (R) across samples ranged from 0.984 to 0.997 (Supplementary Table S10), indicating reliable sequencing data for subsequent small-RNA analysis. Finally, clean reads of 18–30 nt for the three group libraries were similar in size distribution and frequency, and most sequences were 21–24 nt; sequences of 22 nt were most abundant (Figure 1). These results confirmed the homogeneity and uniformity of the sequencing data in the nine libraries. Three comparison datasets were set based on time after P. multocida challenge—P0-vs.-P1, P0-vs.-P3, and P3-vs.-P1. Raw reads were deposited to the NCBI database (SRA: SRP315150).

To determine the miRNA sequences in the reads, miRDseep2 software package was applied. The mean miRNA reads accounted for 41.20, 40.02, and 37.78% of clean reads in the P0, P1, and P3 datasets, respectively. Additionally, unknown small RNAs made up 1.35, 1.47, and 1.23% of clean reads (Table 3), and their functions should be studies.

Bioinformatic analysis of sequencing data across all samples yielded 571 miRNAs, consisting of 509 known and 62 new miRNAs (Supplementary Table S1). Notably, expression of DE-miRNA at 0, 1, and 3 dpi following P. multocida challenge showed 19 (10 known and 9 novel), 16 (12 known and 4 novel), and 18 (9 known and 9 novel) specifically expressed in lung tissues (Figure 2A), respectively. These miRNAs may have critical functions for modulating responses to infection at different stages of infection.

We determined FPKM values at three times based on deep sequencing data. miRNA expression profiles are illustrated in Figure 2B. The miRNA expression changed over time, demonstrating that DE-miRNAs play distinct roles from onset to progression of P. multocida infection. Further, the top six highly expressed miRNAs were identified (FPKM > 100000) (Supplementary Table S2): ocu-miR-451-5p, ocu-miR-26a-5p, ocu-miR-143-3p, ocu-let-7a-5p, ocu-miR-22-3p, and ocu-let-7f-5p.

We carried out paired comparative analysis based on the three FPKM datasets, using |log2(Fold Changes)| ≥ 1, Q value ≤ 0.001, and FPKM > 1 as strict selection criteria, and identified 32 DE-miRNAs (Figure 3A, Supplementary Table S3). All DE-miRNAs were previously recognized sequences. Pairwise comparisons among P0, P1, and P3 datasets indicated 13 miRNAs (6 upregulated and 7 downregulated) for the P0-vs.-P1, 15 (13 upregulated and 2 downregulated) for the P0-vs.-P3, and 20 (6 upregulated and 14 downregulated) for the P3-vs.-P1. A Venn diagram was created to depict overlaps among these comparisons. Seven, four, and six time-specific DE-miRNAs were noted, respectively.

Notably, the expression of ocu-miR-107-3p (log2 ^FC = 1.02 and 1.60) and ocu-miR-29b-5p (log2 ^FC = 1.07 and 2.13) in P1 and P3 datasets were significantly increased compared with P0. At P1, miRNA upregulated to the greatest extent was ocu-miR-215-5p (log2 ^FC = 1.47), and miRNA ocu-miR-122-3p (log2 ^FC = −1.71) was the most downregulated. However, miRNAs ocu-miR-21-3p (log2 ^FC = 2.83) and ocu-miR-128a-3p (log2 ^FC = −4.66) replaced these miRNAs, respectively, at P3 (Figures 3B,C). These miRNAs may be involved in the complex regulation of responses to P. multocida infection.

To further clarify the role of DE-miRNAs, we predicted the DE-miRNA's target genes. The 32 DE-miRNAs may have 7,345 non-redundant targets (Supplementary Table S4). DE-miRNAs' target genes are associated with 65 GO terms (26 biological process, 19 cellular component and 20 molecular function) (Figure 4, Supplementary Table S5). The primary enriched biological process, cellular component, and molecular function GO terms were macromolecule metabolic process, intracellular, and binding, respectively. Moreover, a remarkable 130 KEGG pathways were associated with target genes of DE-miRNAs (Figure 5, Supplementary Table S6). Ubiquitin-mediated proteolysis, MAPK signaling, and proteoglycans in cancer were the most notable KEGG pathways (Figure 6). miRNAs have a range of functions and may participate in P. multocida defense by targeting a variety of pathways and genes.

Eight DE-miRNAs were randomly selected, and RT-qPCR data were compared with the deep high-throughput sequencing data (Figure 7). Concordant results in RNA-seq and RT-qPCR analysis confirmed the reliability of sequencing for identifying DE-miRNAs.

Among the DE-miRNAs, ocu-miR-29b-5p was upregulated at both time points and was selected for intensive exploration for a more comprehensive understanding of predicted targets. EMP2 and TBX4 are targets of ocu-miR-29b-5p, and the expression and function were further clarified. Next, changes in the expression of ocu-miR-29b-5p in infected rabbit tissues, including the lung, spleen, liver, and kidney, were quantified.

The results indicated that expression levels of ocu-miR-29b-5p in all four tissues were markedly elevated after P. multocida challenge. Specifically, ocu-miR-29b-5p was rapidly upregulated and reached a peak at 1 dpi in the lung and spleen samples (Figures 8A,B). Expression of ocu-miR-29b-5p in the liver was significantly increased and reached a peak at 3 dpi (Figure 8C). Consistently, expression of ocu-miR-29b-5p in kidney samples was also a time-dependent upregulation (Figure 8D). Moreover, P. multocida challenge stimulated RK13 cells and induced ocu-miR-29b-5p expression, which was significantly upregulated over time, with ocu-miR-29b-5p expression reaching a peak value at 12 h post-infection (Figure 8E).

Due to P. multocida induced enhanced ocu-miR-29b-5p expression, we wondered whether ocu-miR-29b-5p would affect the infectivity of P. multocida. To explicit this doubt, RK13 cells transfected with or without synthetic ocu-miR-29b-5p mimic and inhibitor were challenged with P. multocida. Subsequent P. multocida invasion analysis indicated that bacterial loads in ocu-miR-29b-5p mimic-transfected cells were significantly lower than in control cells, while inhibited expression of ocu-miR-29b-5p promoted the bacterial loads in PK13 cells. Hence, ocu-miR-29b-5p does reduce P. multocida propagation in RK13 cells (Figure 8F).

To clarify the regulatory relationship between ocu-miR-29b-5p and the targets, expression levels of EMP2 and TBX4 in the rabbit lung samples were quantified. EMP2 was markedly decreased in P. multocida-infected rabbit lung, opposite to the expression of ocu-miR-29b-5p (Figure 9A). Similar changes were observed for the relative expression of TBX4 (Figure 9B). Furthermore, the transfection of ocu-miR-29b-5p significantly reduced expression levels of EMP2 and TBX4, whereas knockdown expression of ocu-miR-29b-5p significantly promoted the expression of EMP2 and TBX4 (Figures 9C,D). Following, we confirmed that EMP2 and TBX4 are direct targets of ocu-miR-29b-5p indicated by the vector psiCheck2 linked to wild or mutant 3′ UTRs of rabbit EMP2 and TBX4 genes, respectively. Co-transfection of luciferase constructs with the ocu-miR-29b-5p mimic showed that the ocu-miR-29b-5p mimic markedly reduced luciferase activity from constructs harboring the wild type 3′ UTR of both EMP2 and TBX4, but not the mutant 3′ UTR of these genes (Figures 9E,F). Collectively, ocu-miR-29b-5p negatively regulates the expression of the EMP2 and TBX4 genes through directly specific binding to their 3′ UTR (Figures 9G,H).