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.

Background

Animal diseases are one of the main threats to modern livestock and poultry production. According to reports, in developed countries, economic losses caused by diseases account for about 17% of the total output value of animal husbandry, while in developing countries, this proportion reaches 35–50% [1]. Sheep are economically important livestock species globally. Pneumonia is a respiratory disease that is widely prevalent in sheep farms. The incidence rate of pneumonia in ruminant animals ranges from 10 to 40%, with higher percentages in childhood [2]. Pneumonia can cause growth retardation and increase the risk of secondary pleurisy in sheep, resulting in major production losses [3]. Bacteria-induced infectious pneumonia is most common, and it has been confirmed that pathogens such as Mycoplasma ovipneumoniae, Mycoplasma filamentosa, Pasteurella multocida, Cryptobacterium septicum, Mannheimia haemolytica, Klebsiella pneumoniae, Escherichia coli, Streptococcus, and other pathogens can all cause pneumonia individually or in combination [4,5,6].

Next generation sequencing (NGS) technology is a promising method for pathogen identification, including rare and newly identified viruses [7]. NGS is independent of cultivation and can efficiently and unbiasedly detect pathogenic microorganisms in clinical samples, providing comprehensive analyses [8, 9]. Previous studies have generally focused on pulmonary microbiota composition in healthy sheep [10, 11]; while the data on the pulmonary microbiota related to diseases in sheep remains relatively scarce [6]. In microbiome research, comparing 16S rRNA gene sequences has become a reliable approach for identifying microorganisms associated with pathogenicity and infection [12]. Compared to pathogen isolation and identification in veterinary medicine, comprehensive analysis of multi-omics can help to quickly and accurately identify pathogens and discover pathogens that are difficult to cultivate.

Mycoplasma ovipneumoniae (M. ovipneumoniae) colonization has been detected in approximately 90% of pneumonia-associated clinical surgery cases [13]. Moreover, it has been identified as the main pathogen causing severe epidemic pneumonia in bighorn sheep [14]. M. ovipneumoniae has also caused severe respiratory disease outbreaks in goats across multiple regions [15, 16]. With the increasing proportion of large-scale breeding in recent years, Mycoplasma pneumonia has become the highest incidence rate infectious disease in the sheep industry. The main manifestations of M. ovipneumoniae infections include wheezing, coughing, fever, weight loss, and pulmonary interstitial hyperplasia inflammation [17]. It is characterized by an unusually wide host range, high phenotype, strong infectivity and high incidence rate, and is prone to secondary or mixed infection of other bacteria [18]. Currently, antibiotics remain the main means of treating M. ovipneumoniae infections. However, the long-term use of antibiotics can affect meat quality, produce antibiotic residues, and even lead to the development of drug resistance. Therefore, starting from the genetic basis and screening for resistance genes of Mycoplasma pneumoniae for disease-resistant breeding is an important way to fundamentally address this problem.

In farm animal breeding programs, disease resistance and susceptibility—the natural, unique mechanisms by which a host responds to infectious pathogens—are often used as main references. Disease resistance depend on the host’s biological and immune responses, so livestock’s ability to resist multiple diseases varies among individuals [19, 20]. Studies on sheep and cattle have reported potential genetic components involved in the differences in pneumonia susceptibility observed among animals [3, 21]. Heritable variations may aid in selecting animals that can resist and recover from infections of multiple pathogenic factors [22]. In addition, omics techniques, bioinformatics, and Genome-wide association studies (GWAS) can identify genomic regions associated with potential causal mutations that affect susceptibility to infectious diseases [23,24,25], and the discovery of these regions may also lead to the establishment of new disease diagnostic tools and alternative therapies. Therefore, using molecular biology and molecular genetics techniques to search for disease-resistance genes for disease-resistance breeding and improve genetic disease resistance in livestock has become an effective supplement to conventional disease prevention and control methods. The screening and cultivation of disease-resistant animal strains have also become a major research direction.

In the current study, we: (1) showed the severe impact of pneumonia on the growth and productivity of sheep, (2) integrated multi-omics technologies to quickly and accurately identify the pathogen of sheep pneumonia, M. ovipneumoniae, and (3) identified genomic regions associated with Mycoplasma pneumonia. In particular, by combining GWAS, transcriptomics, and association analysis with pneumonia, two effective SNPs related to Mycoplasma pneumonia and a candidate FOXF1 gene for anti-Mycoplasma pneumonia were screened and validated.

Methods

Animal parameters and sample collection

All the sheep in this study were from the Minqin experimental farm of Lanzhou University (N38°43′41′′, E103°013′), and we first evaluated their respiratory disease status. At the age of 180 days and after fasting for 12 h, whole blood for measurement of blood physiological indicators, DNA extraction and whole genome sequencing (WGS) was obtained from each sheep via the jugular vein. Subsequently, on the premise of ensuring animal welfare and reducing animal suffering, each animal was slaughtered according to standard commercial procedures following the requirements of the China Council on Animal Care. This study received approval from the Lanzhou University Research Ethics Committee (approval number: 2021-02). During two rounds of slaughter, the lungs of 912 Hu sheep (905 male and 7 female) from this studied sheep farm were examined grossly and photographed. After identifying the lung phenotype of each Hu sheep, lung tissue block samples were collected in triplicate from normal areas of healthy lungs and diseased areas of inflammatory diseased lungs, each for subsequent sequencing, slice staining, and molecular experiments.

All included sheep maintain similar birth ages in each slaughter (i.e. each batch, as detailed in Table 1). All batches of sheep were raised in single pens in identical environments, with ad libitum access to food and water and no differences in management methods or nutritional levels.

Table 1 Incidence of pulmonary lesions caused by Mycoplasma pneumonia in sheep population

Clinical examination, histopathological observation, and phenotype data collection

