This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
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.
Bovine respiratory disease (BRD) is the most common and costly infectious disease in the beef and dairy cattle industry. It causes 70–80% morbidity and 40–50% mortality in feedlot cattle in the United States (
Functional genomic methods such as RNA-sequencing-based transcriptomics can provide a global gene expression profile, and their use in BRD studies can accelerate the understanding of disease mechanisms and the development of diagnosis (
Systems biology is one of the suitable methods to better understand the mechanism of diseases (
In the present study, we used previously published RNA-seq data (
RNA sequencing data from feedlot cattle with and without BRD were obtained from the Gene Expression Omnibus (GEO) database at the National Center for Biotechnology Information (NCBI) under the accession number of GSE162156. Moreover, clinical traits of BRD were obtained from the supplementary material section of the original paper (
Quality control of the raw sequencing data was performed using FastQC1 (version 0.11.9). Raw reads were then trimmed by removing low-quality bases and adaptor sequences using the Trimmomatic software (version 0.39) (
Raw gene expression matrix obtained from the previous steps was normalized to log-counts per million (log-cpm) using the “voom” function of the limma package (version 3.46.0) (
To identify significant highly correlated modules with clinical traits of BRD, all 43 samples (18 healthy and 25 BRD) were used for module–trait relationships (MTRs) analysis. Because the gene coexpression analysis is very sensitive to outliers, the distance-based adjacency metrics of samples was calculated and samples with a standardized connectivity < −2.5 were removed, considered as an outlier. In addition, samples and genes with >50% missing entries and genes with zero variance were identified and excluded from the WGCNA analysis. In this study, a signed weighted coexpression network was constructed in which correlation values between 0 and 1 are considered and values <0.50 are considered as negative correlation, and values >0.50 are considered as positive correlation (
In this method, based on the assumption that BRD may cause a topological change in the coexpression patterns of the healthy samples, and that nonpreserved modules between the healthy and disease samples may be biologically related to BRD, the healthy samples (
To determine which modules are biologically related to BRD, all genes in each module were analyzed for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using the Enrichr web tool (
Highly connected genes (hub genes) in a coexpression module are suitable candidates for explaining behavior and biological function of that module. In other words, highly connected intramodular hub genes have the highest degree of connection in a module and are central to modules in a network, and compared with other genes, they have more biological relevance to the module functions (
For detection of the highly connected and central genes in the PPI network based on the hub genes (hub–hub genes), cytoHubba application (version 0.1) was used (
A summary for the RNA-seq data analysis pipeline and the steps for constructing the weighted gene coexpression network is presented in
Schematic pipeline for RNA-seq data analysis and weighted gene coexpression network construction in this study.
To prevent the negative effects of outlier samples on gene coexpression network analysis, after identifying the outliers, two samples (GSM4943645 and GSM4943656) with a standardized connectivity score < −2.5 were removed (
Sample clustering to detect outliers and network topology analysis.
Module–trait relationships analysis.
In order to understand the biological performance of the significant highly correlated modules with clinical traits of BRD, functional enrichment analysis was performed and a total of 356 biological process and 129 KEGG pathways were significantly enriched in the respective modules. The turquoise module had the highest number of enriched terms and pathways, including 305 biological processes and 116 KEGG pathways. The most significant GO term and KEGG pathway in the turquoise module were “neutrophil-mediated immunity” (GO:0002446, adjusted
Functional enrichment analysis results.
In this study, coexpressed genes in both turquoise and purple modules (as significant highly correlated as well as biologically related modules with BRD) were further evaluated. The MM versus GS plots of these modules are presented in
Scatterplots of module membership (MM) versus gene significance (GS) plots.
Protein–protein interaction (PPI) network based on the hub genes of the turquoise module (identified by module–trait relationships method). Larger nodes and orange octagons represent hub–hub genes and transcription factors, respectively.
List of the hub–hub genes in the turquoise and purple modules as significant highly correlated and biologically related modules with bovine respiratory disease (BRD) along with their module memberships (MM) and gene significance (GS) for rectal temperature (identified by module–trait relationships analysis).
Module | |||||
---|---|---|---|---|---|
Turquoise | Purple | ||||
Genes | MM | GS | Genes | MM | GS |
GAPDH | 0.85 | 0.56 | PRF1 | 0.97 | −0.54 |
IL10 | 0.93 | 0.68 | KLRK1 | 0.92 | −0.71 |
STAT3 | 0.89 | 0.65 | IL2RB | 0.95 | −0.62 |
MAPK1 | 0.81 | 0.66 | LCK | 0.85 | −0.55 |
CASP3 | 0.81 | 0.51 | ITK | 0.82 | −0.60 |
LAMP1 | 0.86 | 0.56 | EOMES | 0.82 | −0.64 |
FN1 | 0.72 | 0.45 | KLRD1 | 0.88 | −0.52 |
MAPK14 | 0.91 | 0.64 | CD40LG | 0.81 | −0.52 |
TLR4 | 0.93 | 0.67 | NCR1 | 0.91 | −0.58 |
TLR2 | 0.89 | 0.50 | CCL5 | 0.92 | −0.49 |
CD68 | 0.90 | 0.65 | LOC618565 | 0.74 | −0.36 |
RAC1 | 0.76 | 0.37 | TBX21 | 0.93 | −0.64 |
CD44 | 0.77 | 0.49 | CD8A | 0.88 | −0.62 |
JAK2 | 0.87 | 0.55 | RUNX3 | 0.83 | −0.68 |
IL1B | 0.86 | 0.50 | XCL2 | 0.74 | −0.30 |
CTSD | 0.85 | 0.65 | CCR8 | 0.95 | −0.64 |
MMP9 | 0.93 | 0.62 | CX3CR1 | 0.90 | −0.57 |
GRB2 | 0.82 | 0.42 | SH2D1A | 0.86 | −0.65 |
ANXA5 | 0.81 | 0.46 | GPR55 | 0.82 | −0.76 |
TYROBP | 0.96 | 0.66 | CTSW | 0.92 | −0.57 |
PTPN6 | 0.79 | 0.51 | KIR2DS1 | 0.72 | −0.34 |
RAB5C | 0.90 | 0.57 | NKG7 | 0.93 | −0.54 |
SOD2 | 0.93 | 0.77 | CD96 | 0.95 | −0.61 |
BECN1 | 0.81 | 0.41 | — | — | — |
WAS | 0.94 | 0.68 | — | — | — |
LAMTOR2 | 0.87 | 0.75 | — | — | — |
CYBA | 0.86 | 0.56 | — | — | — |
SOCS3 | 0.94 | 0.66 | — | — | — |
ALDOA | 0.93 | 0.67 | — | — | — |
SPI1 | 0.96 | 0.61 | — | — | — |
RETN | 0.78 | 0.47 | — | — | — |
MYD88 | 0.93 | 0.68 | — | — | — |
VWF | 0.77 | 0.63 | — | — | — |
PKM | 0.97 | 0.66 | — | — | — |
ORM1 | 0.80 | 0.62 | — | — | — |
STOM | 0.91 | 0.52 | — | — | — |
BCL2L1 | 0.81 | 0.70 | — | — | — |
HSP90AA1 | 0.89 | 0.72 | — | — | — |
MMP25 | 0.96 | 0.67 | — | — | — |
LAMP2 | 0.91 | 0.57 | — | — | — |
PSEN1 | 0.82 | 0.47 | — | — | — |
IL18 | 0.81 | 0.58 | — | — | — |
Note that the rectal temperature is one of the most important and widely used clinical signs of BRD.
For module preservation analysis, healthy samples were used as a reference set to construct the weighted gene coexpression network. Two samples, GSM4943661 and GSM4943667 were identified as outliers based on the distance adjacency metrics of samples and then removed (
Module preservation analysis.
To understand the biological significance and the functional difference between highly preserved, semipreserved, and nonpreserved modules, all the 36 modules were subjected to GO terms and KEGG pathway analysis, and a total of 521 biological processes and 158 KEGG pathways were significantly enriched. Enrichment analysis of four highly preserved modules that included magenta, pink, grey60, and midnightblue revealed a total of 108 and 27 biological processes and KEGG pathways, respectively. Among them, the magenta module with 37 biological processes and 17 KEGG pathways was identified as the most significant enriched module. The other two highly preserved modules, including lightgreen and darkred, showed no biological process or KEGG pathway enrichment. Moreover, one biological process was enriched in the yellowgreen module as the only semipreserved module. The complete information of the functional enrichment analysis for the highly preserved and semipreserved modules is provided in
Six nonpreserved modules were identified as potential candidate modules and biologically related to BRD based on the module preservation and functional enrichment analysis, respectively. Then, MM measures were used to identify intramodular hub genes as important and central genes in these modules. A total number of 326, 251, 216, 145, 131, and 99 hub genes were detected in blue, yellow, red, purple, greenyellow, and salmon modules, respectively. The complete list of the hub genes in the nonpreserved modules can be found in
PPI network based on the hub genes of the purple module (identified by module preservation method). Larger nodes and orange octagons represent hub–hub genes and transcription factors, respectively.
List of the hub–hub genes of the nonpreserved and biologically BRD-related modules that were identified by module preservation analysis.
Modules | |||||
---|---|---|---|---|---|
Blue | Greenyellow | Purple | Red | Salmon | Yellow |
LAMP1 | STAG2 | ISG15 | RAB11FIP2 | IL2RB | C5AR1 |
LAMP2 | PIK3CA | STAT1 | TNF | KLRK1 | IL1B |
TLR4 | SMURF2 | IFIH1 | MAU2 | PRF1 | UBA52 |
PSAP | PTEN | RTP4 | UNKL | ITGAL | PTAFR |
ANXA5 | ADAM10 | USP18 | TAF1 | CCL5 | IL15 |
C3 | TSPAN13 | DDX58 | NR2C2 | CCR5 | GRB2 |
PSEN1 | SNX27 | IRF7 | LCK | GZMA | CD68 |
RAP1B | KLHL11 | MX1 | CD2 | NCR1 | PLEK |
CTSB | PIP5K1B | PARP9 | ARIH2 | RUNX3 | SOCS3 |
STOM | SOS2 | DHX58 | ZFYVE20 | S1PR5 | NFKBIA |
ACTR2 | SNX13 | IFI35 | TSC1 | CX3CR1 | MYD88 |
BECN1 | MGAT4A | RSAD2 | ITK | NKG7 | ANXA1 |
CAT | YES1 | IFI44 | SMC5 | IL12RB2 | FCAR |
NPC2 | GSK3B | IFI44L | LCP2 | GZMB | CYTH4 |
CD86 | SGK3 | UBA7 | PHF8 | CCR8 | CRKL |
ATG7 | UBE2R2 | UBE2L6 | JAK3 | PDCD1 | HECW2 |
VPS35 | TMEM30A | EIF2AK2 | RIC1 | TRPM2 | RAF1 |
PRCP | HECTD1 | IRF9 | SETDB1 | TMEM63A | IL18RAP |
RAB1A | TSPAN33 | ISG20 | FMR1 | GZMH | GADD45B |
RAB5A | RNF217 | MX2 | PLCG1 | — | FAS |
GAA | ROCK1 | XAF1 | CHD3 | — | CNR2 |
MGST1 | TAB2 | PARP14 | FBXO4 | — | CCR1 |
CCR2 | PDP1 | PARP12 | FBXL20 | — | HCAR3 |
CD59 | WAPAL | IFIT5 | PIK3CD | — | SNX18 |
CTSC | WAC | STAT2 | FBXO21 | — | SGK1 |
ACTR1A | PI4K2B | TRIM21 | TIA1 | — | IFNAR2 |
TLR7 | SERPINE2 | OAS1X | TNRC6A | — | VNN2 |
LYZ | — | HERC6 | TLR3 | — | DDIT3 |
GM2A | — | CMPK2 | POGZ | — | WDFY3 |
CST3 | — | ZBP1 | PDPR | — | IFNGR2 |
TNFRSF1B | — | DTX3L | CRAMP1L | — | PPP2R5A |
CAPZA1 | — | ZNFX1 | RERE | — | SELL |
RAB6A | — | IFITM3 | RAB11FIP4 | — | KDM4B |
ARF1 | — | RNASEL | CD6 | — | NCF1 |
FTL | — | SAMD9 | ZC3H11A | — | LOC407171 |
APLP2 | — | GBP4 | TNFRSF10D | — | ARRDC4 |
PECAM1 | — | TRIM25 | UBP1 | — | TARM1 |
ATF6 | — | MB21D1 | SLC37A3 | — | VASP |
CSF2RB | — | OAS1Y | IKBKE | — | SLC2A3 |
RAB18 | — | HERC5 | STAT5A | — | BST1 |
SPTLC1 | — | ADAR | NAA16 | — | MCL1 |
KTN1 | — | GBP7 | ZBTB43 | — | NFAM1 |
SHISA5 | — | OAS2 | CAMSAP1 | — | TNIP1 |
MOSPD2 | — | IFI6 | TWF1 | — | GADD45A |
CD163 | — | FOXS1 | ABI2 | — | NUDT3 |
ATP6V1A | — | PML | VAMP5 | — | MEFV |
IFNAR1 | — | CD53 | CDC7 | — | MXD1 |
CAPN2 | — | CHMP5 | UPB1 | — | LONRF3 |
— | — | GPR97 | — | — | KDM6B |
— | — | — | — | — | ICAM3 |
— | — | — | — | — | TCN1 |
BRD is a multifactorial disease that results from the interaction of environmental stressors and infectious agents of BRDC (
Module–trait relationship analysis identified four significant highly correlated modules including turquoise, purple, blue, and brown with clinical traits of BRD. Functional enrichment analysis indicated that three of the four (75%) identified modules including blue, purple, and turquoise modules had significant enriched terms and pathways. This shows the power and high accuracy of signed networks in separating the modules from each other and identifying more significantly enriched terms or pathways in the modules. The significant enriched terms in the blue module were mostly related to basic cellular activities including “tRNA transport,” “transcription DNA-templated,” “DNA metabolism,” “macromolecule biosynthetic,” and “cell cycle.” On the other hand, the turquoise and purple modules had many biological processes and KEGG pathways closely related to the BRD mechanisms. So, our focus on the turquoise and purple modules (identified by MTRs method) as significant highly correlated and key biologically related modules to BRD.
Coregulated genes in the turquoise module were closely related to the mechanisms of the innate immune system as the first defense line against BRD and this shows the great importance of this module during BRD development. We also found a number of enriched terms related to the adaptive immune system in this module. After infection with various viral or bacterial pathogens, extensive complex interactions begin between the host and the pathogen (
Toll-like receptors (TLRs) are the most important signaling maker PRRs and act as the primary sensors of pathogens (
C-type lectin receptor signaling pathway is another pattern recognition receptor related pathway that is activated by CLRs membrane receptors by identifying different carbohydrates such as mannose, glucan, and fucose in viruses, bacteria, and fungi and activates MAP kinases, the transcription factor
The NF-kappa B signaling pathway is a key pathway that acts as a major mediator in inflammatory responses. Activation of the
Furthermore, other important immune-related terms identified in the turquoise module include “JAK-STAT signaling pathway,” “TNF signaling pathway,” “PI3K-Akt signaling pathway,” “regulation of interferon-gamma production,” and “positive regulation of interleukin-8 secretion.” The JAK-STAT signaling pathway is one of the most important intracellular signaling pathways that regulates communication between cytokine transmembrane receptors and the nucleus and is involved in many biological processes in the body, such as immune regulation, cell differentiation/proliferation, apoptosis, and keep homeostasis in inflammatory conditions (
The PI3K–Akt signaling pathway is a key pathway in all mammalian cells that is involved in several processes, such as cell growth, migration, proliferation, and metabolism as well as the role of this pathway in regulatory T-cell development and memory CD8 T-cell differentiation, has been reported (
Interleukin-8 is a chemokine, that is expressed in various cells especially in macrophages.
In response to signaling proinflammatory cytokines and chemokine such as
We also identified several microbicidal mechanisms, including “Phagocytosis,” “Fc gamma R-mediated phagocytosis,” “Phagosome,” “Lysosome,” and “phagosome maturation” in the turquoise module. In the immune system of multicellular organisms, phagocytosis is a cellular process for the elimination of pathogens and dead cell debris, and is an essential for maintaining tissue homeostasis (
The turquoise module also identified some of the major pathways associated with programmed cell death, including “Apoptosis” and “Necroptosis.” Apoptosis is the first type of programmed cell death that clears cells infected with pathogens (especially viruses) to prevent them from replication and spread (
Interestingly, we identified the terms “regulation of nitric oxide biosynthetic process,” “positive regulation of nitric oxide metabolic process,” and “positive regulation of nitric oxide biosynthetic process” in the turquoise module. Nitric oxide (NO) is a natural molecule with antimicrobial properties that is produced by most mammalian cells (
We also identified some pathways in the turquoise module associated with the adaptive immune system such as “B-cell receptor signaling pathway”, “positive regulation of activated T-cell proliferation,” “Th1/Th2 cell differentiation,” and “Th17 cell differentiation.” B cells are vital cells for the humoral response, and research shows that the B-cell receptor signaling pathway is one of the most important pathways in the immune system of animals to respond to infection with viral agents of BRDC (
Intramodular hub genes (especially hub–hub genes) are highly correlated with the biological function of the module. In this regard, hub–hub genes identified in the turquoise module, such as
Moreover, other important hub–hub genes in the turquoise module included
Among the detected TFs,
Several functional terms identified in the purple module included “Natural killer cell-mediated cytotoxicity,” “T-cell receptor signaling pathway,” “cellular defense response,” “Chemokine signaling pathway,” and “T-cell activation.” Interestingly, in addition to the antiviral activity of T cells that have been discussed above, the T-cell receptor signaling pathway has also been observed in challenges with
Module preservation analysis and functional enrichment showed that six modules (blue, greenyellow, purple, red, salmon, and yellow) first changed their connectivity pattern and network density due to BRD and second, were biologically related to the BRD development. As expected, highly preserved and semipreserved modules were more active in basic and general cellular activities such as “Ribosome,” “rRNA processing,” “DNA packaging,” “peptide biosynthetic process,” and “translation.” Therefore, based on highly preserved and semipreserved modules, it is not possible to explain the molecular mechanisms involved in BRD. On the other hand, the six mentioned nonpreserved modules were closely related to the immune system and the pathogenesis of BRD.
Functional terms in the blue module were mostly related to the innate immune system, including “lysosome,” “phagosome,” “Fc gamma R-mediated phagocytosis,” “leukocyte transendothelial migration,” “Toll-like receptor signaling pathway,” “Chemokine signaling pathway,” “neutrophil activation involved in immune response,” “activation of MAPK activity,” and “positive regulation of NF-kappaB import into nucleus.” The leukocyte transendothelial migration is one of the important steps of the innate immune system in triggering the inflammatory immune response and the migration of the first immune response cells such as neutrophils to the sites of infection (
Coexpressed genes in the greenyellow module were enriched in KEGG pathways and biological process such as “T-cell receptor signaling pathway,” “focal adhesion,” “leukocyte transendothelial migration,” and “regulation of cell motility.” These pathways have also been observed in previous transcriptomic studies in response to infection with the agents of BRDC (
Functional enrichment analysis revealed that the purple module was enriched in the “NOD-like receptor signaling pathway,” “C-type lectin receptor signaling pathway,” “RIG-I-like receptor signaling pathway,” and “necroptosis” KEGG pathways as well as “inflammatory response,” “type I interferon signaling pathway,” “positive regulation of type I interferon production,” “interferon-gamma-mediated signaling pathway,” and “cellular response to interferon-gamma” biological processes. RIG-I-like receptors are important cytoplasmic PRRs that detect intracellular viruses through their genomic RNA (
Several recent studies have shown that the type I interferon signaling pathway has been identified as a key pathway in cows with BRD atfeedlot entry, which indicates that animals show antiviral responses at the entry stage (
In agreement with the previous studies, functional terms including “VEGF signaling pathway,” “T-cell receptor signaling pathway” (
The yellow module showed that its genes were enriched in some important biological processes and KEGG pathways associated with BRD such as “neutrophil activation involved in immune response,” “regulation of inflammatory response,” “positive regulation of T-cell proliferation,” “MAPK signaling pathway,” “NF-kappa B signaling pathway,” “apoptosis,” “TNF signaling pathway,” “JAK-STAT signaling pathway,” and “Cytokine-cytokine receptor interaction.” Furthermore, examination of the relationship between the hub–hub genes of this module and BRD showed that many of these genes, such as
These findings demonstrate the relevance of the mentioned modules as well as their genes, especially hub–hub genes and TFs as important candidates in the development of BRD, helping us to better understand the molecular mechanisms responsible for the immune response to BRD. Further researches are needed to more closely examine the biological behavior and functions of these modules and their genes during BRD.
Given that BRD is the main cause of morbidity and mortality in beef and dairy cattle and has a potential impact on economic losses in the livestock industry, a systems biology approach was used to further investigate the molecular mechanisms of BRD as well as to identify diagnosis biomarkers and therapeutic targets for BRD. In this study by using WGCNA distinct methods (MTRs and MP) and functional enrichment analysis, we identified eight candidate modules that are involved in the immune response and BRD pathogenesis. It is noteworthy that both WGCNA methods showed a similar ability to identify candidate modules during BRD, confirming each other results. Integrated coexpressed hub genes of eight candidate modules with PPI networks, allowed us to identify hub–hub genes that act as central genes in both coexpression and PPI networks and may be important candidates during BRD development. In total, we identified 307 hub–hub genes in eight candidate modules, most of which were potentially involved in BRD. These genes along with other members of the eight candidate modules could be important targets in the pathogenesis of BRD for future researches. Therefore, more research is needed to validate the hub–hub genes reported in this study, especially those whose role in the immune system in response to BRD is still unclear.
The original contributions presented in the study are included in the article/
AH conceived the ideas. AH designed the study. AH and NS analyzed the data. AH, AB, and FF interpreted the data. AH wrote the main manuscript. GJ, MB, HK, and SI helped in writing the manuscript. AB and RA reviewed and edited the manuscript. HM, AB, and HB reviewed the manuscript. All authors read and approved the final version of manuscript.
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.
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.
We would like to thank the University of Tehran for supporting this research.
The Supplementary Material for this article can be found online at:
BCV, bovine coronavirus; BHV-1, bovine herpes virus type 1; BPIV-3, bovine parainfluenza type 3 virus; BRD, bovine respiratory disease; BRSV, bovine respiratory syncytial virus; BVDV, bovine viral diarrhea virus; GEO, Gene Expression Omnibus; GO, Gene Ontology; GS, gene significance; KEGG, Kyoto Encyclopedia of Genes and Genomes; ME, module eigengene; MM, module membership; MTRs, module–trait relationships; MP, module preservation; PPI, protein–protein interaction; RSV, respiratory syncytial virus; RNA-seq, RNA-sequencing; TF, transcription factor; TOM, topological overlap matrix; WGCNA, weighted gene coexpression network analysis.