Application of multi-omics technology in pathogen identification and resistance gene screening of sheep pneumonia

Abstract

Background: Pneumonia constitutes a major health challenge in sheep, severely compromising growth rates and overall productivity, and resulting in considerable economic losses to the sheep industry. To address this issue, the development of disease-resistant breeding programs based on the identification of genetic markers associated with pneumonia susceptibility is of critical importance. This study investigated a sheep population on a farm where pneumonia was endemic. The purpose was to use multi-omics methods to rapidly identify the principal pathogens responsible for pneumonia outbreaks, and to screen for genetic loci and key genes related to pneumonia resistance, thereby providing a scientific basis for the implementation of targeted breeding strategies for pneumonia resistance.

Results: Here, we assessed the impact of pneumonia on sheep growth by evaluating the pneumonia phenotypes of 912 sheep. High-throughput transcriptome sequencing of 40 lungs was conducted to obtain exogenous RNA fragments for microbial sequence alignment. Additionally, 16S rRNA sequencing was performed on lung tissues from 10 healthy and 10 diseased sheep to identify biomarkers associated with phenotypic differences. Mycoplasma ovipneumoniae was identified as the primary pneumonia pathogen, and its presence was further validated by load quantification and immunohistochemical analysis. Integration of genome-wide association study (GWAS) data from 266 lung pathological scores with transcriptome-based differentially expressed genes analysis enabled the identification of five single nucleotide polymorphisms (SNPs) and three potential candidate genes associated with Mycoplasma pneumonia. Subsequent genotyping and phenotype association analyses confirmed the significance of two SNPs and established a strong association between the FOXF1 gene and resistance to Mycoplasma pneumonia.

Conclusions: High-throughput sequencing technologies have enabled the rapid and accurate identification of the causative pathogen of sheep pneumonia. By integrating multi-omics data, two genomic loci significantly associated with Mycoplasma pneumonia were screened, as well as an anti-Mycoplasma pneumonia key gene, FOXF1.

Keywords: Mycoplasma ovipneumoniae; Disease-resistant breeding; Genome-wide association study; Multi-omics; Pathogen identification; Sheep.

Fig. 1

Phenotype identification and histopathological observation of sheep pneumonia lesions. (a) Representative gross photographs of lungs from Hu sheep in the healthy (up) and diseased (down) lung phenotype groups. The lungs from the healthy group were elastic and presented a homogeneous pink appearance, whereas those from the diseased group demonstrated severe inflammatory lesions, including extensive red and gray hepatization (indicated using black arrows). (b) Representative lung tissue hematoxylin–eosin staining images (magnification: 4×, 10×, and 40×) of healthy and diseased groups. (c) Example of pathological scoring of lung tissue (hematoxylin–eosin; magnification: 4×). Lesions are scored as follows: +1, minor; +2, mild; +3, moderate; +4, severe; and + 5, extremely severe. In the picture, red triangle denotes partial alveoli retain their original structure; purple red pentagram denotes partial bronchial lumens contain large numbers of neutrophils, lymphocytes (+ 1); green arrow denotes partial peribronchial lymphoid follicular hyperplasia (+ 1); blue arrow denotes partial alveolar compensatory dilation (+ 2); orange arrow denotes inflammatory cell infiltration at the junction of lung interstitium and lung parenchyma (+ 2); yellow arrow denotes partial collapse and substantial lesions of lung tissue, mainly filled with monocytes (+ 2)

Fig. 2

High-throughput sequencing analysis–based identification of pathogenic microorganisms involved in sheep pneumonia. (a) Differential analysis of lung microbial expression between the healthy and diseased groups. (b, c) Linear regression between expression of candidate pathogenic microorganisms M. ovipneumoniae (M.O., b) and P. multocida (P.M., c) and pulmonary pathological scores. (d) Mapping of sequencing results on the M. ovipneumoniae reference genome. (e) Phylogenetic tree constructed using LDH of M. ovipneumoniae obtained through sequencing. Red triangle indicates M. ovipneumoniae strain in this study. ** p < 0.01, ns indicates no significant difference between groups, unpaired Student’s t test

Fig. 3

16S rRNA sequencing confirms M. ovipneumoniaeas a biomarker between healthy and diseased lung groups. (a) Microbial community α-diversity between the healthy and diseased groups, including the Chao1, Simpson, Pielou’s evenness, Faith’s PD, and Good’s coverage indexes. (b, c) PCoA (b) and NMDS (c) of the ASVs level based on Bray–Curtis distances between the healthy and diseased groups. (d) Bar chart of microbial composition analysis at the phylum level (relative abundance > 1%). The top 10 microbial taxa in terms of abundance are displayed by different color scales. Abundance in a group is the average sample abundance in the group. (e) Bar chart of microbial composition analysis at the genus level (relative abundance > 1%). The top 10 microbial taxa in terms of abundance are displayed by different color scales. Abundance in a group is the average sample abundance in the group. (f) Heatmap of the distribution and variability of bacteria in the lungs. (g) LEfSe analysis. The significance threshold is set to an LDA score of 4