The lungs health or pathological conditions, including the lung lesions and tactile perception of lung elasticity, were comprehensively evaluated. Next, we unfolded the ventral sides of the lungs horizontally and photographed them. The lung tissue blocks used for slice staining were fixed in 4% formaldehyde, embedded in paraffin according to standard procedures, then sectioned, and stained with hematoxylin and eosin for histopathological examination under an optical microscope. Moreover, we assessed the pathological changes in lung tissue and scored them on a 5-point system, ranging from 1 to 5, indicating increasing lesion severity. Additional file 1 summarizes the general histopathological observations and scores of the lung slices.

Combining clinical dissection and histopathological observations, we assessed the lung lesion phenotype and collected pathological score data. Based on these data, we divided 912 Hu sheep into healthy and diseased groups for subsequent analyses.

Collection and analysis of growth traits and blood parameters

The growth performance measurement of Hu sheep starts from 80 days and ends at 180 days, during which the body weight (BW) and feed intake (FI) are measured every 20 days. The BW of lambs at 80 days is recorded as the initial weight, and the BW at 180 days is recorded as the final weight. At the same time, the live weight at slaughter and carcass weight were recorded. In addition, calculate the Total FI of each Hu sheep throughout the entire experimental period (80 d-180 days) and based on the corresponding data, calculate the average daily feed intake (ADFI), average daily weight gain (ADG), and feed conversion rate (FCR): ADFI (kg/day) = Total FI (kg) / 100 (day); ADG (kg/day) = (final weight (kg) - initial weight (kg)) / 100 (day); FCR = ADFI (kg/day) / ADG (kg/day). Sheep whole blood samples collected through the jugular vein before slaughtered were measured on the fully automatic five-classification animal blood cell analyzer Aldex ProCyto Dx* for their blood physiological indicators.

Based on the phenotype grouping of the lungs mentioned above, the growth performance and blood parameter data of all sheep were combined and analyzed for inter group differences to evaluate the impact of pneumonia related respiratory diseases on sheep growth performance and blood parameters.

Pathogen identification of sheep pneumonia by high-throughput transcriptome sequencing

Among the 912 Hu sheep with lung phenotype data, 20 Hu sheep in the healthy group (10 each from batches 1 and 2) and 20 Hu sheep in the diseased group (10 each from batches 1 and 2) were random selected to collect samples from normal or typical lesion areas of the lungs (n = 40) for total RNA extraction. After RNA quality testing, the kit was used to remove ribosomal RNA (rRNA) and retain mRNA and other non-coding RNA. cDNA was synthesized by fragmenting and reverse transcription of RNA, and then a series of steps such as end repair and linker ligation were performed to construct the lncRNA sequencing library. Finally, the library was sequenced on the sequencing platform Illumina HiSeq (Illumina, San Diego, CA, USA) to obtain 150 bp paired-end (PE) reads.

We used Fastp to perform quality control of the offline data based on our lncRNA sequencing results and Picard to convert the results into BAM format files. For pathogen analysis, Genomic Analysis Toolkit (GATK) PathSeq was used according to the gatk-best-practices process. The tables of 40 individual pathogen analyses were merged, and the sum of the unambiguous sequences was obtained. We then screened microbial species with a high number of unambiguous sequences for between-group difference analysis. Bowtie2 was used to compare sequences with the reference genome of the aforementioned microbial species, with a filtering condition of mq ≥ 30, and the mapped sequences were counted. Next, we calculated reads mapped per million reads (RPM) for each microorganism and analyzed the between-group differences. Linear regression analysis was performed on the sequence reads of pathogenic microorganisms and tissue pathological scores. Coverage statistics were performed using Bamdst for the sequenced sequence comparison on the reference genome of M. ovipneumoniae. Finally, the M. ovipneumoniae lactate dehydrogenase (LDH) sequence was extracted and subjected to multiple sequence alignment using MAFFT. Then, a phylogenetic tree was constructed using MEGA.

16S rRNA sequencing analysis for microbial community diversity in lung tissues

From 40 lung tissue samples used for transcriptome sequencing, 20 lung tissue samples (10 each from the healthy and diseased groups) were random selected for 16S rRNA sequencing to determine the colonized microbial community.

