Jinyan Yang1
Yongjie Tang1
Xueqin Liu1Jinning Zhang1
Muhammad Zahoor Khan2Siyuan Mi1Chuduan Wang1
Ying Yu1*- 1Laboratory of Animal Genetics and Breeding, Ministry of Agriculture and Rural Affairs of China, National Engineering Laboratory of Animal Breeding, College of Animal Science and Technology, China Agricultural University, Beijing, China
- 2Department of Animal Sciences, Faculty of Veterinary and Animal Sciences, University of Agriculture, Dera Ismail Khan, Pakistan
Subclinical bovine mastitis is a pathogenic infection of the breast characterized by a marked decrease in milk production and quality. As it has no obvious clinical symptoms, diagnosis and treatment are challenging. Therefore, searching for biomarkers in cows’ peripheral white blood cells is valuable for preventing and treating subclinical mastitis. Thus, in this study, the transcriptome of peripheral blood from healthy and subclinical mastitis cows was characterized to find the regulatory signatures of bovine subclinical mastitis using RNA-seq. A total of 287 differentially expressed genes (DEGs) and 70 differentially expressed lncRNAs (DELs) were detected, and 37 DELs were documented near known Quantitative Trait Loci (QTL) associated with the mastitis of cows. Bioinformatic analysis indicated that lncRNAs MSTRG25101.2, MSTRG.56327.1, and MSTRG.18968.1, which are adjacent to the SCS QTL and SCC QTL, may be candidate lncRNAs that influence the pathogenesis of mastitis in cows by up-regulating the expression of genes TLR4, NOD2, CXCL8, and OAS2. Moreover, the alternative splicing (AS) pattern of transcriptional sequence differences between healthy cows and subclinical mastitis cows suggested a molecular mechanism of mastitis resistance and susceptibility. A total of 2,212 differential alternative splicing (DAS) events, corresponding to 1,621 unique DAS genes, were identified in both groups and significantly enriched in immune and inflammatory pathways. Of these, 29 DAS genes were subject to regulation by 32 alternative splicing SNPs, showing diverse and specific splicing patterns and events. It is hypothesized that the PIK3C2B and PPRPF8 splice variants associated with AS SNPs (rs42705933 and rs133847062) may be risk factors for susceptibility to bovine subclinical mastitis. Altogether, these key blood markers associated with resistance to subclinical mastitis and SNPs associated with alternative splicing of genes provide the basis for genetic breeding for resistance to subclinical mastitis in cows.
1 Introduction
Bovine mastitis is the foremost production and major economic burden confronted by the global dairy industry (Hoque et al., 2019). The disease may be asymptomatic or symptomatic (subclinical or clinical mastitis), persistent (chronic), or nonpersistent. Infectious agents such as Staphylococcus aureus (S. aureus), Streptococcus agalactiae, and Mycoplasma spp as well as environmental pathogens such as Escherichia coli, Klebsiella pneumoniae, and Enterobacter aerogenes are the main causes of mastitis in cattle (Gao et al., 2017; Hoque et al., 2020, 2019). In countries such as China, 10–40% of mastitis cases are caused by S. aureus (Wang et al., 2018). These bacteria enter the mammary gland and are recognized by the interaction of their pathogen-associated molecular patterns (PAMP), and the resulting induced inflammatory response leads to the shedding of somatic cells (SC) into the milk. It has been demonstrated that both the number of somatic cells count (SCC) and changes in their gene expression are associated with physiological processes in the mammary gland and bacterial infection (Leitner et al., 2000), with 200,000 cells/mL being used as the best cut-off point to distinguish between infected and uninfected (Schepers et al., 1997; Goncalves et al., 2018).
Blood is a mixture of various immune cells, such as lymphocytes, neutrophils, and monocytes, which have an advantage in reflecting traits such as immunity and disease resistance in cows (Chaussabel, 2015). Also, circulating leukocytes in the blood play a key role in the onset, development, and regression of mastitis, as they are the primary source of immune cells attracted to the mammary gland during infection (Cheng et al., 2021). Most importantly, blood is an easily accessible and least damaging tissue sample compared to bovine mammary epithelial cells, which can accurately reflect the physiological status and health of cows during lactation (Bai et al., 2016). Thus, gene expression profiles from blood offer new opportunities to clarify the mechanisms underlying the complex traits of cows. Blood biomarkers play an important role in characterizing the disease state of animals, and many studies have reported altered mRNA and long non-coding RNA abundance in the mammary gland after mastitis, with some stable markers such as lncRNA XIST (Ma et al., 2019a), lncRNA-TUB (Wang H. et al, 2019), LRRC75A-AS1 (Wang et al., 2020b) and lncRNAs PRANCR (Mi et al., 2021) were identified, but the role of gene transcription and complex networks in the blood remains unclear.
The diversity of functional protein products is mostly attributed to gene alternative splicing. Alternative splicing is a fundamental mechanism by which introns in pre-mRNAs are clipped and exons are bonded together in different configurations, leading to alterations in the structure of major transcripts (Baralle and Giudice, 2017). The previous findings showed that well over 60% of human genes undergo substitutive splicing (Modrek et al., 2001) which has a strong association with diseases (Faustino and Cooper, 2003; Urbanski et al., 2018; Yang et al., 2019). In cows, one study found that 4,567 of 21,755 bovine genes are alternatively spliced, and the most common AS event is exon skipping (Chacko and Ranganathan, 2009). These AS events are strongly associated with disease resistance in cows. For example, two differential alternative splicing (DAS) events in the genes EXOC7 and KIF2C affect protein functional domains and are associated with the susceptibility to Mycobacterium avium subspecies paratuberculosis (Li et al., 2022). Furthermore, the specific alternative splicing patterns exhibited in genes SLAMF7 (Ju et al., 2012a), BOLA-DQA2 (Hou et al., 2012), and CD46 (Wang et al., 2016) in mammary gland tissues of cows were associated with mastitis resistance in cows.
Moreover, 15% of point mutations on genomic DNA that lead to genetic diseases impact pre-mRNA splicing (Wang et al., 2016). Several splicing-related single nucleotide polymorphisms (SNP) may directly alter the coding region, leading to aberrant alternative splicing that affects the disease phenotype (Wang et al., 2014; Ren et al., 2021). In cattle, a study found that an independent spontaneous splice site variant in COL2A1 (g.32473300 G > A) was most likely responsible for chondrodysplasia during early fetal development in cows. Similarly, the SNP g.5880C > T, SNP g.18174A > G, and SNP g.10766T > C might affect the binding with splicing-related factors, subsequently causing the production of aberrant splice variant HMGB3-TV1 (Li et al., 2013; Wang et al., 2014), NCF4-TV (Ju et al., 2015), and NCF1-TV1 (Zhang et al., 2015), thus increasing somatic cell scores in cows.
In the current study, we used RNA-seq to characterize gene expression patterns in the peripheral white blood cells transcriptome of healthy cows and subclinical mastitis cows. Furthermore, we obtained several candidate blood-based transcriptional biomarkers based on the comprehensive analysis of gene function annotation and gene expression patterns. Afterward, we further analyzed alternative splicing events at the transcription level and searched for splicing-associated mutations to elucidate the molecular regulatory mechanisms of bovine subclinical mastitis resistance.
2 Materials and methods
2.1 Sample collection and preparation
Three healthy cows and three subclinical mastitis cows were chosen according to three consecutive months of SCC records and other related phenotypes. In this study, cows with milk SCC values below 100,000 cells/mL were regarded as healthy cows. In contrast, cows with SCC values in the range of 200,000 cells/mL and 500,000 cells/mL were considered to be suffering from subclinical mastitis cows. For similarity in the biological background, the cows in this study had similar lactation days and all calved one litter (Supplementary Table S1).
About 20 mL of anticoagulant blood sample was collected from each cow. About 10 mL of anticoagulant blood sample was centrifuged at 3,500 r/min for 15 min and the middle leukocyte layer was collected and stored in 1 mL of TRIzol at -80 °C for RNA extraction. In addition, about 2 mL blood sample was sent to Jinhaikeyu company for HP testing (Beijing, China). A further amount of blood was used for genotyping at GGP Bull 150K BeadChip (Neogene, Lansing, MI, United States).
2.2 Hematological parameters
A Sysmex K-4500 automated hematology analyzer (Sysmex Corporation, Kobe, Japan) was used to test 24 hematological parameters, including White Blood Cells (WBC), Red Blood Cells (RBC), Hemoglobin (HGB), Red blood cell-specific volume (HCT), Mean Corpuscular Volume (MCV), Mean Corpuscular Hemoglobin (MCH), Mean Corpuscular Hemoglobin Concentration (MCHC), Platelet count (PLT), Neutrophil Ratio (NETU%), Neutrophil count (NETU#), Lymphocyte Ratio (LYMPH%), Lymphocyte count (LYMPH#), Monocyte Ratio (MONO%). Monocyte count (MONO #), Eosinophil Ratio (EO%), Eosinophil count (EO#), Basophil Ratio (BASO%), Basophil count (BASO#), Platelet Distribution Width (PDW), Mean Platelet Volume (MPV), Red cell distribution width - stand error (RDW-SD), Red blood cell distribution width - coefficient of variation (RDW- CV), Platelet-Large Cell Ratio (P-LCR), Platelet cubic metric distribution width (PCT).
2.3 RNA extraction and RNA sequencing
Total RNA of peripheral white blood cells for each cow was extracted with the TRIzol reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer’s procedures. The concentration of RNA was estimated using the NanoDrop 2000 (ThermoFisher Scientific, Waltham, MA, United States) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2,100 system was used to assess RNA integrity. The RNA quality was checked for contamination and degradation by a 1% agarose gel. Then, qualified RNA was used to construct RNA-seq libraries. Finally, cDNA libraries were sequenced with 150 bp paired-end reads from Illumina NovaSeq 6,000 platform (Novogene, Beijing, China).
2.4 RNA-Seq data analysis
2.4.1 Reads quality control and mapping
The quality of the raw data was evaluated with FastQC version 0.11.8. The clean data were obtained by using Trimmomatic version 0.38 to filter out reads with adapter sequences and low-quality reads from the raw reads (Bolger et al., 2014). Specifically, remove splice sequences, bases with mass less than three at the 5′ and 3′ ends of reads, and all bases with average mass less than 15 in the window; trim reads with sequence length less than 36. Subsequently, paired-end reads from each sample were aligned to the reference genome using STAR version 2.7.7a (Dobin et al., 2013). The generation of SAM files was sorted into BAM files using SAMtools version 1.9 (Li et al., 2009). Subsequently, each sample was assembled and merged using StringTie version 1.3.5 (Pertea et al., 2015).
2.4.2 Prediction novel lncRNA
In the current study, strict filtering conditions were set on the annotated transcripts to obtain novel lncRNA transcripts. Firstly, transcripts shorter than 200bp and less than two exons were discarded. Then, transcripts with the class codes “i" (intronic lncRNA, ilncRNA), “u" (Intervening lncRNAs, lincRNA), and “x" (antisense lncRNA, lncNAT) were preserved. Finally, protein-coding potential predictions were performed using the Coding-Non-Coding Index (CNCI) (Sun et al., 2013), ORF Length and GC content (LGC) (Wang G. et al., 2019), Coding Potential Calculator (CPC) (Kang et al., 2017), and Predictor of Long non-coding RNAs and messenger RNAs based on an improved K-mer scheme (PLEK) (Li et al., 2014), and only those transcripts at the intersection of the four software tools 132 were selected as novel lncRNA transcripts.
2.4.3 Quantification and Identification of Differentially Expressed Genes
Transcripts were quantified using HTSeq-counts version 2.0.1 and results were expressed as read counts (Anders et al., 2015). Subsequently, normalization of reads and differential expression analysis of reads were performed using the DESeq2 (Love et al., 2014). Differentially expressed genes (DEGs) and differentially expressed lncRNAs (DELs) were screened with the criteria of p-value < 0.05 and log2 (fold change) > 1.5.
2.4.4 Prediction of the cis and trans-target genes of lncRNAs
To explore the functions of lncRNAs, we simultaneously predicted the cis- and trans-target genes of lncRNAs. Protein-coding genes near the transcript positions of lncRNAs (100k upstream and downstream) were selected as cis-target genes of lncRNAs with the BEDTools version2.1.2 (Quinlan and Hall, 2010). Trans-target genes were predicted based on Pearson correlation coefficients between DELs and DEGs calculated by the R package-Hmisc, with screening criteria of r > 0.98 and p-value < 0.05.
2.4.5 Gene Alternative Splicing Analysis
In the alternative splicing analysis, alternative splicing events were detected using rMATS version 4.1.1 (Shen et al., 2014). In this study, five alternative splicing patterns, namely exon skipping (ES), retained intron (RI), mutually exclusive exon (MXE), alternative 5′ splice site (A5SS), and alternative 3’ splice site (A3SS), were detected from healthy cows and subclinical mastitis cows using rMATS software version 4.1.1. Statistical verification (FDR <0.05) was performed when the differences in isoform ratio (IncLevelDifference) of different alternative splice genes (DAS genes) exceeded the defined threshold of 0.01%. Static images of sarcograms were generated outside of IGV using sarcograms, which is a python tool that is part of the MISO package.
2.4.6 Functional annotation of SNPs
In this study, all animals were genotyped using the GGP Bull 150K BeadChip (Neogene, Lansing, MI, United States) and the effects of all SNPs were annotated by ensemble Variant Effect Predictor (VEP) to screen for SNPs associated with alternative splicing. Subsequently, SNPs were localized to target genes based on positional information to identify candidate SNPs associated with alternative splice genes.
2.5 Functional enrichment analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis and gene ontology (GO) analysis of DEGs, target genes of DELs, and DAS genes were performed on the KOBAS 3.0 online website (http://kobas.cbi.pku.edu.cn/kobas3/).
2.6 Construction interaction network of DELs, DEGs, and pathways
A regulatory network of DELs, target genes of DELs, and pathways were constructed using Cytoscape version 3.5.1 (Yin et al., 2020)
2.7 Statistical analysis
Linear regression analysis was performed using GraphPad Prism (version nine; GraphPad Software, San Diego, CA, United States) and Student’s t-test was used to examine significant differences in the expression levels of lncRNAs and mRNAs between healthy and subclinical mastitis cows. The ggpubr was used to analyze correlations between routine blood test parameters and standardized read counts of genes.
3 Results
3.1 Overview of RNA-seq data
RNA-seq was performed on six cDNA libraries using the Illumina NovaSeq 6,000 platform. The total reads mapped to the bovine genome have been summarized in Table 1. A total of 12.79 million raw reads per library and 12.59 million clean-read pairs were generated from all cDNA libraries. The GC contents of the reads ranged from 50 to 52%. The average alignment rate of clean reads to the bovine genome (version: ARS-UCD1.2) was 95.89%, while 90.91% of clean reads were uniquely mapped reads. Reads that uniquely mapped to the reference genome were used for further analysis (Table 1).
TABLE 1. Summary of reads mapped to the bovine genome.
3.2 Identification and characterization of lncRNAs in subclinical mastitis cows vs. healthy cows
The identification of novel lncRNAs transcripts was followed by a series of criteria. Upon transcript filtering and protein-coding potential prediction, 5,510, 5,664, 5,296, and 5,165 lncRNAs transcripts were predicted by CPC, CNCI, LGC, and PLEK, respectively (Figure 1A). Of particular note, 4,424 lncRNAs transcripts were identified as novel lncRNAs transcripts, which were evaluated by four prediction tools. Of the 4,424 identified lncRNAs transcripts, 76.29% had two or three exons, while 0.52% had more than six exons (Figure 1B).
FIGURE 1. Identification and characterization of lncRNAs in cows with subclinical mastitis and healthy cows. (A) The Venn diagram of novel lncRNAs. The intersection part showed the novel lncRNAs predicted by the coding potential prediction software CPC, PLEK, CNCI, and LGC together. (B) Exon number of identified lncRNAs transcripts. (C) Distribution of identified lncRNAs transcripts on all chromosomes. (D) Among 4,421 lncRNA transcripts, the proportion of intronic lncRNA (ilncRNA), Intervening lncRNAs (lincRNA), and antisense lncRNA (lncNAT).
In the peripheral blood of healthy cows and subclinical mastitis cows, a total of 6,365 lncRNAs were screened, including 1,941 known lncRNAs transcripts. These transcripts were widely distributed in all chromosomes of the bovine genome, and the highest number of transcripts was on chromosome five (Figure 1C). lncRNAs were classified according to their genomic location, with intervening lncRNAs having the highest proportion (47.11%) of all novel lncRNAs, followed by antisense lncRNAs (26.92%) and the lowest percentage of intronic lncRNAs (18.85%) (Figure 1D).
3.3 lncRNAs and mRNA expression profile changes in subclinical mastitis cows vs. healthy cows
A total of 287 differentially expressed genes (DEGs) and 70 differentially expressed lncRNAs (DELs) were identified between the two groups (p-value < 0.05, log(fold change) > 1.5) (left panel in Figures 2A,B). Compared with healthy cows, 253 genes and 57 lncRNAs were significantly up-regulated, while 34 genes and 13 lncRNAs were significantly down-regulated in subclinical mastitis cows. The clustering results of the heatmaps showed that the expression patterns of lncRNAs and mRNAs were quite different between healthy and subclinical mastitis cows (right panel in Figures 2A,B).
FIGURE 2. lncRNAs and mRNA expression profile changes in subclinical mastitis cows vs. healthy cows. The number (A) and the heatmap (B) differentially expressed genes (DEGs) and differentially expressed lncRNAs (DELs) in subclinical mastitis (SM) cows vs. healthy (HY) cows (p-value < 0.05, log2 (fold change) > 1.5).
Next, the functions of 287 DEGs were investigated. As shown in Figures 3A, B, 287 DEGs were involved in 301 GO terms (p-value < 0.05), which contained 55 molecular compositions (CC), 190 biological processes (BP), and 56 molecular functions (MF). Furthermore, KEGG results of DEGs showed 17 enriched signaling pathways (p-value < 0.05) (Figure 3B). Among them, the NOD-like receptor signaling pathway was enriched with the highest number of DEGs (Supplementary Table S2). The genes FCAR, PTAFR, and SELPLG are involved in the S. aureus infection pathway, among which SELPLG is also involved in the cell adhesion signaling pathway (Supplementary Table S2). Notably, the genes SLC40A1 and SLC11A2 were involved in the ferroptosis signaling pathway, and SLC40A1 was up-regulated (Figure 3B) in subclinical mastitis cows vs. healthy cows.
FIGURE 3. Results of GO and KEGG analysis of 287 differentially expressed genes (DEGs) in subclinical mastitis (SM) cows vs. healthy (HY) cows. (A) REVIGO summary and visualization of the significant Go terms (p-value < 0.05), including 190 GO terms of biological process (BP), 50 GO terms of cellular components (CC), and 56 GO terms of molecular function (MF). (B) KEGG pathways of DEGs (p-value < 0.05), and expression differences of several DEGs involved in the S. aureus infection, cell adhesion molecules, and ferroptosis pathway.
We also investigated the relationship between DEGs and routine blood test parameters and found that the standardized read counts of genes DCK, ARGE, VNN2, and OAS2 were significantly correlated with White Blood Cell (WBC) count (109/L) (p-value < 0.05) (Figure 4A). Similarly, the standardized read counts of genes GBP6 and TBXAS1 were significantly correlated with the Monocyte Ratio (MONO %) in peripheral blood (p-value < 0.05) (Figure 4B).
FIGURE 4. Linear analysis of normalized read counts for several differentially expressed genes (DEGs). and routine blood test parameters. (A) DEGs are significantly associated with White Blood Cell (WBC) count (109/L). (B) DEGs are significantly associated with Monocyte Ratio (MONO %).
3.4 Combined DELs and DEGs of healthy cows and subclinical mastitis cows for exploring their co-expression
Potential target genes of DELs were predicted by co-localization and co-expression. In this study, 40 cis-target genes and 124 trans-target genes were predicted (Supplementary Table S2). The results of GO analysis showed that target genes of DELs were involved in 197 BP, 44 CC, and 55 MF (Figure 5A). Most of the target genes were located in cytoplasmic vesicles, outside of the plasma membrane and were involved in the positive regulation of interleukin-17 production, inflammatory response, and cell adhesion (Figure 5A).
FIGURE 5. Results of GO and KEGG analysis of lncRNAs cis-and trans-target genes. (A) Significant GO terms (p-value < 0.05) for cis- and trans-target genes of lncRNAs, including 197 GO terms for biological processes (BP), 44 GO terms for cellular components (CC) and 55 GO terms for molecular functions (MF). (B) Significant KEGG pathways for cis and trans target genes of lncRNAs (p-value < 0.05), and expression differences of the differentially expressed genes (DEGs) involved in the enrichment of the NOD-like receptor signaling, S. aureus infection, and cell adhesion pathways. (C) Annotation results of DELs in subclinical mastitis cows compared to healthy cows, with red bars showing DELs associated with reproductive traits and mastitis.
The results of the KEGG analysis showed 12 pathways (p-value < 0.05) (Figure 5B). Most of these pathways were related to immune responses and inflammation, like NOD-like receptor, p53, IL-17, and S. aureus infection pathway. Besides, the genes FCAR, PTAFR, and SELPLG are involved in the S. aureus infection pathway, among which SELPLG is also involved in the cell adhesion signaling pathway (Supplementary Table S2). Notably, most DELs and their target genes showed the same expression trends (Supplementary Table S2).
3.5 Identification of candidate biomarkers associated with subclinical mastitis
In this study, the functions of DELs were further annotated by comparing the genomic positions of DELs and QTL-related regions in Cattle QTLdb. The results showed that 245 DELs were located adjacent to 481 QTLs in Cattle QTLdb (Figure 5C). These QTLs surrounding DELs were associated with milk production performance and reproductive traits, with more than 75% of QTLs being associated with milk (Supplementary Table S3). More importantly, 35 DELs were located adjacent to three QTLs related to mastitis, such as SCC QTL, somatic cell score (SCS) QTL, and white blood cell count QTL (Supplementary Table S3). What’s more, it was found that 32 lncRNAs were adjacent to stillbirth and calving ease QTLs (Supplementary Table S3), of which 15 DELs were also adjacent to white blood cell QTL, SCS QTL, and SCC QTL, such as MSTRG.10101.2, MSTRG.28749.24, MSTRG.46076.1, MSTRG.67501.5, and MSTRG.23805.2, MSTRG.28197.1 (Supplementary Table S3), and so on.
3.6 Co-expression networks revealed the regulatory mechanism of lncRNAs in subclinical mastitis
Previous studies have shown that lncRNAs may impact disease resistance of livestock and poultry by regulating the expression levels of their cis-and trans-target genes. In this study, five signaling pathways were shared by DEGs and DELs, including NOD-like receptor, S. aureus infection, cell adhesion, B cell receptor, and IL-17 signaling pathways, involving 18 DEGs, 37 DELs, and 7 DEGs. As shown in Figure 6, the network included 117 nodes, 133 edges, and lncRNAs MSTRG.25101.1, MSTRG.25098.1, MSTRG.25418.3, MSTRG.56318.1, MSTRG.15193.2, MSTRG.56327.1, MSTRG.67146.1, and MSTRG.54896.1 may play key nodes to regulate the expression levels of core genes OAS1Y, GBP4, OAS2, FCAR, SELPLG, DCK, and ICAM3. This result implied the important role of lncRNAs in bovine subclinical mastitis.
FIGURE 6. Co-expression networks of lncRNAs cis- and trans- target gene pathways; In particular, these five pathways were shared by differentially expressed genes (DEGs) and differentially expressed lncRNAs (DELs), including the NOD-like receptor, S. aureus infection, and cell adhesion (CAMs), B-cell receptor, and IL-17 signaling pathways. Purple rectangles represents DELs, blue diamonds represents DEGs, red ellipses represents critical paths, and different target relationships in the same regulatory network are distinguished by different colored lines.
3.7 Identification and characterization of total alternative splicing events and differential alternative splicing events
A comprehensive analysis of alternative splicing events in peripheral blood was also performed. As shown in Figure 7A, a total of 65,839 alternative splicing (AS) events were found between healthy and subclinical mastitis cows (FDR < 0.05), including 2,212 differential alternative splicing (DAS) events. Five patterns of AS consist of exon skipping (ES), retained intron (RI), mutually exclusive exon (MXE), alternative 5′ splice site (A5SS), and alternative 3’ splice site (A3SS). A total of 2,212 events were segregated into 82.59% ES events, 11.66% MXE events, 2.44% A3SS and A5SS events, and 3.30% RI events (Figure 7A). More than two types of AS events were presented in 139 genes, 89 genes showed both MXE and SE events, 18 genes showed both SE and RI events (Figure 7B), and the genes CCDC51, PIK3CG, MCEMP1, STRADA, and ATP2A3 presented three types of AS events.
FIGURE 7. Identification and characterization of total alternative splicing (AS) events and differential AS (DAS) events. (A) Distribution of five AS types including skipping exon (SE), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), mutually exclusive exon (MXE), and retained intron (RI) among 65,839 AS events and 2,212 unique DAS events (FDR < 0.05). (B) Statistical results of shared and specific AS events in two groups. The bar graph on the left shows the number of raw AS events. The connections between points indicated the presence of intersections between corresponding data sets. (C) The number of five AS event genes in the healthy and subclinical mastitis cows.
Further analysis revealed that 60 and 13 AS genes presented RI events, 16 and 17 AS genes presented A5SS events, 13 and 8 AS genes presented A3SS events, and the same number of genes presented MXE events in healthy and subclinical mastitis cows, respectively (Figure 7C). Notably, the number of SE events increased two-fold in subclinical mastitis cows compared to healthy cows (Figure 7C).
3.8 Several DAS genes were mastitis-specific and involved in immune, inflammatory response, and reproduction-related pathways
Through the integration of DAS events that belong to the same gene, 2,212 DAS events were found to correspond to 1,621 unique genes (Figure 8A). Gene Ontology and KEGG pathway analysis were used to characterize the primary functions of AS genes. In all three functional categories, AS genes were annotated to different GO terms, including 142 GO terms for CC, 455 GO terms for BP, and 172 GO terms for MF (Figure 8B), affecting insulin binding, regulation of bacterial entry into host cells, and response to ionizing radiation. KEGG pathway analysis further revealed that DAS genes were enriched in 154 pathways, which were classified into five categories in the KEGG database.
FIGURE 8. Functional characterization of differential alternative splicing (DAS) genes. (A) Distribution of the five essential AS events detected from 1,621 DAS genes, and a schematic illustration of the five AS models, including skipping exon (SE), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), mutually exclusive exon (MXE), and retained intron (RI). (B) Significant GO terms for 1,621 DAS genes (p-value < 0.05), including 455 GO terms for biological processes (BP), 142 GO terms for cellular components (CC) and 172 GO terms for molecular functions (MF). (C) Significant KEGG pathways of 1,621 DAS genes (p-value < 0.05), divided into five types, with different colors representing different types.
At one hierarchical level of the KEGG results (Figure 8C), a total of 42 signaling pathways associated with the organism system were enriched (Supplementary Table S4). In particular, the highest proportion of signaling pathways related to the immune system were reported, such as chemokines, Toll-like receptor, NOD-like receptor, hematopoietic cell lineage, natural killer cell-mediated cytotoxicity, IL-17, and leukocyte transendothelial migration signaling pathways (Figure 8C), which comprised 134 unique genes, representing 62.62% of the catalog (Supplementary Table S4). As a matter of note, 32 DAS genes were found involved in reproductive trait-related signaling pathways, such as thermogenesis, oxytocin, estrogen, and progesterone-mediated oocyte maturation signaling pathways (Supplementary Table S4). Hierarchical level two showed that 18 signaling pathways are involved in the regulation of mastitis, including MAPK, Ras, cGMP-PKG, NF-κb, phosphatidylinositol, and mTOR, wnt, and TNF pathways (Figure 8C). At the third level, a total of 19 pathways were identified as being related to cellular processes.
Worthy of note, the genes SLC40A1, ACSL6, GCLC, ACSL5, CYBB, VDAC3, and FTH1 were also enriched in the ferroptosis signaling pathway (Figure 8C). Of these, SLC40A1 was up-regulated and also enriched in the ferroptosis signaling pathway as DEGs (Figure 3B).
3.9 SNPs rs42705933 and rs133847062 may impact mastitis resistance by altering the alternative splicing patterns of PIK3C2B and PRPF8
From the results of the annotation of ensemble Variant Effect Predictor (VEP) on all SNPs, 51.97% of the SNPs were intronic variants, 34.00% were intergenic variants, and only 0.02% of the variants were associated with alternative splicing (Supplementary Figure S1), which included spliced donor region variants and spliced polypyrimidine tract variants. Ultimately, a total of 32 SNPs associated with alternative splicing were screened, corresponding to 29 unique alternative splicing genes (Table 2), including 23 genes that exhibited SE events, three genes that exhibited MXE events, and three genes that exhibited A5SS events, A3SS events, and RI events (Table 2).
TABLE 2. AS SNPs of all SNPs and the AS events of their target genes.
To explore the mechanism by which alternative splicing variants regulate subclinical mastitis resistance, it was further investigated the genotype of these 32 AS SNPs between healthy and subclinical mastitis cows. Two AS SNPs were associated with subclinical mastitis resistance, and the genotypes of these AS SNPs were completely dissimilar between the healthy and subclinical mastitis cows (Table 3). Specifically, SNP rs42705933 on chromosome 16 at 2,210,365 bp of gene PIK3C2B had the genotype AA in healthy cows and AG in subclinical mastitis cows. Furthermore, SNP rs133847062 on chromosome 19 at 22,777,382bp of gene PRPF8 had the genotype AC in healthy cows and AA in subclinical mastitis cows (Table 3). Notably, the two DAS genes both presented exon skipping compared to the healthy cows (Figure 9A).
TABLE 3. Genotypes of SNP rs42705933 and rs133847062 in two groups.
FIGURE 9. Interacting networks of alternative splicing genes, differentially expressed genes (DEGs), and the target genes of differentially expressed lncRNA (DELs). (A) differential alternative splicing (DAS) genes PIK3C2B and PRPF8 exhibited skipping exon (SE) events compared to healthy cows. Red represents healthy cows (HY1, HY2, and HY3) and yellow represents subclinical mastitis cows (SM1, SM2, and SM3). Inclevel represents inclusion level for sample replicates, and gene structure graph represents alternative splicing isoforms imputed from GTF files. (B) PPI network of subclinical mastitis resistance-associated genes. Purple nodes represents DEGs, red nodes represents the target genes of lncRNAs that were close to somatic cell count (SCC) and somatic cell score (SCS), and white cell count QTLs, blue nodes represents two DAS genes regulated by SNPs, and gray nodes represents the second shell of interactors. The thickness of the line represents the strength of the interaction, and the arrow represents the up-regulated gene.
To gain more insight into the interactions of these two mastitis resistance-related DAS genes with DEGs, a PPI network was constructed for functional association analysis. As a result, 110 interactions were found to be present in 37 genes (interaction score >0.4). As shown in Figure 9B, some proteins were able to interact, such as GBP4 could interact with OAS1Y and OAS2, MYD88 could interact with S100A9, and TLR4 could interact with CXCL8, NOD2, and PTAFR. On the other hand, one protein interacts with another protein through a third protein. For example, OAS2 could interact with TBXAS1 via PTGS1, PIK3C2B could interact with AREG via EGFR, and CXCR1 could interact with S100A8 via CXCL6. The above results suggest that these proteins play a synergistic role in the regulation of bovine subclinical mastitis.
4 Discussion
In this study, the expression profiles of mRNAs and lncRNAs from healthy cows and subclinical mastitis cows were characterized, and the functions of 287 DEGs, 70 DELs, and 1,621 DAS genes were annotated, obtaining a number of significantly enriched pathways. Notably, the NOD-like receptor signaling pathway, cytokine-cell receptor interaction, Toll-like receptor signaling pathway, and S. aureus infection pathway were also significantly enriched in the available studies of bovine mammary epithelial (BME) cells (Wang et al., 2016; Wang et al., 2020a; Lin et al., 2021), suggesting that these pathways play important regulatory roles in the pathogenesis of mastitis in blood and breast tissue. Subsequently, based on target gene prediction and QTL database annotation, the results showed that the lncRNAs MSTRG25101.2, MSTRG.56327 0.1, and MSTRG.18968.1, which are adjacent to the SCS QTL and SCC QTL, regulate the expression levels of genes NOD2, CXCL8, OAS2, and TLR4, which are involved in the NOD-like receptor signaling pathway, Cytokine-cytokine receptor interaction, and Toll-like receptor signaling pathway.
The innate immune system is the host’s first line of defense against pathogenic microorganisms and relies on pattern recognition receptors (PRRs) (Kawai and Akira, 2010; Celhar et al., 2012; Kawasaki and Kawai, 2014). These PRRs recognize molecular features expressed on microorganisms and the interaction of PRRs with microbe-associated molecular patterns (MAMP) can lead to the expression of pro-inflammatory cytokines and other immunomodulatory molecules (Gilbert et al., 2013). Among several different types of PRRs in cows, Toll-like receptors (TLRs) and NOD-like receptors (NLRs) have been extensively studied (Bhattarai et al., 2018). In the present study, genes NOD2, CXCL8, OAS2, and TLR4 were differentially expressed in the blood transcriptome of subclinical mastitis cows and involved in the NOD-like receptor, Toll-like receptor pathway and were identified as key candidates for this study. TLR4 and NOD2, as key PRRs, and their recognition of MAMP have been shown to be a key event in the development of mammary inflammation (Goldammer et al., 2004; Aitken et al., 2011; Porcherie et al., 2012). CXCL8 is an important neutrophil chelator in BME cells (Wellnitz and Bruckmaier, 2012), induces chemotaxis in the target cells, and has a significant change in expression levels upon TLR4 activation, which could serve as a potential biological for improving the outcome of mastitis markers (Islam et al., 2020).
Ferroptosis is a regulatory form of iron-dependent cell death (Dixon et al., 2012; Zuo et al., 2022), and is characterized by abnormal iron metabolism (Yang and Stockwell, 2016), lethal lipid peroxidation, and reactive oxygen species accumulation. Several studies have demonstrated that ferroptosis plays a vital regulatory role in the pathology of cancers such as lung cancer (Yang and Stockwell, 2016), gastrointestinal cancer (Zhu et al., 2020), and breast cancer (Xu et al., 2021; Zhang et al., 2021). Yet, few studies reported the relationship between ferroptosis and subclinical mastitis of cows, and the precise mechanism of programmed cell death in bovine mastitis is not fully understood. In our study, KEGG results showed that SLC40A1 (Deng et al., 2021), SLC11A2 (Weijiao et al., 2021), ACSL6, GCLC (Chang et al., 2021; Wang et al., 2021), ACSL5, CYBB (Ping et al., 2022), VDAC3 (Yang et al., 2020; Zhu et al., 2021), and FTH1 (Liu et al., 2022) were involved in the ferroptosis signaling pathway (Supplementary Table S4). Among them, SLC40A1 was up-regulated in subclinical mastitis cows compared to healthy cows (Figure 3B). Moreover, SLC40A1, ACSL6, GCLC, ACSL5, CYBB, VDAC3, and FTH1 were also DAS genes. Of these, the SLC40A1 presented SE and RI events and the FTH1 showed four alternative splicing patterns, including SE and RI, MXE, and A3SS. It is hypothesized that alternative splicing may influence the development of subclinical mastitis by altering the alternative spliced pattern of genes associated with ferroptosis.
RNA splicing is an important regulatory mechanism that links trait-associated variants and complex traits. Alternative splicing generates diverse transcripts with significant roles in disease resistance (Ju et al., 2012a; Asselstine et al., 2022) and the metabolic process (Sun et al., 2021) of the cow. In this study, a total of 2,212 differential splicing events were detected, affecting 1,621 different genes, and GO enrichment analysis showed that they mainly affected some genes that are thought to be involved in transcriptional regulation, proteolytic processes, and neurodevelopment (Figure 8C). It was furthermore reported that the two groups exhibited different AS profiles, with a two-fold increase in the number of SE events in subclinical mastitis cows compared to healthy cows (Figure 7C), whereas there was no significant increase in the number of MXE events, A5SS events and A3SS events, or even a decrease in the number of RI events, demonstrating a link between specific alternative splicing patterns and subclinical mastitis in cows, and similar results were obtained in mammary gland tissues infected with S. aureus (Sun et al., 2021). Consequently, we hypothesized that mastitis-specific AS patterns may increase susceptibility to mastitis and that gene-specific AS patterns in the mammary gland of healthy cows may positively influence the disease resistance of cows, which requires further study.
Mutations affected by pre-mRNA splicing account for more than 15% of disease-causing mutations (Caminsky et al., 2014). Previous studies have documented numerous loci of disease susceptibility genes associated with certain specific splicing patterns in cows. For instance, SNPs rs39652267 and rs39631044 in the 3′ flanking region of the gene JAK2 were dramatically associated with SCC and SCS, and SNP rs43046497 on the intron nine of the gene STAT5A was significantly linked to IL-6 (Usman et al., 2014). These results indicated that the alternative splicing process of mastitis resistance-associated genes is mediated by SNPs. In our study, a total of 32 AS SNPs were identified (Table 2), of which two SNPs rs42705933 and rs133847062 had completely different genotypes between healthy cows and subclinical mastitis cows, which corresponded to SE events of PIK3C2B and PRPF8 (Table 3). Notably, neither of the candidate mastitis genes affected by splice variants in this study, PIK3C2B and PRPF8 were simultaneously regulated at the transcriptional level, indicating independent regulation of genes by alternative splicing, while the PPI interaction network showed that these two alternative splice genes interacted with subclinical mastitis-related genes at the protein level (Figure 9B).
Among them, PIK3C2B belongs to class II of the PI3Ks family and performs a pivotal role in the control of membrane trafficking and intracellular signaling (Vanhaesebroeck et al., 2010; Margaria et al., 2019). The splicing factor PRPF8 is critical for breast cell survival and has potential prognostic value in breast cancer (Cao et al., 2022). What is more, these two SNPs are variants of spliced polypyrimidine fragments (Supplementary Table S5), which are recognized by the polypyrimidine fragment binding protein of the spliceosome complex and are essential for the initial recognition of introns in mammals (Riolo et al., 2021), and the expansion of the polypyrimidine fragment enhances the efficiency of splicing (Romfo and Wise, 1997). Therefore, we hypothesized that the SNPs rs42705933 and rs133847062 may act as pathogenic mutations that cause the genes PIK3C2B and PRPF8 to exhibit exon skipping, which in turn affects the resistance of cows to subclinical mastitis, but the detailed regulatory mechanism remains to be verified in subsequent future studies.
5 Conclusion
Our data suggest that genes TLR4, NOD2, CXCL8, and OAS2 are key components involved in the host immune response to bovine mastitis and that the lncRNAs MSTRG25101.2, MSTRG.56327.1, and MSTRG.18968.1 adjacent to the SCC QTL and SCS QTL may potentially regulate expression of these candidate genes. In addition, characterization of peripheral white blood cell transcriptome alternative splicing events indicated that specific types of alternative splicing events appear to be associated with mastitis resistance in cows and that single nucleotide polymorphisms in the DAS genes PIK3C2B and PPRF8 may be risk factors for mastitis susceptibility in cows. The genes and lncRNAs highlighted in this study could serve as expression biomarkers for mastitis that can be used for genetic improvement of dairy cattle for resilience to mastitis. In conclusion, our study provides a basis for further investigation of the molecular mechanisms of resistance and susceptibility to subclinical mastitis in cows.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA839517.
Ethics statement
The animal study was reviewed and approved by the Animal Welfare Committee of China Agricultural University, Beijing, China.
Author contributions
YY and JY contributed to the conception and design of the study. YT, XL, and SM contributed to the data collection. JY, JZ, and XL analyzed and interpreted the data. JY contributed to the writing of the manuscript. YT, MK, and CW made major contributions to the revision and refinement of the article. YY supervised the research and provided funding. All authors contributed to the article and approved the submitted version. All authors contributed to manuscript revision, read, and approved the submitted version.
Funding
This article was financially supported by the National Key R&D Program of China (2021YFD1200900, 2021YFD1200903), NSFC-PSF Joint Project (31961143009), Beijing Dairy Industry Innovation Team (BAIC06), China Agriculture Research System of MOF and MARA, Beijing Natural Science Foundation (6182021), and the Program for Changjiang Scholar and Innovation Research Team in University (IRT-15R62).
Acknowledgments
The authors thank all members of the Laboratory of Animal Molecular and Quantitative Genetics (China Agricultural University) for sample collection work and Gerile Dari, Wenlong Li, Yibing Liu, and others for helpful discussions.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.949850/full#supplementary-material
References
Aitken, S. L., Corl, C. M., and Sordillo, L. M. (2011). Immunopathology of mastitis: Insights into disease recognition and resolution. J. Mammary Gland. Biol. Neoplasia 16 (4), 291–304. doi:10.1007/s10911-011-9230-4
Anders, S., Pyl, P. T., and Huber, W. (2015). HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics 31 (2), 166–169. doi:10.1093/bioinformatics/btu638
Asselstine, V., Medrano, J. F., and Canovas, A. (2022). Identification of novel alternative splicing associated with mastitis disease in Holstein dairy cows using large gap read mapping. BMC Genomics 23 (1), 222. doi:10.1186/s12864-022-08430-x
Bai, X., Zheng, Z., Liu, B., Ji, X., Bai, Y., and Zhang, W. (2016). Whole blood transcriptional profiling comparison between different milk yield of Chinese Holstein cows using RNA-seq data. BMC Genomics 17 (7), 512. doi:10.1186/s12864-016-2901-1
Baralle, F. E., and Giudice, J. (2017). Alternative splicing as a regulator of development and tissue identity. Nat. Rev. Mol. Cell Biol. 18 (7), 437–451. doi:10.1038/nrm.2017.27
Bhattarai, D., Worku, T., Dad, R., Rehman, Z. U., Gong, X., and Zhang, S. (2018). Mechanism of pattern recognition receptors (PRRs) and host pathogen interplay in bovine mastitis. Microb. Pathog. 120, 64–70. doi:10.1016/j.micpath.2018.04.010
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30 (15), 2114–2120. doi:10.1093/bioinformatics/btu170
Caminsky, N., Mucaki, E. J., and Rogan, P. K. (2014). Interpretation of mRNA splicing mutations in genetic disease: Review of the literature and guidelines for information-theoretical analysis. F1000Res. 3, 282. doi:10.12688/f1000research.5654.1
Cao, D., Xue, J., Huang, G., An, J., and An, W. (2022). The role of splicing factor PRPF8 in breast cancer. Technol. Health Care. 30 (S1), 293–301. doi:10.3233/THC-THC228028
Celhar, T., Magalhaes, R., and Fairhurst, A. M. (2012). TLR7 and TLR9 in SLE: When sensing self goes wrong. Immunol. Res. 53 (1-3), 58–77. doi:10.1007/s12026-012-8270-1
Chacko, E., and Ranganathan, S. (2009). Genome-wide analysis of alternative splicing in cow: Implications in bovine as a model for human diseases. BMC Genomics 10, S11. Suppl 3. doi:10.1186/1471-2164-10-S3-S11
Chang, K., Yuan, C., and Liu, X. (2021). Ferroptosis-Related gene signature accurately predicts survival outcomes in patients with Clear-Cell renal cell carcinoma. Front. Oncol. 11, 649347. doi:10.3389/fonc.2021.649347
Chaussabel, D. (2015). Assessment of immune status using blood transcriptomics and potential implications for global health. Semin. Immunol. 27 (1), 58–66. doi:10.1016/j.smim.2015.03.002
Cheng, Z., Buggiotti, L., Salavati, M., Marchitelli, C., Palma-Vera, S., Wylie, A., et al. (2021). Global transcriptomic profiles of circulating leucocytes in early lactation cows with clinical or subclinical mastitis. Mol. Biol. Rep. 48 (5), 4611–4623. doi:10.1007/s11033-021-06494-8
Deng, S., Zheng, Y., Mo, Y., Xu, X., Li, Y., Zhang, Y., et al. (2021). Ferroptosis suppressive genes correlate with immunosuppression in glioblastoma. World Neurosurg. 152, e436–e448. doi:10.1016/j.wneu.2021.05.098
Dixon, S. J., Lemberg, K. M., Lamprecht, M. R., Skouta, R., Zaitsev, E. M., Gleason, C. E., et al. (2012). Ferroptosis: An iron-dependent form of nonapoptotic cell death. Cell 149 (5), 1060–1072. doi:10.1016/j.cell.2012.03.042
Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). Star: Ultrafast universal RNA-seq aligner. Bioinformatics 29 (1), 15–21. doi:10.1093/bioinformatics/bts635
Faustino, N. A., and Cooper, T. A. (2003). Pre-mRNA splicing and human disease. Genes Dev. 17 (4), 419–437. doi:10.1101/gad.1048803
Gao, J., Barkema, H. W., Zhang, L., Liu, G., Deng, Z., Cai, L., et al. (2017). Incidence of clinical mastitis and distribution of pathogens on large Chinese dairy farms. J. Dairy Sci. 100 (6), 4797–4806. doi:10.3168/jds.2016-12334
Gilbert, F. B., Cunha, P., Jensen, K., Glass, E. J., Foucras, G., Robert-Granie, C., et al. (2013). Differential response of bovine mammary epithelial cells to Staphylococcus aureus or Escherichia coli agonists of the innate immune system. Vet. Res. 44, 40. doi:10.1186/1297-9716-44-40
Goldammer, T., Zerbe, H., Molenaar, A., Schuberth, H. J., Brunner, R. M., Kata, S. R., et al. (2004). Mastitis increases mammary mRNA abundance of beta-defensin 5, toll-like-receptor 2 (TLR2), and TLR4 but not TLR9 in cattle. Clin. Diagn. Lab. Immunol. 11 (1), 174–185. doi:10.1128/cdli.11.1.174-185.2004
Goncalves, J. L., Cue, R. I., Botaro, B. G., Horst, J. A., Valloto, A. A., and Santos, M. V. (2018). Milk losses associated with somatic cell counts by parity and stage of lactation. J. Dairy Sci. 101 (5), 4357–4366. doi:10.3168/jds.2017-13286
Hoque, M. N., Istiaq, A., Clement, R. A., Gibson, K. M., Saha, O., Islam, O. K., et al. (2020). Insights into the resistome of bovine clinical mastitis microbiome, a key factor in disease complication. Front. Microbiol. 11, 860. doi:10.3389/fmicb.2020.00860
Hoque, M. N., Istiaq, A., Clement, R. A., Sultana, M., Crandall, K. A., Siddiki, A. Z., et al. (2019). Metagenomic deep sequencing reveals association of microbiome signature with functional biases in bovine mastitis. Sci. Rep. 9 (1), 13536. doi:10.1038/s41598-019-49468-4
Hou, Q., Huang, J., Ju, Z., Li, Q., Li, L., Wang, C., et al. (2012). Identification of splice variants targeted microRNAs and functional single nucleotide polymorphisms of the BOLA-DQA2 gene in dairy cattle. DNA Cell Biol. 31 (5), 739–744. doi:10.1089/dna.2011.1402
Islam, M. A., Takagi, M., Fukuyama, K., Komatsu, R., Albarracin, L., Nochi, T., et al. (2020). Transcriptome analysis of the inflammatory responses of bovine mammary epithelial cells: Exploring immunomodulatory target genes for bovine mastitis. Pathogens 9 (3), E200. doi:10.3390/pathogens9030200
Ju, Z., Wang, C., Li, Q., Hou, M., Gao, S., Hou, Q., et al. (2012a). Alternative splicing and mRNA expression analysis of bovine SLAMF7 gene in healthy and mastitis mammary tissues. Mol. Biol. Rep. 39 (4), 4155–4161. doi:10.1007/s11033-011-1198-z
Ju, Z., Wang, C., Wang, X., Yang, C., Sun, Y., Jiang, Q., et al. (2015). Role of an SNP in alternative splicing of bovine NCF4 and mastitis susceptibility. PLoS One 10 (11), e0143705. doi:10.1371/journal.pone.0143705
Kang, Y. J., Yang, D. C., Kong, L., Hou, M., Meng, Y. Q., Wei, L., et al. (2017). CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 45 (W1), W12–W16. doi:10.1093/nar/gkx428
Kawai, T., and Akira, S. (2010). The role of pattern-recognition receptors in innate immunity: Update on Toll-like receptors. Nat. Immunol. 11 (5), 373–384. doi:10.1038/ni.1863
Kawasaki, T., and Kawai, T. (2014). Toll-like receptor signaling pathways. Front. Immunol. 5, 461. doi:10.3389/fimmu.2014.00461
Leitner, G., Yadlin, B., Glickman, A., Chaffer, M., and Saran, A. (2000). Systemic and local immune response of cows to intramammary infection with Staphylococcus aureus. Res. Vet. Sci. 69 (2), 181–184. doi:10.1053/rvsc.2000.0409
Li, A., Zhang, J., and Zhou, Z. (2014). Plek: A tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinforma. 15, 311. doi:10.1186/1471-2105-15-311
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25 (16), 2078–2079. doi:10.1093/bioinformatics/btp352
Li, H., Huang, J., Zhang, J., Gao, Y., Han, B., and Sun, D. (2022). Identification of alternative splicing events associated with paratuberculosis in dairy cattle using Multi-Tissue RNA sequencing data. Genes (Basel) 13 (3), 497. doi:10.3390/genes13030497
Li, L., Huang, J., Ju, Z., Li, Q., Wang, C., Qi, C., et al. (2013). Multiple promoters and targeted microRNAs direct the expressions of HMGB3 gene transcripts in dairy cattle. Anim. Genet. 44 (3), 241–250. doi:10.1111/age.12007
Lin, C., Zhu, Y., Hao, Z., Xu, H., Li, T., Yang, J., et al. (2021). Genome-Wide analysis of LncRNA in bovine mammary epithelial cell injuries induced by Escherichia coli and Staphylococcus aureus. Int. J. Mol. Sci. 22 (18), 9719. doi:10.3390/ijms22189719
Liu, J., Ren, Z., Yang, L., Zhu, L., Li, Y., Bie, C., et al. (2022). The NSUN5-FTH1/FTL pathway mediates ferroptosis in bone marrow-derived mesenchymal stem cells. Cell Death Discov. 8 (1), 99. doi:10.1038/s41420-022-00902-z
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 (12), 550. doi:10.1186/s13059-014-0550-8
Ma, M., Pei, Y., Wang, X., Feng, J., Zhang, Y., and Gao, M. Q. (2019a). LncRNA XIST mediates bovine mammary epithelial cell inflammatory response via NF-κB/NLRP3 inflammasome pathway. Cell Prolif. 52 (1), e12525. doi:10.1111/cpr.12525
Margaria, J. P., Ratto, E., Gozzelino, L., Li, H., and Hirsch, E. (2019). Class II PI3Ks at the intersection between signal transduction and membrane trafficking. Biomolecules 9 (3), E104. doi:10.3390/biom9030104
Mi, S., Tang, Y., Dari, G., Shi, Y., Zhang, J., Zhang, H., et al. (2021). Transcriptome sequencing analysis for the identification of stable lncRNAs associated with bovine Staphylococcus aureus mastitis. J. Anim. Sci. Biotechnol. 12 (1), 120. doi:10.1186/s40104-021-00639-2
Modrek, B., Resch, A., Grasso, C., and Lee, C. (2001). Genome-wide detection of alternative splicing in expressed sequences of human genes. Nucleic Acids Res. 29 (13), 2850–2859. doi:10.1093/nar/29.13.2850
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33 (3), 290–295. doi:10.1038/nbt.3122
Ping, S., Wang, S., Zhao, Y., He, J., Li, G., Li, D., et al. (2022). Identification and validation of a ferroptosis-related gene signature for predicting survival in skin cutaneous melanoma. Cancer Med. doi:10.1002/cam4.4706
Porcherie, A., Cunha, P., Trotereau, A., Roussel, P., Gilbert, F. B., Rainard, P., et al. (2012). Repertoire of Escherichia coli agonists sensed by innate immunity receptors of the bovine udder and mammary epithelial cells. Vet. Res. 43, 14. doi:10.1186/1297-9716-43-14
Quinlan, A. R., and Hall, I. M. (2010). BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 26 (6), 841–842. doi:10.1093/bioinformatics/btq033
Ren, P., Lu, L., Cai, S., Chen, J., Lin, W., and Han, F. (2021). Alternative splicing: A new cause and potential therapeutic target in autoimmune disease. Front. Immunol. 12, 713540. doi:10.3389/fimmu.2021.713540
Riolo, G., Cantara, S., and Ricci, C. (2021). What's wrong in a jump? Prediction and validation of splice site variants. Methods Protoc. 4 (3), 62. doi:10.3390/mps4030062
Romfo, C. M., and Wise, J. A. (1997). Both the polypyrimidine tract and the 3' splice site function prior to the first step of splicing in fission yeast. Nucleic Acids Res. 25 (22), 4658–4665. doi:10.1093/nar/25.22.4658
Schepers, A. J., Lam, T. J., Schukken, Y. H., Wilmink, J. B., and Hanekamp, W. J. (1997). Estimation of variance components for somatic cell counts to determine thresholds for uninfected quarters. J. Dairy Sci. 80 (8), 1833–1840. doi:10.3168/jds.S0022-0302(97)76118-6
Shen, S., Park, J. W., Lu, Z. X., Lin, L., Henry, M. D., Wu, Y. N., et al. (2014). Rmats: Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc. Natl. Acad. Sci. U. S. A. 111 (51), E5593–E5601. doi:10.1073/pnas.1419161111
Sun, H. Z., Zhu, Z., Zhou, M., Wang, J., Dugan, M., and Guan, L. L. (2021). Gene co-expression and alternative splicing analysis of key metabolic tissues to unravel the regulatory signatures of fatty acid composition in cattle. RNA Biol. 18 (6), 854–862. doi:10.1080/15476286.2020.1824060
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41 (17), e166. doi:10.1093/nar/gkt646
Urbanski, L. M., Leclair, N., and Anczukow, O. (2018). Alternative-splicing defects in cancer: Splicing regulators and their downstream targets, guiding the way to novel cancer therapeutics. Wiley Interdiscip. Rev. RNA 9 (4), e1476. doi:10.1002/wrna.1476
Usman, T., Yu, Y., Liu, C., Wang, X., Zhang, Q., and Wang, Y. (2014). Genetic effects of single nucleotide polymorphisms in JAK2 and STAT5A genes on susceptibility of Chinese Holsteins to mastitis. Mol. Biol. Rep. 41 (12), 8293–8301. doi:10.1007/s11033-014-3730-4
Vanhaesebroeck, B., Guillermet-Guibert, J., Graupera, M., and Bilanges, B. (2010). The emerging mechanisms of isoform-specific PI3K signalling. Nat. Rev. Mol. Cell Biol. 11 (5), 329–341. doi:10.1038/nrm2882
Wang, D., Wei, G., Ma, J., Cheng, S., Jia, L., Song, X., et al. (2021). Identification of the prognostic value of ferroptosis-related gene signature in breast cancer patients. BMC Cancer 21 (1), 645. doi:10.1186/s12885-021-08341-2
Wang, G., Yin, H., Li, B., Yu, C., Wang, F., Xu, X., et al. (2019a). Characterization and identification of long non-coding RNAs based on feature relationship. Bioinformatics 35 (17), 2949–2956. doi:10.1093/bioinformatics/btz008
Wang, H., Wang, X., Li, X., Wang, Q., Qing, S., Zhang, Y., et al. (2019b). A novel long non-coding RNA regulates the immune response in MAC-T cells and contributes to bovine mastitis. FEBS J. 286 (9), 1780–1795. doi:10.1111/febs.14783
Wang, W., Lin, X., Jiang, T., Peng, Z., Xu, J., Yi, L., et al. (2018). Prevalence and characterization of Staphylococcus aureus cultured from raw milk taken from dairy cows with mastitis in beijing, China. Front. Microbiol. 9, 1123. doi:10.3389/fmicb.2018.01123
Wang, X. G., Ju, Z. H., Hou, M. H., Jiang, Q., Yang, C. H., Zhang, Y., et al. (2016a). Deciphering transcriptome and complex alternative splicing transcripts in mammary gland tissues from cows naturally infected with Staphylococcus aureus mastitis. PLoS One 11 (7), e0159719. doi:10.1371/journal.pone.0159719
Wang, X., Su, F., Yu, X., Geng, N., Li, L., Wang, R., et al. (2020a). RNA-Seq whole transcriptome analysis of bovine mammary epithelial cells in response to intracellular Staphylococcus aureus. Front. Vet. Sci. 7, 642. doi:10.3389/fvets.2020.00642
Wang, X., Wang, H., Zhang, R., Li, D., and Gao, M. Q. (2020b). LRRC75A antisense lncRNA1 knockout attenuates inflammatory responses of bovine mammary epithelial cells. Int. J. Biol. Sci. 16 (2), 251–263. doi:10.7150/ijbs.38214
Wang, X., Zhong, J., Gao, Y., Ju, Z., and Huang, J. (2014). A SNP in intron 8 of CD46 causes a novel transcript associated with mastitis in Holsteins. BMC Genomics 15, 630. doi:10.1186/1471-2164-15-630
Weijiao, Y., Fuchun, L., Mengjie, C., Xiaoqing, Q., Hao, L., Yuan, L., et al. (2021). Immune infiltration and a ferroptosis-associated gene signature for predicting the prognosis of patients with endometrial cancer. Aging (Albany NY) 13 (12), 16713–16732. doi:10.18632/aging.203190
Wellnitz, O., and Bruckmaier, R. M. (2012). The innate immune response of the bovine mammary gland to bacterial infection. Vet. J. 192 (2), 148–152. doi:10.1016/j.tvjl.2011.09.013
Xu, Z., Jiang, S., Ma, J., Tang, D., Yan, C., and Fang, K. (2021). Comprehensive analysis of Ferroptosis-Related LncRNAs in breast cancer patients reveals prognostic value and relationship with tumor immune microenvironment. Front. Surg. 8, 742360. doi:10.3389/fsurg.2021.742360
Yang, Q., Zhao, J., Zhang, W., Chen, D., and Wang, Y. (2019). Aberrant alternative splicing in breast cancer. J. Mol. Cell Biol. 11 (10), 920–929. doi:10.1093/jmcb/mjz033
Yang, W. S., and Stockwell, B. R. (2016). Ferroptosis: Death by lipid peroxidation. Trends Cell Biol. 26 (3), 165–176. doi:10.1016/j.tcb.2015.10.014
Yang, Y., Luo, M., Zhang, K., Zhang, J., Gao, T., Connell, D. O., et al. (2020). Nedd4 ubiquitylates VDAC2/3 to suppress erastin-induced ferroptosis in melanoma. Nat. Commun. 11 (1), 433. doi:10.1038/s41467-020-14324-x
Yin, J., Zeng, X., Ai, Z., Yu, M., Wu, Y., and Li, S. (2020). Construction and analysis of a lncRNA-miRNA-mRNA network based on competitive endogenous RNA reveal functional lncRNAs in oral cancer. BMC Med. Genomics 13 (1), 84. doi:10.1186/s12920-020-00741-w
Zhang, K., Ping, L., Du, T., Liang, G., Huang, Y., Li, Z., et al. (2021). A Ferroptosis-Related lncRNAs signature predicts prognosis and immune microenvironment for breast cancer. Front. Mol. Biosci. 8, 678877. doi:10.3389/fmolb.2021.678877
Zhang, Z., Wang, X., Li, R., Ju, Z., Qi, C., Zhang, Y., et al. (2015). Genetic mutations potentially cause two novel NCF1 splice variants up-regulated in the mammary gland, blood and neutrophil of cows infected by Escherichia coli. Microbiol. Res. 174, 24–32. doi:10.1016/j.micres.2015.03.005
Zhu, J. H., De Mello, R. A., Yan, Q. L., Wang, J. W., Chen, Y., Ye, Q. H., et al. (2020). MiR-139-5p/SLC7A11 inhibits the proliferation, invasion and metastasis of pancreatic carcinoma via PI3K/Akt signaling pathway. Biochim. Biophys. Acta. Mol. Basis Dis. 1866 (6), 165747. doi:10.1016/j.bbadis.2020.165747
Zhu, T., Liu, B., Wu, D., Xu, G., and Fan, Y. (2021). Autophagy regulates VDAC3 ubiquitination by FBXW7 to promote Erastin-Induced ferroptosis in acute lymphoblastic leukemia. Front. Cell Dev. Biol. 9, 740884. doi:10.3389/fcell.2021.740884
Keywords: bovine subclinical mastitis, lncRNA, mRNA, alternative splicing, single nucleotide polymorphism, peripheral blood transcriptome
Citation: Yang J, Tang Y, Liu X, Zhang J, Zahoor Khan M, Mi S, Wang C and Yu Y (2022) Characterization of peripheral white blood cells transcriptome to unravel the regulatory signatures of bovine subclinical mastitis resistance. Front. Genet. 13:949850. doi: 10.3389/fgene.2022.949850
Received: 21 May 2022; Accepted: 26 August 2022;
Published: 20 September 2022.
Edited by:
Ann-Kristin Östlund Farrants, Stockholm University, SwedenReviewed by:
Peipei Ma, Shanghai Jiao Tong University, ChinaM. Nazmul Hoque, Bangabandhu Sheikh Mujibur Rahman Agricultural University, Bangladesh
Copyright © 2022 Yang, Tang, Liu, Zhang, Zahoor Khan, Mi, Wang and Yu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ying Yu, yuying@cau.edu.cn