Correction
1 Dec 2016: Wang XG, Ju ZH, Hou MH, Jiang Q, Yang CH, et al. (2016) Correction: Deciphering Transcriptome and Complex Alternative Splicing Transcripts in Mammary Gland Tissues from Cows Naturally Infected with Staphylococcus aureus Mastitis. PLOS ONE 11(12): e0167666. https://doi.org/10.1371/journal.pone.0167666 View correction
Figures
Abstract
Alternative splicing (AS) contributes to the complexity of the mammalian proteome and plays an important role in diseases, including infectious diseases. The differential AS patterns of these transcript sequences between the healthy (HS3A) and mastitic (HS8A) cows naturally infected by Staphylococcus aureus were compared to understand the molecular mechanisms underlying mastitis resistance and susceptibility. In this study, using the Illumina paired-end RNA sequencing method, 1352 differentially expressed genes (DEGs) with higher than twofold changes were found in the HS3A and HS8A mammary gland tissues. Gene ontology and KEGG pathway analyses revealed that the cytokine–cytokine receptor interaction pathway is the most significantly enriched pathway. Approximately 16k annotated unigenes were respectively identified in two libraries, based on the bovine Bos taurus UMD3.1 sequence assembly and search. A total of 52.62% and 51.24% annotated unigenes were alternatively spliced in term of exon skipping, intron retention, alternative 5′ splicing and alternative 3ʹ splicing. Additionally, 1,317 AS unigenes were HS3A-specific, whereas 1,093 AS unigenes were HS8A-specific. Some immune-related genes, such as ITGB6, MYD88, ADA, ACKR1, and TNFRSF1B, and their potential relationships with mastitis were highlighted. From Chromosome 2, 4, 6, 7, 10, 13, 14, 17, and 20, 3.66% (HS3A) and 5.4% (HS8A) novel transcripts, which harbor known quantitative trait locus associated with clinical mastitis, were identified. Many DEGs in the healthy and mastitic mammary glands are involved in immune, defense, and inflammation responses. These DEGs, which exhibit diverse and specific splicing patterns and events, can endow dairy cattle with the potential complex genetic resistance against mastitis.
Citation: Wang XG, Ju ZH, Hou MH, Jiang Q, Yang CH, Zhang Y, et al. (2016) Deciphering Transcriptome and Complex Alternative Splicing Transcripts in Mammary Gland Tissues from Cows Naturally Infected with Staphylococcus aureus Mastitis. PLoS ONE 11(7): e0159719. https://doi.org/10.1371/journal.pone.0159719
Editor: Marinus F.W. te Pas, Wageningen UR Livestock Research, NETHERLANDS
Received: February 22, 2016; Accepted: July 6, 2016; Published: July 26, 2016
Copyright: © 2016 Wang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Data were deposited in the National Center for Biotechnology Information (NCBI) SRA database (SRX1447227) and within the paper and supporting information.
Funding: This study was supported by grants from Youth Talents Training Program of Shandong Academy of Agricultural Sciences (SAAS-YTTP-2014), the National Natural Science Foundation of China (31371255; 31271328; 31401049), and the Cow Innovation Team of the Shandong Province Modern Agricultural Industry Technology System (SDAIT-12-011-02). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Bovine mastitis is an inflammation of the mammary gland invaded and infected by bacteria. This disease results in considerable economic loss and engenders food safety and animal welfare concerns in the dairy industry [1]. The major microorganisms responsible for mastitis are Staphylococcus aureus (Staph. aureus), Streptococcus and Escherichia coli (E. coli). Among these three pathogens, Staph. aureus is the most frequent cause of udder infection [2]. Obtaining insights into the processes of bovine defense and immune response to mastitis could provide new solutions to mastitis infection. Moreover, a genetic strategy based on the molecular mechanism of cow mastitis demonstrates positive effects on the reduction of antibiotic use in dairy cow breeding and improves the safety of milk products [3].
The alternative splicing (AS) of genes is a common phenomenon in mammalian tissues and cell types in response to stimulations in the eukaryon [4–6]. A gene can produce multiple mRNA transcripts and diverse protein isoforms through this process; subsequently, the gene differentiates proteins to display various binding properties, intercellular localizations, and expression regulations, resulting in related, distinct, or even opposing functions [7–9]. Next-generation sequencing technology is a rapid and cost-effective approach to screen functional candidate genes, differentially expressed genes (DEGs), and important signal pathways that preliminarily explain the molecular mechanism in various tissues. This technology is also used to further identify gene AS related to important economic traits of interest [10,11]. Recent genome-wide association studies have reported that approximately 27.22% of genes in the bovine embryo [12], 38.85% of genes in the adipose tissue [9], and >90% of genes in human tissues [13] undergo AS. More importantly, splice variants of many immune-related genes are associated with various diseases, such as bovine mastitis [14–16]. These variants are one of the major determinant factors of diseases [17]. Currently knowledge on the molecular mechanism underlying the individual differences in immune response to bovine mastitis, especially at the later stages of natural infection with pathogens, is still limited. The interaction between mastitis pathogens and the host immune system is extremely complex [18]. We hypothesize that the differences induced by the AS of genes can be used by the immune system of cows to process complex information to initiate host response to invading pathogens. Therefore, distinguishing the transcriptomic characteristics and differential patterns of AS in bovine mammary glands between healthy and mastitic cows naturally infected with Staph. aureus is important.
To investigate the relevant genes involved in bovine mastitis susceptibility and their regulatory patterns, we initially selected mammary glands from healthy and mastitis-infected cow groups to perform transcriptome sequencing using Illumina HiSeqTM 2000 platform. We obtained a number of candidate genes and signal pathways related to inflammation, defense, and immune responses according to the gene functional annotation and comprehensive analysis of patterns of gene expression. Our findings can provide a foundation for further research on the specific functions of candidate genes related to bovine mastitis susceptibility. The results can also elucidate the molecular process and potential mechanism of cow response to natural infection with Staph. aureus.
Materials and Methods
Ethics Statement
All experiments were carried out according to the Regulations for the Administration of Affairs Concerning Experimental Animals published by the Ministry of Science and Technology, China in 2004 and approved by the Animal Care and Use Committee from the Dairy Cattle Research Center, Shandong Academy of Agricultural Sciences, Shandong, P. R. China.
Mammary gland tissues and RNA isolation
Six 3- to 5-year-old Holstein cows were obtained from the standardized dairy farm of Shandong Province, China. Three of these cows were healthy, whereas the other three were mastitic cows infected with Staph. aureus. Briefly, a cow was defined as healthy if the clinical symptoms such as swelling, redness, hardness or pain were not observed in the udder and no main pathogens was examined from the cow’s mammary tissues using culture and PCR methods. The mastitis group used for this study was referred to as those cows with Staph. aureus and milk somatic cell count per mL above 1 million. At the slaughterhouse, a part of fresh mammary gland tissue were used for pathogen identification, another samples were cut, cleaned with RNase-free water, and then immediately frozen in liquid nitrogen until further use. After pathological evaluation, the three healthy and three mastitis-infected mammary glands were pooled as HS3A and HS8A groups, respectively. The total RNA from the two pool samples was extracted using Trizol reagent (Invitrogen) according to the manufacturer’s instructions. The quality of RNA samples was assessed with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The RNA Integrity Number was 7.6–7.9, and the 28S/18S ratio was 1.9, and the OD260/280 ratio was 2.0.
Construction of cDNA library and sequencing
The cDNA libraries of the two groups (HS3A and HS8A) were constructed as follows. First, total RNAs were isolated from the two samples, and the mRNAs were purified and enriched using magnetic beads with oligo (dT). Second, mRNA sequences were fragmented, and first-strand cDNAs were synthesized using the cleaved segments as a template with six base random primers (Illumina). Subsequently, second-strand cDNA synthesis was performed by adding the buffer, dNTPs, RNase H, and DNA polymerase I. Then, the synthesized cDNA was subjected to end-repair and phosphorylation using T4 DNA polymerase, Klenow DNA polymerase and T4 PNK. Third, Poly (A) and sequencing joints were added to the cDNA. Fourth, the products of ligation reaction were purified on a 2% TAE-agarose gel. A range of cDNA fragments (200 ± 25 bp) were selected from the agarose gel. Fifteen rounds of PCR amplification were performed to enrich the cDNA template using PCR Primer PE 1.0 and PE 2.0 (Illumina) with Phusion DNA Polymerase. Finally, the two constructed libraries were sequenced using the PE technology (2 × 75 bp read length) Illumina HiSeqTM 2000 platform (BGI, Shenzhen, China). The randomness of mRNA fragmentation was evaluated with the distribution of reads in the reference genes (Fig A in S1 File).
Data filtering and transcriptome assembly
Raw reads were cleaned by removing adapter sequences as well as reads with too many unknown base calls (N), low complexity, and low-quality bases, and data quality was controlled by a stringent process to improve the accuracy of the transcriptome analysis results. Reads were filtered following three criteria. (1) Reads with adapter contaminant were first removed. (2) Reads in which the percentage of unknown nucleotides was higher than 5%, the corresponding reads were discarded. (3) Reads in which the percentage of bases with a quality score of ≤5 was greater than 50% were eliminated. After clean data from the HS3A and HS8A were generated, transcriptome assemblies were performed using the SOAP2 software [19]. The genome assembly Bos taurus UMD3.1 deposited in Ensembl (http://www.ensembl.org/info/data/ftp/index.html) was used as the reference genome. We gathered the paired-end (PE) reads with one end mapped on the unique contig and the other end located in the gap region to further shorten the remaining gaps. We also performed local assembly with the unmapped end to fill in the small gaps within the scaffolds. Such sequences that contain the least “N” and do not extend on either end were defined as unigenes.
Identification of differentially expressed genes
The expression levels of genes were calculated using the RPKM method to identify the DEGs between the healthy HS3A and mastitic HS8A groups. This method can eliminate the influence of deviations in transcript lengths and sequencing levels [20]. The corresponding formula is RPKM = (106C)/(NL/103), where C is the number of fragments aligned to the exons of the gene, N is the number of total fragments aligned to all the genes, and L is the base number of the gene CDS. In our analysis, the DEGs were screened using a false discovery rate threshold value of ≤0.001 and an absolute value of log2 ratio of ≥1. All fold changes of differentially expressed genes are log2 values.
Enrichment analysis of differentially expressed genes
The potential functions of assembled DEGs were predicted and annotated against the Gene ontology (GO) and KEGG protein databases. The GO terms especially enriched in DEGs were defined using hypergeometric distribution testing and Bonferroni correction with p-value ≤0.05 as the threshold. For the KEGG protein database, detailed information about each gene can be shown in a signal pathway to reveal its molecular regulatory network and metabolic pathway. The enrichment of a pathway is obtained by multiple testing and is considered significant if q-value ≤ 0.05.
Alternative splicing analysis of genes
We used SOAPsplice software with a default setting and RNA-Seq data to detect splice junctions and identify the potential AS patterns of genes (Fig B in S1 File). SOAPsplice uses a novel approach comprising the following steps: (1) identifying as many reasonable splice junction candidates as possible and (2) filtering the false positives with two effective filtering strategies [21]. First, junction sites, which provide information about boundaries and combinations of different exons in a transcript, are detected by SOAPsplice. Then, all junction sites of the same gene are used to distinguish the type of AS event. A brief introduction of the algorithms used to detect the four AS events in this study is provided in Fig C in S1 File.
Prediction of novel transcripts
Transcripts with reads were assembled using Cufflink [22]. Reads that are continuous and overlapping according to the distribution of reads on the reference genome would form a transcriptional activity area (TAR). Different TARs were linked to form an assembled transcript (Fig D in S1 File). If the gene models were found in intergenic regions (200 bp away from upstream or downstream genes), the transcripts were regarded as candidates for novel transcripts.
Results
Illumina paired-end sequencing and transcriptome assembly
RNA was extracted from the healthy and Staph. aureus-infected mammary tissues to generate a comprehensive survey of genes associated with mastitis infection. Each sequencing feature obtained using Illumina PE sequencing technology yielded 2 × 75 bp independent reads from either end of a cDNA fragment. In the present study, a total of 80% and 78.65% sequencing reads in the HS3A and HS8A were mapped onto the reference Bos taurus genome assembly (UMD3.1), respectively. Finally, the read assembly yielded approximate 16k unigenes with an average length of 2k bp in two libraries (Table 1).
Determining subtle transcriptomic changes in the alveolar mammary tissues can provide insights into molecular mechanisms underlying biological processes, such as immune and defense responses against mastitic infection. In this study, 1,352 DEGs were identified, including 602 up-regulated genes and 750 down-regulated genes, and found specific to the HS8A mammary gland tissue (Table A in S2 File). The top 100 DEGs from the aforementioned genes were selected by the reads per kilobase of the exon model per million mapped reads (RPKM) calculation and further ranked based on the expression fold change value [log2(HS8A/HS3A)] of the genes (Table B in S2 File). Out of the 100 top DEGs, 66% (66/100) genes were up-regulated in the HS8A, indicating that these genes play important roles bovine response to mastitis.
A total of 1,352 DEGs were matched and classified into three functional categories, namely, molecular function (MF), biological process (BP), and cellular component (CC), according to the GO classification system. Among these DEGs, 1,028 matched genes were involved in molecular functions and were clustered into 367 classifications. The top 10 significant enrichment GO terms, such as the protein binding, cytokine activity, receptor binding, chemokine receptor activity, and chemokine binding were also significantly enriched (Table A in S3 File). For the BP, the “response to stimulus” term (32.5%, 332 out of 1,023 DEGs) was the most significant enrichment (Table B in S3 File). Moreover, 124, 74, 17, and 55 DEGs were significantly enriched (p < 0.05) and respectively classified into the GO terms immune system process, defense response, regulation of cell migration, and regulation of immune system process. In addition, 38, 30, 28, 22, 69, and 6 DEGs were respectively classified into the GO terms response to wounding, immune response, response to bacterium, inflammatory response, programmed cell death, and innate immune response in spite of their non-significant enrichments (p > 0.05, Table B in S3 File). The 60 DEGs participating in the “immune system process,” “defense response,” “immune response,” and “inflammatory response” GO terms and expressing with above 1.5-fold change in two groups are listed in S4 File. For the CC analysis, “extracellular region” was the most abundant GO term (Table C in S3 File).
To further understand the biological functions of the unigenes, we mapped DEGs to the signal pathways described in KEGG. Consequently, 1,349 genes were assigned to the KEGG database, and these genes were classified into 221 biological pathways based on the reference canonical pathways. Thirty-three highly ranked KEGG pathway were significantly enriched (p<0.05,Table D in S3 File) and seven pathways are associated with immune, such as the “cytokine–cytokine receptor interaction” pathway (q< 0.05) is the most significant pathway with immunological function (Fig 1). These related genes could act as important regulatory knots during mastitis infection.
Genes enclosed in red boxes are up-regulated, whereas those in green boxes are down-regulated in the mammary gland tissues of mastitis-infected with respect to those in the mammary gland tissues of healthy cows.
Novel transcripts and mastitis-related quantitative trait locus (QTLs)
We identified new transcripts because of the imperfection of the current annotation of the available transcriptome database. A sum of 118,045 and 108,279 candidates for novel transcripts from the 29 autosomes, as well as chromosome MT and X, were identified in the healthy and mastitic tissues, respectively. The transcripts of the mastitic group had an average length of 304 bp, 5 reads per bp, and 2 exons. Similarly, the transcripts of the healthy group possessed approximately 306 bp length, 5 reads per bp, and 2 exons. Using the location of the transcripts in the AnimalQTL (thttp://www.animalgenome.org/cgi-bin/QTLdb/BT/search), we found 4,321 (3.66% of total novel transcripts) novel transcripts of the healthy group and 5,848 (5.4%) novel transcripts of the mastitic group from BTA2, 4, 6, 7, 10, 13, 14, 17, and 20, which were located in the QTLs involved in clinical mastitis (Fig 2; S5 File). A total of 154 candidate novel transcripts (TU45009–TU45164) were specific to the mastitic group (Table B in S5 File). These transcripts were located in the BTA2 (45.57–50.89 Mb), harboring the QTL (46.2–53.0 Mbp) associated with clinical mastitis in Norwegian Red cattle [23] and Staph. aureus incidence rate of clinical mastitis in Holstein cows [24]. Another HS8A-specific novel transcript harbored the QTL associated with the duration of mastitis in Polish Holstein–Friesian cows [25]. Specific novel transcripts were identified in the same chromosomes and QTL regions in the two groups. For instance, 677 novel transcripts expressed only in the mastitic group were located in the region BTA14:1.89–16.51 Mb, harboring the QTL that affects somatic cell score and mastitis in Finnish Ayrshire cattle [26]. Moreover, a few transcripts, such as TU23644 and TU25507, were common in the two groups (Tables A and B in S5 File). The results suggest that mastitis-specific novel transcripts are candidate novel genes for mastitis susceptibility in cows.
BTA: Bovine autosome.
Alternative splicing analysis between healthy and mastitic bovine mammary gland tissues
Seven patterns of gene AS include exon skipping, intron retention, alternative 5′ splice site, alternative 3′ splice site, alternative first exon, alternative last exon, and mutually exclusive exon [7], as shown in the Fig A in S1 File. In the present study, four patterns of gene AS, namely, exon skipping, intron retention, alternative 5′ splice site, and alternative 3′ splice site in the mammary gland transcriptome were detected using the SOAP splice program (Fig B in S1 File). A total of 8,549 (52.62%) and 8,325 (51.24%) unigenes were alternatively spliced, exhibiting 34,523 and 30,410 AS events (Fig 3). In total, 1,317 alternatively spliced unigenes were HS3A-specific, whereas 1,093 alternatively spliced unigenes were HS8A-specific. Further classification of these AS events reveal that a total of 3,497 and 3,221 unigenes demonstrated exon skipping in HS3A and HS8A, respectively. Moreover, 4,146 and 3,818 unigenes presented intron retention, 4,979 and 4,662 unigenes displayed alternative 5′ splicing, 5,364 and 5,090 unigenes showed alternative 3′ splicing in HS3A and HS8A, respectively (Fig 3).
The red and green histograms show the numbers of genes and the corresponding numbers of alternative splicing events, respectively.
The DAVID functional clustering analysis showed that 51 co-expressed genes that participate in the immune, defense, and inflammatory responses showed intron retention (Table A in S6 File). Of these genes, 12 genes (MYD88, C2, TFE3, VDAC1, ENPP2, TLR6, OAS2, IL27R, CD163, CD8A, ENPP1, and SAA3) and eight genes (NCF1, ADA, ITGB6, MAP2K3, MLF2, CSF3, POLR3C, and TNFRSF1B) were healthy-group-specific and mastitic-group-specific, respectively (Table 2). The number of intron retention events of each gene ranged from 1 to 14. For example, the alpha-S2-casein Casocidin-1 (ENSBTAG00000005005) gene has the highest number (14) of intron retention events. Alpha-S2-casein Casocidin-1 is an antibacterial peptide that exhibits antibacterial activity and inhibits the growth of E. coli and Staph. aureus [27]. These genes and their variants probably play important roles in mastitis susceptibility or resistance in dairy cattle.
In addition to genes with retained introns, 21 genes involved in the immune, defense, and inflammatory responses were enriched and displayed exon skipping in both the healthy and mastitic groups (Table B in S6 File). Out of these genes, six genes (CD4, CCL28, C1S, CFB, IL1R1, and MIF) were specific to the healthy group and six genes (ADA, CCL5, COLEC12, CYBA, POLM, and THBS1) were specific to the mastitic group, respectively (Table 3). In addition, four genes (BOLA, SWAP70, NFKB2, and TLR4) were co-expressed; and they exhibited the same exon skipping events in the HS3A and HS8A groups. Five genes (CSN1S2, LTF, TMEM173, PSMB8, and CD46) displayed different exon-skipping patterns, generating distinct transcripts (Table B in S6 File). For example, eight exon skipping events of the bovine LTF gene were detected in the healthy group, whereas five exon skipping events were detected in the mastitic group.
In the healthy and mastitic groups, 17 genes (BMPR1A, IL18, SMAD1, CD4, NCF1, TNFRSF1A, CD46, CCL25, TLR4, CARD9, S1PR3, ENPP2, PSEN2, MLF2, LYST, SAMHD1, and SAA3) and 9 genes (PNLIPRP2, NLRP3, IL1R-I, ITIH4, C7, COTL1, HIF1A, IL17D, and C4A) that were involved in immune, defense, and inflammation responses underwent alternative 5′ splicing site pattern (Table C in S6 File). Eight genes shared the same alternative 5′ splicing events, whereas 17 genes were expressed in different alternative 5ʹ splicing events (Table D in S6 File).
Twenty genes displayed healthy-group-specific alternative 3′-splicing events, whereas nineteen genes expressed mastitic-group-specific alternative 3′-splicing site patterns (Table E in S6 File). Between the two groups, 26 genes shared the same AS patterns, whereas 17 genes had different AS patterns (Table F in S6 File).
In total, 109 genes that are associated with immune, inflammation and defense were alternatively spliced. Further, we found eight genes involved in immune, inflammation and defense were differentially expressed (≥1.5-fold change) and alternatively spliced in the healthy and mastitic tissues, such as C7, CCL19, CCL5, CD4, CSN1S2, ACKR1, ITIH4 and PNLIPRP2 genes (Table E in S4 File).
Discussion
In the present study, we performed transcriptome sequencing of mammary glands from healthy cows and mastitis-infected cows, for which we respectively generated comprehensive transcriptome data of approximately 30 million and 29 million clean reads. We analyzed and identified 1,352 DEGs and complex AS patterns of genes between the healthy and mastitic tissues. Many of the genes have various transcripts involved in the inflammation, defense, and immune processes against infection. We also investigated the potential relationships between candidate novel transcripts and known QTLs associated with clinical bovine mastitis.
Using the clustering analyses of GO protein databases for DEGs, we identified a considerable number functional molecular processes related to the immune defense system, such as the immune system process, defense response, immune response, and inflammatory response terms. Among these terms, the immune system process and defense response terms were significantly enriched. We performed KEGG pathway analysis and found that several important pathways, such as cytokine–cytokine receptor interaction (64 DEGs), chemokine signaling pathway (53 DEGs), and complement and coagulation cascade pathways (32 DEGs), were also enriched. The aforementioned pathways are related to immunological functions. Interestingly, the cytokine–cytokine receptor interaction pathway is only a part of the chemokine signaling pathway. Moreover, the results showed that the two enriched signal pathways occurred in the eosinophils, neutrophils, macrophages, and T lymphocytes. In recent years, cytokine–cytokine receptor interaction has been considered a response to infectious agents and has been reported to be involved in the resolution of inflammation, which is linked to clinical mastitis [28]. A total of 39 DEGs participated in the Fc gamma R-mediated phagocytosis pathway, ranking fifth in the list of pathways. Phagocytosis is an essential element of the immune response permitting the elimination of pathogens, cellular debris, and apoptotic cells. This result suggests the importance of the recruitment and activation of macrophages and neutrophils in infectious sites during the later stage of Staph. aureus mastitis.
Many DEGs play important roles in host defense, inflammation, and tissue damage and healing. In a susceptible cow host, the persistence of bacteria pathogens, such as Sta. aureus, results in aberrant and extended inflammation and subsequent destruction of the mammary gland structures. Therefore, the 10 DEGs of the collagen family are important to heal the damaged tissue. For instance, Collagen Type I Alpha 1 (COL1A1) is up-regulated in the mastitic mammary gland and associated with tissue damage and repair during the later stages of infection [29]. Furthermore, the components of the complement system, such as C1S, C2, C3, C4A, C6, C7, and C8 genes, were differentially expressed; they underwent alternative splicing. The complement system is an essential component of the innate immune response, which is activated after molecular patterns associated with microorganisms, abnormal host cells, and modified molecules in the extracellular environment have been recognized. The consequent proteolytic cascade tags the complement activator for elimination and elicits a proinflammatory response, resulting in the recruitment and activation of immune cells from both the innate and adaptive immune systems [30]. These complement genes are good candidates for improving resistance to mastitis.
Several inflammation and immune-related genes were also highly expressed and alternatively spliced in the udder tissue. Duffy (Fy) antigen/receptor for chemokines (ACKR1, also named DARC) exhibited a 1.73-fold change and retained intron, which can bind to most inflammatory chemokines, activate and recruite leukocytes in the immune system [31]. Serin peptidase inhibitors SERPINA1 and SERPINA11 were highly expressed in the mammary glands, exhibiting 2.52-fold and 1.66-fold changes, respectively (Table B in S3 File). SERPINA1 encodes a serum protease inhibitor that can protect tissues against neutrophil attack and is associated with decreased bovine somatic cell score as well as increased milk yield, milk fat yield, and productive life [32]. Milk fat globule-EGF factor 8 (MFGE8) transcript, which has a 2.44-fold change, was highly expressed in HS3A (RPKM = 218) and HS8A (RPKM = 1183). Mfge8, a soluble milk glycoprotein, is known to regulate inflammation and immunity by mediating the clearance of apoptotic lymphocytes and epithelial cells [33]. The expression of Fc fragment of IgG, low affinity IIb, receptor (FCGR2A, also named CD32) in the mastitis-infected bovine tissue was up-regulated by 1.87-fold compared with that in the healthy tissue. CD32 can participate in biological processes, such as immune response, regulation of immune response, and anti-inflammation [34]. Programmed cell death 4 (neoplastic transformation inhibitor, PDCD4) encoding a tumor suppressor, a highly expressed transcript, was up-regulated in HS8A. PDCD4 is a proinflammatory protein that promotes the activation of the transcription factor NF-κB and suppresses interleukin 10 [35]. A 1.53-fold down-regulation of immunoglobulin lambda-like polypeptide 1 (IGLL1) was found in the mastitic tissue. Bovine IGLL1 is homologous to the variable and the constant domain of the immunoglobulin light chain, which is related to immune function, particularly for B cells and γδ T cells [36]. In the mammary gland, ENSBTAG00000031160 (IGLL1) was identified as a novel abundant transcript similar to immunoglobulin lambda-like polypeptide 1 precursor (immunoglobulin-related 14.1 protein), thereby warranting further functional investigation.
Bovine mastitis, which is an inflammation-driven disease of the mammary gland, occurs in response to physical damage or pathogen infection. Staph. aureus infection tends to cause mild, subclinical inflammation and persistent infection. Ineffective pathogen clearance frequently leads to chronic and persistent infection [37]. Clinical mastitis is a complex trait, and the different genes regulating immune responses are known to be pathogen specific [28]. The mounting of inflammation, defense, and activation of the innate immune system are complex processes involving many different genes, regulatory pathways, and networks, as well as a multitude of environmental factors, which all provide defense lines against pathogens. The complexity of this system possibly necessitates very different mechanisms that provide and shape effective, rapid, and reversible responses with limited resources [38]. AS, which is a versatile mechanism of gene expression regulation, is one of the most significant contributors to the functional complexity of the mammalian genome [8,20]. This process may indicate that different pathogens exploit different host mechanisms invasion strategies. In this study, we performed the first analysis of AS complexity in the mammary gland tissues of cows by combining RNA-Seq data and bioinformatics. We also identified 8,549 (52.62%) and 8,325 (51.24%) unigenes were alternatively spliced, exhibiting 34,523 and 30,410 AS events in the healthy and mastitic tissues, respectively. Many of these genes, which are immune-, defense-, and inflammation-related, underwent different and specific splicing patterns and events. For example, 12 healthy-group-specific genes (MYD88, C2, TFE3, VDAC1, ENPP2, TLR6, OAS2, IL27R, CD163, CD8A, ENPP1, and SAA3) and 8 mastitic-group-specific genes (NCF1, ADA, ITGB6, MAP2K3, MLF2, CSF3, POLR3C, and TNFRSF1B) retained introns. Bovine Myeloid Differentiation Primary Response 88 (MYD88), which retained intron 1 in the healthy group, is a critical protein in the lipopolysaccharide (LPS)-induced NF-κB and apoptotic signaling pathways. This protein plays functional roles in transducing LPS signaling from TLR4 to downstream effector molecules involved in NF-κB activation in endothelial cells [39]. Integrin Subunit Beta 6 (ITGB6), which retained intron 15 in the mastitic tissue, was found to be differentially expressed in mastitis-resistant and susceptible sheep infected with Staph. epidermidis [40]. The predicted protein structure of the novel transcript ITGB6 was altered because of intron retention; as a result, the function of ITGB6 was affected and risk of mastitis was increased. TNFRSF1B, which retained intron 8, was assigned to the inflammatory response and immune response GO terms. ADA was not differentially expressed in the two groups. However, ADA retained intron 5 and deleted exon 2 in the mastitic group. The deficiency of ADA causes a severe combined immunodeficiency disease, in which the dysfunctions of both B and T lymphocytes occur with impaired cellular immunity and decreased production of immunoglobulin (http://www.genecards.org/). We speculate that the domain changes in ADA would alter its function via the AS mechanism, thus leading to different susceptibilities to mastitis.
In our recent studies, we have reported that the novel transcripts of genes produced by AS are significant in bovine mastitis resistance and susceptibility [15,16,41,42]. Among these DEGs, genes that participated in the immune, defense, and inflammatory responses were identified by DAVID clustering to show numerous AS events. In this study, several genes previously reported to initiate antimicrobial and antiviral immune responses, such as C-C Motif Chemokine Ligand 5 (CCL5) [43,44], Colec2 [45], Lactotransferrin (LTF) [46–48], CD46 Molecule (CD46) [16], and Neutrophil Cytosolic Factor 1 (NCF1) [42], showed specific AS events in mastitis-infected mammary glands. Among these genes, CCL5 and CSN1S2 were included in the top 100 highly expressed genes in the HS8A and HS3A mammary glands. We established in our previous study that LTF has two splice variants associated with bovine Sta. aureus mastitis [46]. The BOLA-DQA2-SV1 transcript, which belongs to the BOLA class II genes, plays an important role in mastitis resistance among dairy cattle [41]. In the current study, CD46 exhibited four splicing patterns, namely, exon skipping, intro retention, and alternative 5′ and 3′ splicing. We showed in our previous study that CD46 probably plays a critical role in the Strepotococcus mastitis risk among dairy cows via an alternative splicing mechanism caused by a functional mutation in intron 8 [16]. Two splice variants of NCF1 were sharply up-regulated in the mammary tissues, blood, and neutrophils of mastitis-infected cows compared with those of healthy cows. Moreover, a splicing-related SNP was associated with increased milk somatic cell score in cows [42].
Several QTLs have been found to affect resistance to clinical mastitis in dairy cattle populations [23–26]. We also found many novel transcripts, some of which were located in the known QTLs associated with clinical mastitis. For example, the novel 251 bp TU98975 was identified at Chr8:61,371,653–61,371,903 in the mastitic group. This region of chromosome 8, from 61,042,106 bp to 61,507,067 bp, has been found to be related to clinical mastitis in US Holstein cattle population. This chromosome can be used to explain the 10 highest proportions of variance; however, no annotated genes were found in this region [28]. Our finding indicates that the novel transcript TU98975 may be a candidate gene for mastitis resistance, providing a novel explanation to the finding of the previous study.
Conclusion
We infer that gene-specific AS events in healthy bovine mammary glands may have positive effects on immunoregulatory activity by protecting the host from infection. By contrast, mastitic-specific AS events have negative influences on resistance, thereby increasing the susceptibility to mastitis. Therefore, our findings provide a valuable basis for further understanding of the molecular mechanism underlying mastitis resistance and susceptibility in dairy cows. The diverse and specific splicing patterns and events of genes lead to the loss-of-function and gain-of-function of genes and diversity and contribute to the potential complex genetic resistance against mastitis in dairy cows.
Supporting Information
S1 File.
Fig A. Sequencing randomness assessment and distribution statistics of reads mapped onto the reference gene. Fig B. Schematic diagram of seven kinds of alternative splicing. Fig C. Sketches of the algorithms of the four splicing patterns. Fig D. Schematic diagram of the gene structural optimization method.
https://doi.org/10.1371/journal.pone.0159719.s001
(DOCX)
S2 File.
Table A. A total of 1352 DEGs including 602 up-regulated genes and 750 down-regulated genes were identified in the healthy and mastitic groups. Table B. Top 100 highly expressed genes in HS8A and HS3A mammary gland tissues.
https://doi.org/10.1371/journal.pone.0159719.s002
(XLSX)
S3 File.
Table A. DEGs enriched in molecular function category. Table B. DEGs enriched in biological process category. Table C. DEGs enriched in cellular component category.
https://doi.org/10.1371/journal.pone.0159719.s003
(XLSX)
S4 File.
Table A. A total of 124 DEGs participated in immune system process on the basis of GO term. Table B. A total of 74 DEGs participated in defense response on the basis of GO term. Table C. A total of 30 DEGs participated in immune response on the basis of GO term. Table D. A total of 22 DEGs participated in inflammatory response on the basis of GO term. Table E. List of genes participated in immune, inflammation and defense and expressed with above 1.5-fold change in two groups.
https://doi.org/10.1371/journal.pone.0159719.s004
(XLSX)
S5 File.
Table A. Novel transcripts harboring QTLs associated with clinical mastitis in the healthy mammary gland tissue. Table B. Novel transcripts harboring QTLs associated with clinical mastitis in the mastitic mammary gland tissue.
https://doi.org/10.1371/journal.pone.0159719.s005
(XLSX)
S6 File.
Table A. Genes exhibited retained intron patterns in HS3A and HS8A and participated in immune, defense and inflammatory responses. Table B. Overlapped genes with the same and different exon skipping pattern expressed in HS8A and HS3A and participated in immune, defense and inflammatory responses. Table C. Genes with 5’-splicing events were healthy-group-specific and mastitic- group-specific and participated in immune, defense and inflammatory responses. Table D. Alternative 5' splicing site pattern of genes expressed and participated in immune, defense and inflammatory responses in healthy and mastitis groups. Table E. Genes with 3’-splicing events were healthy-group-specific and mastitic- group-specific and participated in immune, defense and inflammatory responses. Table F. Alternative 3' splicing site pattern of genes expressed and participated in immune, defense and inflammatory responses in healthy and mastitis groups.
https://doi.org/10.1371/journal.pone.0159719.s006
(XLSX)
Author Contributions
Conceived and designed the experiments: JMH. Performed the experiments: XGW ZHJ MHH. Analyzed the data: XGW JMH. Contributed reagents/materials/analysis tools: QJ CHY YZ YS RLL CFW JFZ. Wrote the paper: XGW JMH.
References
- 1. Bar D, Tauer LW, Bennett G, Gonzalez RN, Hertl JA, Schukken YH. et al. (2008) The cost of generic clinical mastitis in dairy cows as estimated by using dynamic programming. J Dairy Sci 91:2205–2214. pmid:18487643
- 2. Vasudevan P, Nair MK, Annamalai T, Venkitanarayanan KS (2003) Phenotypic and genotypic characterization of bovine mastitis isolates of Staphylococcus aureus for biofilm formation. Vet Microbiol 92:179–185. pmid:12488081
- 3. Kulkarni AG, Kaliwal BB (2013) Bovine mastitis: a review. International Journal of Recent Scientific Research 4:543–548.
- 4. Xu Q, Modrek B, Lee C (2002) Genome-wide detection of tissue-specific alternative splicing in the human transcriptome. Nucleic Acids Res 30:3754–3766. pmid:12202761
- 5. Merkin J, Russell C, Chen P, Burge CB (2012) Evolutionary dynamics of gene and isoform regulation in Mammalian tissues. Science 338:1593–1599. pmid:23258891
- 6. Mallory MJ, Allon SJ, Qiu J, Gazzara MR, Tapescu I, Martineza NM, et al. (2015) Induced transcription and stability of CELF2 mRNA drives widespread alternative splicing during T-cell signaling. Proc Natl Acad Sci U S A 112:E2139–E2148. pmid:25870297
- 7. Black DL (2003) Mechanisms of alternative pre-messenger RNA splicing. Annu Rev Biochem 72:291–336. pmid:12626338
- 8. Lareau LF, Green RE, Bhatnagar RS, Brenner SE (2004) The evolving roles of alternative splicing. Curr Opin Struct Biol 14:273–282. pmid:15193306
- 9. Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, et al. (2008) Alternative isoform regulation in human tissue transcriptomes. Nature 456:470–476. pmid:18978772
- 10. Mardis ER (2008) Next-generation DNA sequencing methods. Annu Rev Genomics Hum. Genet 9:387–402. pmid:18576944
- 11. Stapley J, Reger J, Feulner PG, Smadja C, Galindo J, Ekblom R, et al. (2010) Adaptation genomics: the next generation. Trends Ecol Evol 25:705–712. pmid:20952088
- 12. Driver AM, Peñagaricano F, Huang W, Ahmad KR, Hackbart KS, Wiltbank MC, et al. (2012) RNA-Seq analysis uncovers transcriptomic variations between morphologically similar in vivo- and in vitro-derived bovine blastocysts. BMC Genomics 13:118. pmid:22452724
- 13. Zhou Y, Sun J, Li C, Wang Y, Li L, Cai H, et al. (2014) Characterization of transcriptional complexity during adipose tissue development in bovines of different ages and sexes. PLoS One 9:e101261. pmid:24983926
- 14. 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:241–250. pmid:23206268
- 15. Wang X, Huang J, Zhao L, Wang C, Ju Z, Li Q, et al. (2012) The exon 29 c. 3535A> T in the alpha-2-macroglobulin gene causing aberrant splice variants is associated with mastitis in dairy cattle. Immunogenetics 64: 807–816. pmid:22923050
- 16. Wang X, Zhong J, Gao Y, Ju Z, Huang J (2014) A SNP in intron 8 of CD46 causes a novel transcript associated with mastitis in Holsteins. BMC Genomics 15:630. pmid:25070150
- 17. Garcia-Blanco MA, Baraniak AP, Lasda EL (2004) Alternative splicing in disease and therapy. Nat Biotechnol 22:535–546. pmid:15122293
- 18. Thompson-Crispi K, Atalla H, Miglior F, Mallard BA (2014) Bovine mastitis: frontiers in immunogenetics. Front Immunol 5:493. pmid:25339959
- 19. Li R, Yu C, Li Y, Lam TW, Yiu SM, Kristiansen K, et al. (2009) SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25:1966–1967. pmid:19497933
- 20. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B (2008) Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature methods 5:621–628. pmid:18516045
- 21. Huang S, Zhang J, Li R, Zhang W, He Z, Lam TW, et al. (2011) SOAPsplice: genome-wide ab initio detection of splice junctions from RNA-Seq data. Front Gene 2:46.
- 22. Roberts A, Pimentel H, Trapnell C, Pachter L (2011) Identification of novel transcripts in annotated genomes using RNA-Seq. Bioinformatics 27:2325–2329. pmid:21697122
- 23. Sodeland M, Kent MP, Olsen HG, Opsal MA, Svendsen M, Sehested E, et al. (2011) Quantitative trait loci for clinical mastitis on chromosomes 2, 6, 14 and 20 in Norwegian Red cattle. Anim Genet 42:457–465. pmid:21906097
- 24. Verbeke J, Van Poucke M, Peelman L, Piepers S, De Vliegher S (2014) Associations between CXCR1 polymorphisms and pathogen-specific incidence rate of clinical mastitis, test-day somatic cell count, and test-day milk yield. J Dairy Sci 97:7927–7939. pmid:25459910
- 25. Wojdak-Maksymiec K, Szyda J, Strabel T (2013) Parity-dependent association between TNF-α and LTF gene polymorphisms and clinical mastitis in dairy cattle. BMC Vet Res 9:114. pmid:23758855
- 26. Schulman NF, Viitala SM, de Koning DJ, Virta J, Maki-Tanila A, Vilkki JH, et al. (2004) Quantitative trait loci for health traits in Finnish Ayrshire cattle. J Dairy Sci 87:443–449. pmid:14762087
- 27. Kusumaningtyas E (2014) The role of milk peptide as antimicrobial agent in supporting health status. Indonesian Bulletin of Animal and Veterinary Sciences 23:63–75.
- 28. Tiezzi F, Parker-Gaddis KL, Cole JB, Clay JS, Maltecca C (2015) A genome-wide association study for clinical mastitis in first parity US Holstein cows using single-step approach and genomic matrix re-weighting procedure. PLoS One 10:e0114919. pmid:25658712
- 29. Huang J, Luo G, Zhang Z, Wang X, Ju Z, Qi C, et al. (2014) iTRAQ-proteomics and bioinformatics analyses of mammary tissue from cows with clinical mastitis due to natural infection with Staphylococci aureus. BMC Genomics 15:839. pmid:25273983
- 30. Bajic G, Degn SE, Thiel S, Andersen GR (2015) Complement activation, regulation, and molecular basis for complement-related diseases. EMBO J 34:2735–57. pmid:26489954
- 31. Horne K, Woolley IJ (2009) Shedding light on DARC: the role of the Duffy antigen/receptor for chemokines in inflammation, infection and malignancy. Inflamm Res 58:431–435. pmid:19290478
- 32. Khatib H, Heifetz E, Dekkers JC (2005) Association of the protease inhibitor gene with production traits in Holstein dairy cattle. J Dairy Sci 88:1208–1213. pmid:15738254
- 33. Atabai K, Jame S, Azhar N, Kuo A, Lam M, McKleroy W, et al. (2009) Mfge8 diminishes the severity of tissue fibrosis in mice by binding and targeting collagen for uptake by macrophages. J Clin Invest 119:3713–3722. pmid:19884654
- 34. Piccioni M, Monari C, Bevilacqua S, Perito S, Bistoni F, Kozel TR, et al. (2011) A critical role for FcgammaRIIB in up-regulation of Fas ligand induced by a microbial polysaccharide. Clin Exp Immunol 165:190–201. pmid:21605112
- 35. Sheedy FJ, Palsson-McDermott E, Hennessy EJ, Martin C, O′Leary JJ, Ruan Q, et al. (2010) Negative regulation of TLR4 via targeting of the proinflammatory tumor suppressor PDCD4 by the microRNA miR-21. Nat Immunol 11:141–147. pmid:19946272
- 36. Cerri RL, Thompson IM, Kim IH, Ealy AD, Hansen PJ, Staples CR, et al. (2012) Effects of lactation and pregnancy on gene expression of endometrium of Holstein cows at day 17 of the estrous cycle or pregnancy. J Dairy Sci 95:5657–5675. pmid:22884349
- 37. Günther J, Esch K, Poschadel N, Petzl W, Zerbe H, Mitterhuemer S, et al. (2011) Comparative kinetics of Escherichia coli- and Staphylococcus aureus-specific activation of key immune pathways in mammary epithelial cells demonstrates that S. aureus elicits a delayed response dominated by interleukin-6 (IL-6) but not by IL-1A or tumor necrosis factor alpha. Infect Immun 79:695–707. pmid:21115717
- 38. Lynn DJ1, Winsor GL, Chan C, Richard N, Laird MR, Barsky A, et al. (2008) InnateDB: facilitating systems-level analyses of the mammalian innate immune response. Mol Syst Biol 4:218. pmid:18766178
- 39. Cates EA, Connor EE, Mosser DM, Bannerman DD (2009) Functional characterization of bovine TIRAP and MyD88 in mediating bacterial lipopolysaccharide-induced endothelial NF-kappaB activation and apoptosis. Comp Immunol Microbiol Infect Dis 32:477–490. pmid:18760477
- 40. Bonnefont CM, Toufeer M, Caubet C, Foulon E, Tasca C, Aurel MR, et al. (2011) Transcriptomic analysis of milk somatic cells in mastitis resistant and susceptible sheep upon challenge with Staphylococcus epidermidis and Staphylococcus aureus. BMC Genomics 12:208. pmid:21527017
- 41. 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:739–744. pmid:22084936
- 42. 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. pmid:25946326
- 43. Tyner JW, Uchida O, Kajiwara N, Kim EY, Patel AC, O'Sullivan MP, et al. (2005) CCL5-CCR5 interaction provides antiapoptotic signals for macrophage survival during viral infection. Nat Med 11:1180–1187. pmid:16208318
- 44. Nesbeth YC, Martinez DG, Toraya S, Scarlett UK, Cubillos-Ruiz JR, Rutkowski MR, et al. (2010) CD4+ T cells elicit host immune responses to MHC class II− ovarian cancer through CCL5 secretion and CD40-mediated licensing of dendritic cells. J Immunol 184:5654–5662. pmid:20400704
- 45. Hösel M, Broxtermann M, Janicki H, Esser K, Arzberger S, Hartmann P, et al. (2012) Toll-like receptor 2-mediated innate immune response in human nonparenchymal liver cells toward adeno-associated viral vectors. Hepatology 55:287–297. pmid:21898480
- 46. Huang JM, Wang ZY, Ju ZH, Wang CF, Li QL, Sun T, et al. (2011) Two splice variants of the bovine lactoferrin gene identified in Staphylococcus aureus isolated from mastitis in dairy cattle. Genet Mol Res 10:3199–3203. pmid:22194176
- 47. Teran LM, Rüggeberg S, Santiago J, Fuentes-Arenas F, Hernández JL, Montes-Vizuet AR, et al. (2012) Immune response to seasonal influenza A virus infection: A proteomic approach. Arch Med Res 43:464–469. pmid:22960859
- 48. Wojdak-Maksymiec K, Szyda J, Strabel T (2013) Parity-dependent association between TNF-α and LTF gene polymorphisms and clinical mastitis in dairy cattle. BMC Vet Res 9:114. pmid:23758855