Total DNA was extracted using the OMEGA Soil DNA Kit (M5635-02; Omega Bio-Tek, Norcross, GA, USA). Next, the V3-V4 region of the standard bacterial 16S rRNA gene was amplified through polymerase chain reaction (PCR) with the following primers: forward, 5′-ACTCCTACGGGAGGCAGCA-3′, and reverse, 5′-GGACTACHVGGGTWTCTAAT-3′. PCR reactions, containing 5 µL 5Xreaction buffer, 5 µL 5XGC buffer, 2 µL dNTP (2.5mM), 1 µL each primer (10µM), 2 µL DNA template, 8.75 µL ddH₂O and 0.25 µL Q5 DNA Polymerase (NEB) in a volume of 25 µL, were amplified by thermocycling: 2 min at 98 °C for initial denaturation; 30 cycles of 15 s denaturation at 98 °C, 30 s annealing at 55 °C, and 30 s extension at 72 °C; followed by 5 min final elongation at 72 °C. The amplification products were purified and recovered to prepare a sequencing library. Next, we performed dual-end sequencing of the library on Illumina NovaSeq. The PE reads obtained through sequencing were first spliced according to their overlap relationships, using FastQC for sequence quality control. The filtered data were further processed using the DADA2 method [26] on QIIME2 (https://qiime2.org) to generate amplicon sequence variants (ASVs). Finally, we obtained the abundance tables for these ASVs (i.e., the ASV feature tables). After completing the denoising of all libraries, we merged all ASV feature tables and removed singleton ASVs. The Greengenes (http://greengenes.secondgenome.com/) database was selected for ASVs taxonomic annotation. After distinguishing and grouping the samples, we conducted species composition analysis at the phylum and genus levels, and performed a significant difference analysis of the abundance of the Mycoplasma and Pasteurella between groups at the genus level. α-Diversity metrics, including the Chao1, Simpson, Faith’s PD, Pielou’s evenness, and Good’s coverage indexes, were calculated based on the ASV feature tables by using the diversity function in the R package vegan (https://rdrr.io/cran/vegan/). β-Diversity was estimated using the Bray–Curtis distance algorithm and visualized using Principal coordinate analysis (PCoA) and Nonmetric multidimensional scaling (NMDS). Finally, we used the Python software package Linear discriminant analysis (LDA) effect size (LEfSe) to analyze robust biomarkers with LDA scores > 4 and p < 0.05 in both groups.

Verification of M. ovipneumoniae load in experimental sheep lung tissue

Due to the high coverage on the reference genome of M. ovipneumoniae, the expression of LDH gene was used to characterize the load of M. ovipneumoniae. The fluorescence quantitative detection of LDH gene was designed and the detailed information of the relevant primers was listed in Additional file 2. Using lung DNA as a template, fluorescence quantitative PCR (qPCR) was performed on LDH and standardized with β-tubulin to characterize the M. ovipneumoniae load in each sheep lung. qPCR reaction system (10 µL): 0.5 µL each primer (10 µM), 1 µL DNA template, 5 µL TB Green Premix Ex Tag II (Takara), 3 µL ddH₂O. Quantitative program: pre-denaturation at 95 °C for 3 min, denaturation at 95 °C for 10 s, annealing at 60 °C for 10 s, extension at 72 °C for 15 s, 40 cycles. All quantitative measurements were obtained in triplicate, and the qPCR reactions were performed in a Bio-Rad CFX384 sequence detector.

Immunohistochemistry of M. ovipneumoniae

In order to further characterize M. ovipneumoniae in lungs, the lung tissue paraffin sections prepared were sequentially dewaxed, rehydrated, antigen repaired, and quenched with endogenous peroxidase. Then, sections were blocked with 5% BSA for 40 min and incubated with primary antibody (self-made rabbit immune serum, dilution ratio 1:500) and goat anti-rabbit HRP labeled secondary antibody. Simultaneously set up a negative control that was not incubated with the antibody. DAB staining and hematoxylin counterstaining of sections were performed, and finally dehydrated and sealed. The images were examined under an optical microscope and collected.

WGS and data processing

Genomic DNA was extracted from Hu sheep blood samples according to the instructions of the Easy Pure Blood Genomic DNA Kit (TransGen Biotech) and sent to Molbreeding Biotech (Shijiazhuang, China). A PE sequencing library with an insert size of 350 bp was constructed using a Truseq Nano DNA HT Sample Preparation Kit (Illumina), according to the manufacturer’s instructions. The libraries that have passed quality inspection were subjected to PE sequencing with 150-bp reads (PE150) on Illumina HiSeq. Raw sequencing data were converted to fastq files [27]. Here, FastQC (version 0.10.1) was used to filter the raw data; in particular, linker sequences, sequences containing the uncalled base N, and < 60-bp sequences were discarded.

The obtained clean reads were mapped to the sheep reference genome Oar_v1.0 using Burrows-Wheeler Aligner (version 0.7.8) [28]. SAMtools (version 1.12) [29] was used to convert the mapping results to the BAM format, remove duplicate reads, and retain pairs with the highest mapping quality. Then, GATK (version 3.4.0) was used to call SNPs and generate files in the VCF format. Next, we used VCFtools (version 0.1.14) software to correct the GATK results and screen out high-quality SNPs [30]. The screening conditions were as follows: (1) number of supported SNP reads > 5, (2) deletion rate < 20%, (3) minimum allele frequency > 5%, and (4) root mean squared mapping quality > 20. Beagle (version 5.0), which self-fills genotype data, was used to eliminate Not Available (NA) value present in the genotype data. Finally, 23,480,428 SNPs and 912 sheep, which passed the filtering and quality control steps, were used for subsequent analysis.

GWAS and significant SNPs

A GWAS was performed on lung pathological scores using the General Linear Model in the R package rMVP [31]. In addition, the influence of population stratification was corrected through the addition of the first five PCs derived from whole-genome SNPs. Lamb birthplace and BW variables were included as covariates in the analysis model. To avoid potential false positives, the genome-wide significance threshold was adjusted, and the significance threshold of GWAS results for sheep lung pathological scores traits was set to p < 1 × 10^− 6. The R package CMplot (https://github.com/YinLiLin/R-CMplot) was used to draw quantile–quantile (Q-Q) and Manhattan plots, set the parameter ‘threshold = 10^− 6’, which means the gray threshold line in the graph is 6. Before gene annotation, we used PopLDecay to assess for the attenuation of linkage imbalance. Ensembl and NCBI were used to determine the location of the significant SNPs within the sheep reference genome Oar_v1.0 and identify annotated genes within 200 kb upstream and downstream of significant SNPs.

Conventional analysis of RNA-Seq data and screening of candidate key genes

The high-quality sequences (i.e., clean data) obtained by filtering the raw data of transcriptome sequencing were mapped to the sheep reference genome Oar_v1.0, and the mapped reads were assembled and spliced to restore the transcriptional sequence. We then used Cufflinks [32] to calculate each gene’s expression level, standardized to fragments per kilobase of exon per million fragments mapped (FPKM). Next, we used the R package DESeq to perform principal component analysis (PCA) on each sample based on gene expression levels. The stattest function in the Ballgown package was used to search for mRNAs differentially expressed between the healthy and diseased groups, based on the following|log₂(fold change)| > 1 and adjusted p < 0.05. Next, further expression clustering and enrichment analyses, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses, were conducted on the identified differentially expressed genes (DEGs). To validate our RNA-Seq results, total RNA of sheep lung tissues extracted from those performed RNA-Seq were subjected to cDNA synthesis using the RevertAid First Strand cDNA synthesis kit (Thermo Scientific Co., Ltd, China). Next, we performed qRT-PCR to measure the expression levels of 10 randomly selected DEGs (PROM1, PAK5, SLIT1, PGAM1, UBTD1, ADGRB3, FOXF1, PRDX6, PTGS1 and SMAD7) by using a TB Green Kit (Takara) and lung tissue cDNA. Additional file 2 presents the details of our qRT-PCR primers. The expression level of each gene was normalized by β-actin [33]. The reaction system and program of qRT-PCR are the same as those of qPCR mentioned above. All quantitative measurements were obtained in triplicate, and all qRT-PCR assays were performed on the Bio-Rad CFX384 sequence detector. DEGs obtained from transcriptome analysis were cross-analyzed with annotated genes from the above GWAS results to screen for potential candidate functional genes and their associated significant SNPs.

Genotyping and association analysis of significant SNP loci

We randomly selected 854 blood DNA samples from a population of 912 sheep, and used competitive allele-specific-based polymerase chain reaction (KASPar) (LGC Genomics, Hoddesdon, UK) to genotype the significant SNPs, according to a previously published method [34]. The primer pairs designed for our KASPar genotyping are listed in Additional file 2.

To verify the relationship between significant SNPs and screened candidate key genes, as well as their correlation with the anti-Mycoplasma pneumonia mechanism. We grouped the expanded sheep population (n = 203/197/194/189) based on the different genotypes of each significant SNP locus, and then detected the differences in the expression of corresponding candidate key genes, as well as in the lung phenotypes (including MO load and lung pathological score) between groups. The expression level detection of candidate key genes was evaluated by qRT-PCR, as for DEGs. The detection method for M. ovipneumoniae load was qPCR. Furthermore, we studied the association between significant SNP loci and inflammation related blood physiological indicators based on genotype grouping.

Statistical analysis

qRT-PCR and qPCR data were processed using the 2^−ΔΔCt method [35]. Internal reference genes were used to standardize the expression of each gene. Differences between two groups (healthy and diseased groups) were evaluated using Student’s t test, whereas those between three groups (three genotypes groups) were assessed one-way analysis of variance followed by Dunnett’s test for multiple comparisons. These tests were performed using SPSS (version 20.0; IBM, Chicago, IL, USA). Graphics were drawn using R (version 4.2.3) and GraphPad Prism (version 5.0; GraphPad Software, San Diego, CA, USA). All results are expressed as means ± standard errors of the means. Finally, p < 0.05 and p < 0.01 were considered to indicate significance and extreme significance, respectively.

Results

Phenotypic identification of pneumonia lesions and their impact on sheep growth performance

During an investigation in early July 2022, we found that approximately 20% of sheep experienced coughing, runny nose, and elevated body temperature at our study sheep farm in Minqin County, Wuwei City, Gansu Province, China. In some of these sheep, the clinical respiratory symptoms worsened over time. Two batches of sheep fed in this farm were slaughtered. In clinical autopsy, we observed both healthy lungs and diseased lungs (Fig. 1a). Substantial lesions appeared in the pointed lobe of the diseased lung, with bleeding points on the surface of the lung. The trachea contained a lot of foam-like mucus, and hilar lymph nodes were swollen. Histopathological analysis of the lungs and trachea showed proliferation of parabronchial lymphoid follicles or lymphocyte infiltration in diseased lungs. The bronchi were filled with pus, mixed with monocytes, neutrophils, and lymphocytes. The alveoli collapsed to varying degrees, and small abscesses were often seen in the collapsed tissue area. Congestion and bleeding of small and medium-sized arteries in the pulmonary interstitium, as well as congestion of capillaries were visible (Fig. 1b). A 5-point scoring system was established on the basis of the severity of each lesion in lung tissue slices; Fig. 1c presents a standard example. These scores data were collected as phenotypic traits for subsequent research.

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)

Based on the epidemiological, clinical anatomical, and histopathological results mentioned above, we grouped all 912 observed sheep into healthy lung (hereinafter, healthy group) and diseased lung (hereinafter, diseased group) groups, and compiled a pneumonia phenotype statistical table (Table 1). It can be seen that the average incidence of pulmonary lesions in sheep caused by respiratory diseases was about 35%. In addition, from July 2022 to January 2023, pneumonia in sheep farms was gradually spreading.

According to the phenotype grouping, we statistically analyzed the growth traits (Table 2) and blood physiological indicators (Table 3) of 912 sheep for differences between groups. We measured each individual’s BW, FI, live weight at slaughter, and carcass weight, and calculated the corresponding ADFI, BWG, ADG, and FCR. No between-group BW differences were noted in the 80-day-old Hu sheep; however, in the 180-day-old Hu sheep, the diseased group demonstrated a extremely significant decrease in BW compared with the healthy group (p < 0.01); the live weight before slaughter and carcass weight were also significantly lower than those in the healthy group (both p < 0.01). Moreover, throughout the feeding trial period, the ADFI and ADG of the diseased Hu sheep remained significantly lower than those of the healthy group (all p < 0.01). In contrast, the FCR in the diseased group was higher than in the healthy group; however, this difference was nonsignificant. We also measured 11 blood physiological indicators related to inflammation. Although most indicators such as lymphocyte (LYMPH), monocyte (MONO), eosinophil (EO), basophil (BASO), EO% and BASO% did not show statistically significant differences between groups, the total white blood cell count, neutrophil count, and proportion in the diseased group were significantly higher than those in the healthy group (all p < 0.01).

Table 2 Statistical analysis of differences in growth data between healthy and diseased lung groups of sheep

Table 3 Statistical analysis of differences in blood indicators between healthy and diseased lung groups of sheep

Pathogen identification method based on high-throughput sequencing analysis

High-throughput transcriptome sequencing of two batches of sheep lung tissue yielded average sequencing reads of approximately 82.8 M and 106.9 M, respectively (Additional file 3). Sheep genome-related sequences accounted for 74.7% and 71.3%, respectively (Additional file 4). After excluding the host genome sequence, several major microorganisms that colonize the lungs were identified by comparing exogenous RNA fragments and calculating the RPM of the aligned microorganisms (Additional file 5). The inter group difference analysis of these microorganisms related sequence quantities emphasized the importance of M. ovipneumoniae (p < 0.01 in both batches) and Pasteurella multocida (P. multocida, p < 0.01 in batch 2) in diseased groups of all batches. In addition, except for Mycobacterium abscesses (M. abscessus, p < 0.01 in batch 1) which showed higher sequence quantity in the healthy group, the sequence quantity of other microorganisms was very low and there was no significant difference between the groups (Fig. 2a). Therefore, M. ovipneumoniae and P. multocida were selected as candidate pathogenic microorganisms. The linear regression between sequencing results and histopathological scores showed a good positive correlation between the sequencing quantity of M. ovipneumoniae and histopathological scores (R² = 0.6788, p < 0.01; Fig. 2b); while such correlation was not noted for P. multocida (R² = 0.1556, p = 0.0118; Fig. 2c). Moreover, the sequencing results showed good coverage of the reference genome of M. ovipneumoniae (Fig. 2d); and the genes with high expression levels of Mycoplasma were sequenced. Considering the strong conservation of 16S rRNA, which is not conducive to genetic analysis of closely related strains, this study used LDH to construct a phylogenetic tree. In the phylogenetic tree, the pathogenic strain detected in this study is M. ovipneumoniae, and it has high homology with the isolated strain (NXNK2203) from Ningxia, which is geographically close (Fig. 2e).

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

Meanwhile, after further 16S rRNA sequencing analysis of 20 sheep from the 2nd batch mentioned above, we obtained 1,743,066 clean reads (Additional file 6), and then counted ASVs at various taxonomic levels (Additional file 7) along with a preliminary ASVs abundance tables of 10,336 ASVs in lung samples. The average sequencing depth was 169, which helps to obtain reliable microbial community analysis results. These count data could be annotated to 530 genera. These data were then normalized to relative abundance by using the total-sum scaling method. The evaluation indicators of species richness, diversity, and evolutionary diversity, namely the Chao1, Simpson, and Faith’s PD indexes, as well as the Pielou’s evenness and Good’s coverage indexes, revealed no significant differences in lung tissue between groups (Fig. 3a). PCoA based on Bray–Curtis distances revealed differential separation of microbial communities between the healthy and diseased groups (Fig. 3b). NMDS analysis also demonstrated between-group separation, with a stress value of 0.131 (Fig. 3c). We also evaluated pulmonary microbiota composition at different taxonomic levels. The most predominant phylum in both healthy and diseased lungs was Bacteroidetes (41.97% and 30.75%, respectively), followed by Firmicutes (34.5% and 24.68%, respectively), Proteobacteria (11.43% and 17.56%, respectively), and Tenericutes (1.19% and 16.51%, respectively; Fig. 3d). Moreover, both healthy and diseased lungs demonstrated a predominance of the genus Prevotella (34.36% and 17.73%, respectively); in contrast, compared with diseased lungs, healthy lungs demonstrated significantly lower predominance of Mycoplasma (0.82% vs. 16.26%, p = 0.057) and Pasteurella (0.15% vs. 9.03%, p < 0.05; Fig. 3e). A heatmap displayed the distribution and variability of the top 10 bacterial genera in both groups (Fig. 3f). LEfSe analyze, with an LDA threshold set to 4, revealed that the phyla Prevotellaceae and Bacteroidetes were enriched in the healthy group. The genus Mycoplasma, the family Mycoplasmataceae, and the order Mycoplasmatales were highly abundant in the diseased group (Fig. 3g).

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

The correctness and accuracy of our high-throughput sequencing analysis in identifying M. ovipneumoniae as a pathogen have been further supported by the laboratory isolation of M. ovipneumoniae (strain number GSWWX349), but this result is related to another unpublished study and therefore not shown here.

Verification of M. ovipneumoniae load in experimental sheep lung tissue

To verify the accuracy of high-throughput sequencing analysis in identifying pathogenic bacteria causing pneumonia, we extracted lung DNA from an expanded sheep population (n = 270) and then performed fluorescence quantitative detection of LDH, which characterizes the load of M. ovipneumoniae. Inter group analysis showed that the M. ovipneumoniae load in the diseased group was higher than that in the healthy group (Fig. 4a), and reached an extremely significant level (p < 0.01), which is consistent with the results of identifying pathogens using high-throughput sequencing. The immunohistochemical results of M. ovipneumoniae in lung tissue showed that the negative control group had substantial lesions and a decrease in the number of alveoli, but the entire lung tissue section did not show positive staining for M. ovipneumoniae. After incubating with primary and secondary antibodies, the sheep lung slices showed varying degrees of positive staining (yellow/orange/brown), and positive staining was mainly observed in the cytoplasm of type I and type II alveolar cells in the alveolar septum, indicating the widespread presence of M. ovipneumoniae in the alveolar epithelial cells of diseased lungs (Fig. 4b).

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

GWAS of pulmonary pathological scores

We analyzed the relationship between 266 pathological scores of sheep lung and 23,480,428 filtered high-quality SNP markers (Additional file 8). The Q-Q plot (Fig. 5b) demonstrated that some SNPs departed from the expected probability, possibly because of the corresponding changes in the severity scores of pneumonia lesions. In the Manhattan plot, SNPs with p values below the significance threshold were considered significantly correlated with the phenotypic traits. The GWAS results demonstrated that 59 SNPs (p < 10^− 6) were significantly correlated with pneumonia lesion scores at the genomic level, located on chromosomes 1, 3, 6, 9, 11, 12, 13, 14, 17, 18, 20, and 21 (Fig. 5a). Table 4 presents the annotation information of these 59 significantly correlated SNPs; of them, 10 SNPs were located in the intron or exon regions of the gene, whereas the remaining SNPs were located in noncoding RNA or intergenic regions. Genome-wide linkage disequilibrium (LD) decay analysis of 912 Hu sheep (Additional file 9) showed that r² began to exhibit a trend toward stability when the distance reached 200 kb, suggesting that the genes within 200 kb of the significant SNPs were potential candidates. Finally, except for no genes annotated near the significant SNPs on chromosomes 1 and 18, a total of 30 genes were annotated in the 400-kb region near other significant SNPs in the genome.

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

Table 4 Information list of SNPs associated with the pathological score of sheep lung by GWAS

Transcriptomic analysis for screening candidate key genes associated with M. ovipneumoniae infection

To further investigate the mechanisms underlying resistance to naturally prevalent Mycoplasma pneumonia in our sheep, we performed RNA-Seq on lung tissues from both healthy and diseased group sheep. After removing some samples to appropriately cluster them and perform more accurate DEGs analysis between the groups, we retained 19 and 16 samples from the healthy and diseased groups, respectively. We statistically analyzed the distribution of clean reads aligned to the genome (Additional file 10) and calculated gene expression levels. Based on PCA of gene expression, 35 lung samples formed good clusters according to their respective lung phenotype states (Fig. 6a). In addition, by analyzing the genes expressed in the healthy and diseased groups, we identified 2,727 DEGs in the diseased group, and the volcano plot (Fig. 6b) showed the distribution and expression differences of 1,606 upregulated DEGs and 1,121 downregulated DEGs. Bidirectional clustering analysis achieved clustering of samples with the same phenotype and clustering of DEGs (Fig. 6c). We then subjected the DEGs to GO and KEGG pathway enrichment analyses. GO enrichment analysis results (Fig. 6d) demonstrated that the DEGs were highly enriched in GO terms related to various cellular components. Moreover, they were mainly enriched in GO terms related to molecular functions, such as signaling receptor activity, molecular transducer activity, and transmembrane signaling receptor activity. Finally, they were enriched in GO terms related to biological processes, including stimulus response, cell adhesion, and multicellular organismal process regulation. In the KEGG pathway enrichment analysis (Fig. 6e), we noted that the DEGs were highly enriched in crucial pathways including viral carcinogenesis, epidermal growth factor receptor (EGFR) tyrosine kinase inhibitor resistance, axon guidance, fluid shear stress and atherosclerosis, and nucleocytoplasmic transport. The inter group expression of 10 randomly selected DEGs (PROM1, PAK5, SLIT1, PGAM1, UBTD1, ADGRB3, FOXF1, PRDX6, PTGS1, and SMAD7) showed consistent differential trends with our RNA-Seq results (Fig. 6f), indicating the reliability of our RNA-Seq data.

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

Finally, through overlap analysis, we identified three annotated DEGs (ADGRB3, PAK5, and FOXF1) from GWAS that had 5 significant SNPs. Compared with the healthy group, the expression of PAK5 was increased in diseased lung tissue, while the expression of ADGRB3 and FOXF1 was decreased.

Different genotypes of important SNPs affect FOXF1 gene expression and individual pneumonia phenotype

Five SNPs were genotyped through KASPar, and 2 significant SNPs (chr14: g. 13215150 A > G, chr14: g. 13215157 A > G) were successfully validated (Additional file 11). Both SNPs are annotated to the same differentially expressed gene FOXF1. According to the different genotypes of each important SNP locus, the expanded sheep were divided into three groups, and the expression of the candidate key gene FOXF1 and the sheep lung phenotype (including M. ovipneumoniae load and lung pathological score) were detected in each group.

The association analysis between genotype and M. ovipneumoniae load showed that at the mutation site g. 13,215,150 A > G on chromosome 14, the expression level of FOXF1 gradually decreased in individuals with AA, AG, and GG genotypes, while the M. ovipneumoniae load gradually increased (Fig. 7a). At another mutation site g. 13,215,157 A > G on chromosome 14, FOXF1 expression gradually increased in individuals with AA, AG, and GG genotypes, while the M. ovipneumoniae load gradually decreased (Fig. 7b). According to the genotyping results, sheep individuals exhibit AA genotype at chr.14: 13,215,150 site, while most of them also exhibit GG genotype at chr.14: 13,215,157 site. At this time, FOXF1 expression is the highest compared to other genotype individuals, while the M. ovipneumoniae load is the lowest. This is consistent with the high expression of FOXF1 in healthy lung tissue observed in RNA-Seq analysis. This correlation trend was also validated in the association analysis between genotype and another lung phenotype trait, pathological score (Fig. 7c). What’s more, the analysis of blood physiological indicators for sheep expanding population showed that individuals with GG genotype at the mutation site g. 13,215,150 A > G generally had higher inflammation-related indices (Table 5). In summary, we speculate that FOXF1 may serve as a resistance gene for Mycoplasma pneumonia in sheep.

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

Table 5 Significant SNP and their association analysis with blood parameters

Discussion

Recently, sheep pneumonia incidence has increased in China [36]. Hu sheep is a local breed of sheep raised in China, known for its early maturity and high reproductive efficiency [37]. From August 2022 to January 2023, Hu sheep reared at a study sheep farm developed pneumonia with an increased incidence of inflammatory lesions in the lungs, but no deaths. Such a situation urgently requires a rational assessment of the epidemiologic etiology of pneumonia in this sheep farm.

When describing chronic or acute sheep pneumonia, it is necessary to consider mortality and subclinical effects, as chronic pneumonia may gradually recover, while acute pneumonia has severe symptoms and may even lead to death [14, 38]. The detection of chronic nonprogressive pneumonia usually requires autopsy or postslaughter lung examination. Compared with experimental infection, natural pneumonia occurrence more effectively reflects the disease characteristics under commercial sheep-breeding conditions. Our study strictly starts from the macroscopic of clinical autopsy and histopathological microscopic aspects to classify the phenotype of sheep lungs to distinguish between healthy or inflammatory lesions. We observed pathological changes in sheep lungs, such as lung consolidation, inflammatory cell infiltration, and vascular lesions—similar to the M. ovipneumoniae infection. We confirmed the negative impact of pneumonia on sheep growth and development based on phenotype grouping analysis. This is consistent with reports that some sheep infected with pneumonia have lower ADG, FCR, live weight, and carcass yield [39]. Blood physiological indicators can reflect the health status of organisms from an immunological perspective and are important indicators for disease screening and diagnosis [40]. We found that the white blood cell, neutrophil and its proportion in diseased sheep increased, which is consistent with previous reports related to pneumonia [41, 42]. The cost of pneumonia has prompted the search for effective solutions to combat the disease, and research aims to explore potential resistance genes to cultivate advantageous livestock and poultry breeds with considerable resistance to the disease. Identifying the pathogen of pneumonia is the first step in disease resistance research. But the determination of the pathogen cannot rely solely on pathological examination.

Traditional veterinary pathogen identification first involves the isolation, purification, and cultivation of pathogens, and then relies on biochemical testing, Sanger sequencing, serological testing, and qRT-PCR to identify pathogens [38, 43]. This undoubtedly requires longer cycles and specific culture media, as well as consideration for small and difficult to culture pathogens such as Mycoplasma. It seems promising to achieve high-resolution rapid detection of clinical samples. The rapid development of high-throughput sequencing technology has led to the widespread use of WGS to characterize microbial genomes, especially for clinical pathogen identification and bioinformatics analysis [44, 45]. Respiratory infectious diseases are caused by pathogens invading the respiratory system [46], and although their colonization ability varies by tissue, nasal and lung isolates represent similar biological populations [47]. Direct sequencing of lung tissue not only provides information about the lung microbiome [48], but also obtains a large amount of host genome sequences. It is crucial to use high-throughput sequencing technology to design pathogen identification schemes suitable for different studies when dealing with clinical samples of unknown pathogens. In this study, considering the high transcription levels of active pathogens in diseased tissues, we anticipate that the proportion of microbial sequences in RNA sequencing results will be high. Therefore, we performed high-throughput sequencing of sheep lungs following the lncRNA analysis process. After filtering out the host genome, the obtained exogenous RNA fragments were compared with the microbial genome database, and it was determined that the pathogen of farm epidemic pneumonia was M. ovipneumoniae. The physiological and pathological status of the lungs varies depending on the richness and abundance of the microbiota [49, 50]. The increase in the number of dominant pathogenic bacteria in the lungs leads to the occurrence and deterioration of pneumonia [51]. 16S rRNA sequencing confirmed significant enrichment of Prevotella (34.36%, 17.73%) in the lungs of healthy sheep [52]. The species composition and LEfSe analysis showed significant differences in the abundance and relative contribution of Mycoplasma between the healthy and diseased groups, indicating that it is a biomarker for pneumonia lesions. We wanted to verify this conclusion through experiments. Considering that M. ovipneumoniae is difficult to count using traditional agar plates, a fluorescent quantitative assay protocol for the LDH gene was specific designed to characterize the load of M. ovipneumoniae based on the genome sequencing sequence of the pathogenic strain in this study. The load detection and immunohistochemical localization of M. ovipneumoniae confirmed the accuracy of our high-throughput sequencing analysis in identifying the pathogen. High-throughput sequencing characterizes all RNA or DNA present in the sample, enabling analysis of the entire microbiome, transcriptome, or host genome in diseased animal samples. Andriyan et al. [53]. sequenced the DNA of poultry samples infected with West Nile Virus (WNV) collected in the United States between 2006 and 2007, and identified a new genetic variant strain containing 13 nucleotide deletions. Chen et al. [54]. accidentally discovered a new coronavirus during a survey of domestic poultry using RNA-Seq. Sequencing technology has successfully realized qualitative and quantitative analysis of pathogens, especially diagnosis and treatment of difficult and critical infectious diseases and identification of new unknown pathogens, with its advantages of high flux, high accuracy, etc. It can also be applied to drug resistance gene monitoring, epidemiological tracking and investigation, and has incomparable advantages over traditional pathogenic identification technology [8].

In recent years, about 90% of lambs with chronic bronchopneumonia have been noted to be positive for pulmonary Mycoplasma after slaughter [55]. In 2001, M. ovipneumoniae was detected in 88% of 453 sheep farms in the United States [56]. Over 2007–2019, 16.4% of Mycoplasma strains isolated from small ruminants were identified to be M. ovipneumoniae in France [57]. In 2021, the nasal swabs from > 40% of sheep in Xinjiang, a major sheep-breeding area in China, tested positive for M. ovipneumoniae [58]. M. ovipneumoniae is prevalent worldwide, mainly transmitted through close and repeated respiratory contact, and can persist in the environment and is difficult to eliminate [59]. All our studied sheep shared an identical environment (including fence, water, and feed), which suggests that the differences in pathological changes in sheep lungs are not due to differences in environmental microorganisms, possibly indicating that M. ovipneumoniae infection resistance or recovery abilities can differ among sheep populations fed in identical environments.

The resistance and resilience of animals to infectious diseases is a complex trait with varying genetic patterns. Effective research strategies such as GWAS are needed to analyze the genetic basis of complex traits [60], which can analyze the associations of genetic markers representing phenotypes with variations in a set of DNA samples [61]. The GWAS tool has identified polymorphic loci or genes associated with important economic traits in the sheep genome, such as reproductive quality [62], wool [63], milk [64], and meat yield [65]. However, GWAS-based research related to disease resilience in breeding programs has been insufficient thus far [66]. Bacterial load cannot determine the severity of lesions, and uneven bacterial colonization also makes it difficult to obtain accurate load data. We believe that an intuitive judgment of the severity of lesions can better understand the sheep’s response and tolerance to M. ovipneumoniae infection. The GWAS results of our pathological scoring identified genomic regions and loci associated with Mycoplasma pneumonia lesions, and annotated key genes. To screen for the most promising candidate genes for disease resistance and conduct targeted functional studies, we chose a multi-omics research approach. The integrated analysis of multi-omics such as whole genome/transcriptome has successfully identified some disease resistance genes in the past [67, 68]. Our comparative transcriptome analysis of healthy and diseased lung in study sheep population validated population segregation based on lung phenotypes. The GO enrichment results of the 2727 identified DEGs include response to stimulus, cell adhesion, and signaling receptor activity, which is consistent with the pathogenic mechanisms of Mycoplasma, as adhesion is the first step of Mycoplasma invasion into the host, thereby creating conditions for subsequent colonization and infection by Mycoplasma [69]. In addition, respiratory epithelial cells, as the first barrier for defense mechanisms in animals, have the functions of synthesizing and secreting mucin, degrading/inhibiting inflammatory mediators, regulating local immune responses, exerting physiological or pathological regulatory effects, and signaling across membranes [70]. When Mycoplasma crosses respiratory epithelial cells, it triggers innate immune responses by producing various inflammatory mediators, thereby mediating lung inflammation [71]. This is also consistent with the results showing that DEGs are enriched in the immune system and leukocyte migration. DEGs are enriched in KEGG pathways such as Fc epsilon RI signaling pathway, Axon guidance, Toll − like receptor signaling pathway, etc. In particular, there have been many reports on the inflammatory response triggered by Mycoplasma invasion into host cells through the Toll − like receptor signaling pathway [17, 72]. The enriched pathways identified in our study are also worthy of further investigation. Based on GWAS and transcriptome, we further screened three overlapping DEGs: ADRB3, PAK5, and FOXF1, as well as five significant SNPs of these annotated genes. Two SNPs with mutations were successfully identified through genotyping. Further analysis was conducted on the different genotypes of these two SNPs and their annotated FOXF1 gene expression, as well as their relationship with individual pneumonia phenotypes. We found that the M. ovipneumoniae load was lowest and pulmonary pathological score was lowest in genotypes with high FOXF1 expression, which is fully consistent with the results of transcriptome analysis and further indicates the dominant genotype associated with resistance to M. ovipneumoniae. In summary, our study identified the key gene FOXF1 and 2 important SNPs (chr.14: g. 13215150 A > G, chr.14: g. 13215157 A > G) associated with resistance to M. ovipneumonia lesions.

In recent years, the genetic mechanisms and candidate genes related to resistance traits of Mycoplasma pneumonia in sheep have received increasing attention from scholars both domestically and internationally. Cao et al. [73] have conducted an in-depth analysis of genomic data from 3,938 sheep samples and found that the introgressed alleles in a specific region of PADI2 (chr2: 248,302,667–248,306,614) are associated with anti-pneumonia. Kathryn et al. [23] discovered five SNPs linked to pneumonic lesions through GWAS, with 37 more showing possible significance. Four SNPs were tied to pleurisy, and 11 others were of potential significance. These SNPs are located near genes that play roles in DNA repair and immune defense, some of which have known connections to respiratory diseases. Our findings contribute new SNPs and candidate resistance genes to the study of resistance to Mycoplasma pneumonia in sheep. It is speculated that these newly discovered SNPs may affect sheep’s resistance to Mycoplasma pneumonia and their individual blood immune status by regulating the expression of FOXF1 gene. FOXF1, a transcription factor forkhead box (FOX) family member, is a transcription factor with a role in organ development and acute injury. FOXF1 expression is significantly higher in pulmonary endothelial cells than in endothelial cells in other organs, such as the brain, heart, pancreas, liver, and kidneys [74, 75]; however, lung injury can reduce FOXF1 expression in pulmonary endothelial cells [76]. FOXF1 regulates key genes involved in inflammation, extracellular matrix remodeling, and endothelial barrier function, thereby stimulating lung repair and regeneration [77]. FOXF1-mediated BMP9/ACVRL1 signaling activates lung epithelial progenitor cells, promoting neonatal pulmonary angiogenesis and alveolarization [78]. It has also been proven that FOXF1, as an antifibrotic factor, inhibits pulmonary fibrosis by preventing the conversion of CDH2 to CDH11 cadherin in myofibroblasts [79]. Numerous studies on FOXF1 in the lungs have emphasized its regulatory function in lung diseases, which also supports our conclusion of screening FOXF1 as a candidate gene for resistance to Mycoplasma pneumonia in sheep through multi-omics methods. However, our work also has limitations. The functional validation of the FOXF1 gene in anti-Mycoplasma pneumonia at the cellular level and its resistance mechanism and regulatory pathway against M. ovipneumonia infection are still unknown, which deserves further research.

The research strategy based on multi-omics sequencing for pathogen identification has enriched and rapidly achieved direct and accurate detection of pathogens. Integrating multi-omics data to mine key genes can be used for future disease resistant breeding, but different omics methods and experimental designs make it difficult to integrate and extract meaningful information from the data, and strict requirements are placed on the collection of phenotype data. We believe that the standardization of such research strategies will promote the development of disease resistant breeding in the livestock and poultry industry.

Conclusions

Our study identified M. ovipneumoniae as the primary pathogen responsible for epidemic pneumonia in sheep farms through high-throughput sequencing, pathogen load quantification in tissues, and immunohistochemical assays. We identified genomic mutation sites associated with Mycoplasma pneumonia through GWAS of pulmonary pathological scores. Especially by integrating multi-omics analysis and association analysis of pneumonia phenotype, two important SNPs related to Mycoplasma pneumonia and a candidate gene FOXF1 for anti-Mycoplasma pneumonia were screened and validated.