Abstract Infection with Pasteurella multocida (P. multocida) causes severe epidemic diseases in rabbits and is responsible for the pronounced economic losses in the livestock industry. Long non-coding RNAs (lncRNAs) have been proven to exert vital functions in regulating the host immune responses to bacterial attacks. However, little is known about how lncRNAs participate in the rabbit's immune response against P. multocida infection in the lungs. LncRNA and mRNA expression profiles were analyzed by transcriptomics and bioinformatics during P. multocida infection. A total of 336 lncRNAs and 7,014 mRNAs were differentially regulated at 1 day and 3 days post infection (dpi). Nearly 80% of the differentially expressed lncRNAs exhibited an increased expression at 3 dpi suggesting that the P. multocida genes are responsible for regulation. Moreover, GO and KEGG enriched analysis indicated that the immune-related pathways including pattern recognition receptors (PRRs), cytokines, and chemokines were significantly enriched at 3 dpi. These results indicate that the dysregulated immune-related genes may play crucial roles in defending against P. multocida attacks. Overall, these results advance our cognition of the role of lncRNAs and mRNAs in modulating the rabbit's innate immune response against P. multocida attacks, which will offer a valuable clue for further studies into exploring P. multocida-related diseases in human. Keywords: long non-coding RNA, gene expression, Pasteurella multocida attacks, high throughput sequencing, functional enrichment, immune response Introduction Pasteurella multocida (P. multocida) is a notorious Gram-negative opportunistic pathogen that is ubiquitous in the respiratory tracts of different animal species and results in enormous economic losses ([31]1, [32]2). It can be sorted into five serogroups (A, B, D, E, and F) and 16 Heddleston serotypes based on its capsular antigens and lipopolysaccharide (LPS) antigens ([33]2, [34]3). In rabbits, P. multocida is a predominant cause of death with a spectrum of hemorrhagic septicemia and respiratory problems ([35]1, [36]4). Notably, the prevalence of P.multocida has been reported as 94%, and even most adult rabbits are considered to be infected with P. multocida ([37]5). Human cases of P. multocida infections, which have been transmitted via licking or biting from rabbits, have raised substantial concern ([38]6, [39]7). However, the molecular pathogenesis of the disease remains unclear. With advances in the next-generation RNA sequencing (RNA-seq), it is possible to quantify both coding and non-coding RNAs (ncRNAs) across a considerably dynamic range, and it is a forceful tool to accurately identify the differential expression (DE) of mRNAs and ncRNAs ([40]8). Advances in RNA-seq have also revealed that < 2% of the mammalian genome is transcribed into protein-coding mRNAs, while the remaining genome produces numerous ncRNAs ([41]9). ncRNAs were deemed as junk molecules of transcription, and it has recently been discovered and proven that they exert an essential function in diverse cellular physiological and pathological processes. Based on their sequence-length cutoff of 200 nucleotides, ncRNAs are subdivided into small ncRNAs, such as microRNAs (miRNAs) and long ncRNAs (lncRNAs). These ncRNAs regulate the gene expression at the levels of pre-transcription, transcription, and post-transcription ([42]10). To date, RNA-seq technology combined with bioinformatics analysis have yielded countless DE mRNAs and ncRNAs, which act important roles in understanding the complex molecular mechanism of host-pathogen interactions ([43]11, [44]12). LncRNAs are defined as transcribed RNA molecules that are more than 200 nucleotides long and lack protein-coding potential. Generally, they can be divided into five distinct subtypes based on the position and direction to nearby protein-coding genes: sense, antisense, intergenic, intronic, and bidirectional ([45]13). LncRNAs are potent regulators that directly interact with other RNAs, proteins, or chromatin to govern gene expression and function through distinct mechanisms, including signaling, scaffolds, decoys, and competitive endogenous RNAs (CeRNA) ([46]14, [47]15). Accumulating shreds of evidence have clarified that these lncRNAs combined with mRNAs contribute to numerous biological processes, especially in the immune response of host-pathogen interactions against various infectious agents ([48]12). For example, the long intergenic ncRNA (lincRNA)-cyclooxygenase 2 (Cox2) was found to trans-regulate both the activation and repression of distinct classes of immune response genes in bone marrow-derived dendritic cells upon toll-like receptor (TLR) ligand stimuli ([49]16). The lincRNA THRIL (TNFα and hnRNPL related immunoregulatory LincRNA) is capable of modulating the expression of the proinflammatory cytokine TNF-α in trans by forming a complex with the hnRNPL (heterogeneous nuclear ribonucleoprotein L) that functions at the TNF-α promoter upon TLR1/2-stimulated THP1-derived human macrophages ([50]17). Meanwhile, multiple TLR ligands induced the lnc MARCKS (myristoylated alanine-rich protein kinase C) production to the master regulator of inflammatory responses interacted with APEX1 (apurinic/apyrimidinic endodeoxyribonuclease 1) to form a hnRNPL complex at the MARCKS promoter ([51]18). Recently, an increasing number of lncRNAs have been identified, and their functions in the host-pathogen interactions of humans, mice, and pigs have been clarified ([52]12). However, the alteration in the lncRNA expression induced by rabbit P. multocida remains unclear, as do the roles of these transcripts in modulating the lung response to infection remain unclear. In our research, we depicted the lncRNA and mRNA profiles of rabbit lungs infected with P. multocida for 1 day (P1) and 3 days (P3) or without P. multocida (control, P0). This study aimed to obtain the DE transcripts (lncRNAs and mRNAs) during the P. multocida insults with the aim of exploring the primary biological processes and pathways associated with the pathogenesis mechanisms. Thus, the data obtained from this research could be used to establish different cellular and molecular events upon the P. multocida attack, which could contribute to uncovering the molecular pathogenesis implying the P. multocida attack. The findings of this study also provide a valuable reference for further study on lncRNAs concerning P. multocida-related diseases in human. Materials and Methods Pasteurella multocida Isolation, Identification, Pathogenicity, and Infection The strain was isolated from a suspected rabbit lung infected with P. multocida and stored in our lab. Then, the recovered strain was inoculated onto 5% defibrinated sheep blood agar and MacConkey agar and incubated at 37°C for 24 h. A presumptive P. multocida colony was preliminarily identified using routine biochemical testing and API 20 E strips (bioMérieux, Durham, NC, USA), and further confirmed by PCR assay to amplify the inherent gene according to the published articles ([53]19). All animal studies in this article were performed according to the guidelines of the Institutional Animal Care and Use Committee of Shandong Agricultural University and the “Guidelines for Experimental Animals” of the Ministry of Science and Technology (Beijing, China). New Zealand rabbits that were 38-days-old were purchased from Qingdao Kangda Rabbit Industry Development Co., Ltd. They were tested negative for P. multocida by separating them from the nares. The detected samples were collected by rotating a sterile cotton swab three times in both anterior nares. We also tested negative for P. multocida by detecting the antibodies against P. multocida in their serum applying a commercial ELISA assay. Blood samples were obtained from the marginal vein of the rabbit ear and serum was separated from the blood using a centrifuge (1,000 rpm, 10 min). All the rabbits were adapted to the experimental environment for seven days, and then the experiment was conducted. The bacteria were quantified using a plate counting method and diluted in PBS at a final concentration of 1 × 10^9 colony forming units (CFU)/mL. Eight 45-days-old New Zealand rabbits that were confirmed free of P. multocida were divided randomly into two groups of four individuals. In the infected group, the New Zealand rabbits were intraperitoneally injected with 1 mL of the prepared bacteria. In the control group, the New Zealand rabbits were intraperitoneally injected with 1 mL of PBS. The incidence and mortality of the New Zealand rabbits were observed, and lung samples were collected from the New Zealand rabbits for bacterial isolation and identification. The cultured P. multocida was diluted using a 10-fold gradient with sterile saline at the final centrations of 10^9, 10^8, 10^7, 10^6, 10^5, and 10^4 CFU/mL. Fifty-six New Zealand rabbits were randomly divided into seven groups, with eight individuals in each group and intraperitoneally injected with different concentrations of P. multocida solution and 1 mL of PBS. Morbidity and mortality were observed for 15 consecutive days. The median lethal dose (LD 50) was calculated using the SPSS 19.0 software. Twenty-five rabbits were randomly divided into the infected (n = 20) and control (n = 5) groups. Prior to the experimentation, the control group rabbits were euthanized by pentobarbital overdose (100 mg/kg, intravenous). The lung samples were collected from all freshly killed rabbits. The infected group was subcutaneously injected with 10^7 CFU of P. multocida in 1 mL of PBS. After P1 and P3, every five rabbits with significant clinical symptoms of depression, anorexia, snuffles, serous nasal exudate, or dyspnea were euthanized, and the lung of each rabbit was dissected. All samples were collected aseptically, frozen in liquid nitrogen immediately, and transported to the lab, then stored at −80°C. The remaining alive rabbits, post-experiment, were euthanized by the same procedure. Total RNA Isolation and Qualification 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 the rabbits' lung tissues through a specialized Total RNA kit (Tiangen Biotech Co, Beijing, China) following the manufacturer's protocol. RNA concentration and purity were monitored using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA), and RNA integrity was checked using the Agilent 2100 Bioanalyzer (Thermo Fisher Scientific, MA, USA). The samples with RNA Integrity Number (RIN) ≥8 were subjected to further analysis. Library Construction and Sequencing For the mRNA and lncRNA library, we firstly removed the ribosomal RNA (rRNA) using target-specific oligos, and we removed both the cytoplasmic and mitochondrial ribosomal RNA using RNase H reagents from the above total RNA samples. Secondly, the obtained RNA was fragmented into small pieces using divalent cations under elevated temperature to achieve a solid-phase reversible immobilization bead purification. Utilizing the cleaved RNA fragments, we synthesized the first-strand cDNA using random hexamer primers and reverse transcriptase, followed by the second-strand cDNA synthesis using DNA Polymerase I and RNase H. In this reaction, the RNA template was removed, and an alternative strand was synthesized. Moreover, dTTP was replaced by dUTP to generate the ds cDNA. These cDNA fragments had the addition of a single “A” base and subsequent ligation of the adapter. After the uracil-DNA glycosylase treatment, the incorporation of dUTP quenched the second strand during amplification. The products are enriched with PCR to create the final cDNA library. Finally, the quality of the constructed libraries was assessed by checking the distribution of the fragment size using the Agilent 2100 Bioanalyzer, and their quantity was assessed using quantitative real-time PCR (qRT-PCR) (TaqMan Probe). The qualified libraries were sequenced on the BGISEQ-500 System (BGI-Shenzhen, China). Expression and Differential Expression Analysis Raw reads were processed by removing the low-quality reads and adaptor contaminants, and other poly-N using the SOAPnuke (v1.5.2) ([54]20), and the resulting clean reads were used for further study. The clean reads were aligned to the reference genome of Oryctolagus cuniculus 2.0 using both the HISAT2 (v2.0.4) ([55]21) and Bowtie2 (v2.2.5) ([56]22). The fragments per kilobase of transcript per million fragments mapped (FPKM) was used to quantify the expression levels of a gene or lncRNA calculated by the RSEM (v1.2.12) software ([57]23). For the samples with biological replicates, DE analysis of lncRNAs and mRNAs in the three groups was determined by the DESeq2 ([58]24) with a Q ≤ 0.05. The DESeq2 assumes a negative binomial distribution for gene counts, normalizes for read depths, and fits a generalized linear model. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Enrichment Analyses Enrichment analysis of the DE genes was conducted using the Gene Ontology (GO, [59]http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, [60]https://www.kegg.jp/) databases by Phyper ([61]https://en.wikipedia.org/wiki/Hypergeometric_distribution) based on a hypergeometric test. The significant levels of terms and pathways were corrected using a Q value with a rigorous threshold (Q value ≤ 0.05) by Bonferroni ([62]25). The GO terms or KEGG terms meeting this condition were considered as significantly enriched terms. qRT-PCR Validation Total RNA was separated from the lung tissues using the Trizol Reagent (Takara Biotechnology, Dalian, China), and the cDNA was reversed from the RNA using a two-step qRT-PCR Kit (Takara Biotechnology, Dalian, China) following the manufacturer's instructions. The qRT-PCR reaction was conducted using the SYBR Green assay (Takara Biotechnology, Dalian, China) in the CFX96 Real-Time PCR Detection System (Bio-Rad) as previously described ([63]11). The specific quantitative primers of lncRNAs and mRNAs used in this study were either designed using the primer 6.0 software or cited from previously published literature ([64]Table 1). All the primers were synthesized by Sangon Biotech (Sangon Biotech, Shanghai, China). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an endogenous control for the lncRNA and mRNA. Table 1. The primers used in this experiment for qRT-PCR. Gene name Sequence (5^′-3^′) Production size (bp) References