Fig. 4

Verification of M. ovipneumoniae load in experimental sheep lung tissue. (a) Detection of M. ovipneumoniae load between groups in an expanded sheep population (healthy = 86, diseased = 184). LDH gene fluorescence quantitative detection was used to characterize the load of M. ovipneumoniae, with β-tubulin gene as an endogenous control. (b) Immunohistochemical detection of M. ovipneumoniae in the lungs of diseased group. A negative control represents the absence of incubation of M. ovipneumoniae primary antibody, while a positive indicates normal incubation of M. ovipneumoniae primary antibody. The yellow color in the field of view of the positive picture indicates positive staining of M. ovipneumoniae, and the darker the staining, the stronger the positive intensity. The data are presented as means ± standard errors of the means (n = 3). ** p < 0.01, unpaired Student’s t test. LDH: lactate dehydrogenase

Fig. 5

GWAS of pulmonary pathological scores. (a) Manhattan plot of the GWAS results, with the threshold line (gray line) set at 6. The genes annotated in the figure are those annotated within 200 kb upstream and downstream of the significant loci on the threshold line. The genes denoted in red and green are those significantly upregulated and downregulated in the diseased lung group in differentially expressed gene analysis, respectively. (b) Q-Q plot of the GWAS result. The red line represents the expected statistical distribution, the blue dot represents the actual observed statistical distribution, and the light blue area represents the confidence interval

Fig. 6

Transcriptomic analysis for candidate key genes associated with M. ovipneumoniae infection. (a) PCA of RNA-Seq data in the healthy and diseased groups. (b) A volcano plot presenting DEGs between the healthy and diseased groups. (c) Heatmap showing bidirectional clustering of DEGs and samples. (d) GO enrichment analysis of DEGs. GO classification is based on cellular components (CC), molecular functions (MF), and biological processes (BP). The top 10 GO term entries with the smallest p value, i.e., the most significant enrichment, in each GO classification were selected for display. (e) KEGG pathway enrichment analysis of DEGs. The larger the bubble, the greater the degree of enrichment. Closer FDR values to zero indicate more significant enrichment. The top 20 KEGG pathways with the smallest FDR values, i.e., the most significant enrichment, were selected for display. (f) Verification of RNA-Seq results through qRT-PCR. Relative expression levels calculated from standard curves were normalized to the β-actin gene, as the endogenous control. The data are presented as means ± standard errors of the means (n = 3). * p < 0.05, ** p < 0.01, unpaired Student’s t test. PROM1: prominin 1; PAK5: p21 (RAC1) activated kinase 5; SLIT1: slit guidance ligand 1; PGAM1: phosphoglycerate mutase 1; UBTD1: ubiquitin domain containing 1; ADGRB3: adhesion G protein-coupled receptor B3; FOXF1: forkhead box F1; PRDX6: peroxiredoxin 6; PTGS1: prostaglandin-endoperoxide synthase 1; SMAD7: SMAD family member 7

Fig. 7

Different genotypes of important SNPs affect FOXF1 gene expression and individual pneumonia phenotype. (a) The expression level of candidate FOXF1 gene and M. ovipneumoniae load in the lungs of individuals with different genotypes (FOXF1: AA = 116, AG = 78, GG = 9; LDH: AA = 114, AG = 74, GG = 9) of SNP loci (chr.14: 13215150). (b) The expression level of candidate FOXF1 gene and M. ovipneumoniae load in the lungs of individuals with different genotypes (FOXF1: AA = 5, AG = 51, GG = 138; LDH: AA = 5, AG = 50, GG = 134) of SNP loci (chr.14: 13215157). (c) The impact of different genotypes of important SNPs (chr.14: 13215150, AA = 116, AG = 78, GG = 9; chr.14: 13215157, AA = 5, AG = 51, GG = 138) on their pulmonary pathological scores. Relative expression levels of FOXF1 calculated from standard curves were normalized to the β-actin gene, as the endogenous control. LDH gene fluorescence quantitative detection was used to characterize the load of M. ovipneumoniae, with β-tubulin gene as an endogenous control. The data are presented as means ± standard errors of the means (n = 3). Different superscript lowercase and capital letters respectively showed significant difference (P < 0.05) and extremely significant difference (P < 0.01), one-way ANOVA and Dunnett’s test. * p < 0.05, unpaired Student’s t test. FOXF1: forkhead box F1; LDH: lactate dehydrogenase