RNA-Seq Bayesian Network Exploration of Immune System in Bovine.

BACKGROUND: The stress is one of main factors effects on production system. Several factors (both genetic and environmental elements) regulate immune response to stress.
OBJECTIVES: In order to determine the major immune system regulatory genes underlying stress responses, a learning Bayesian network approach for those regulatory genes was applied to RNA-Seq data from a bovine leukocyte model system.
MATERIAL AND METHODS: The transcriptome dataset GSE37447 was used from GEO and a Bayesian network on differentially expressed genes was learned to investigate the gene regulatory network.
RESULTS: Applying the method produced a strongly interconnected network with four genes (TERF2IP, PDCD10, DDX10 and CENPE) acting as nodes, suggesting these genes may be important in the transcriptome regulation program of stress response. Of these genes TERF2IP has been shown previously to regulate gene expression, act as a regulator of the nuclear factor-kappa B (NF-κB) signalling, and to activate expression of NF-κB target genes; PDCD10 encodes a conserved protein associated with cell apoptosis; DDX10 encodes a DEAD box protein and is believed to be associated with cellular growth and division; and CENPE involves unstable spindle microtubule capture at kinetochores. Together these genes are involved in DNA damage of apoptosis, RNA splicing, DNA repairing, and regulating cell division in the bovine genome. The topology of the learned Bayesian gene network indicated that the genes had a minimal interrelationship with each other. This type of structure, using the publically available computational tool, was also observed on human orthologous genes of the differentially expressed genes.
CONCLUSIONS: Overall, the results might be used in transcriptomic-assisted selection and design of new drug targets to treat stress-related problems in bovines. Copyright:

1. Background

The physiological stress-induced immune response in the bovine could happen from many different circumstances, including injury, calving, weaning, dry period, cell-mediated destruction of pathogens, stress, failure of the mammary glands defense mechanism and mastitis, and other sources of stress. Some of this physiological immune response, such as mastitis response, is of vital important in bovine milk production, but is hardly possible to tackle using quantitative genetic theory. For instance, the estimated heritability of mastitis is quite low (0.01 to 0.17 in different references). Therefore, disorders of immune responses affecting many dangerous and costly diseases in cattle should be well addressed. One approach would be integrative systems biology methods, using OMICS data from RNA-Seq transcriptomics to explore the molecular networks underlying immune response mechanisms. RNA-Seq is a novel sequencing technology generating detailed information on gene expression ( 1 ). In the context of bovine transcriptomics studies, RNA-Seq has been used in various areas, e.g. detection of novel splice variants in Zebu cattle to cure horn cancer ( 2 ), transcriptome profiling to study growth and development of muscle in Chinese Luxi an

Angus beef cattle ( 3 ), the stress response to weaning in bovine leukocytes ( 4 ), and transcriptional profiling of peripheral blood leukocytes from cattle infected with mycobacterium bovis ( 5 ). The Bayesian network (BN) is an attractive formalism that could capture gene regulatory properties and conditional probabilistic independence among genes. This formalism reduces the parameter space search over the domain of variables. In this way, considering k variables (X1, …, Xk) then, the notion of:

Leads to a dramatic decrease in the number of

parameters over parameter space. BN is a strong method that is able to learn and capture linear/nonlinear, combinatorial and stochastic relationships among many variables ( 6 ). BN is being used in modelling immune response ( 6 ), mastitis management on dairy farms ( 7 ), presence of claw and digital skin diseases ( 8 ), to study simple gene regulatory networks of immune system candidate genes in dairy cattle ( 9 ), and to estimate inference of gene regulatory network from RNA-Seq time series data ( 10 ). Over several years, much attention has been given to improve production traits in cattle, which might have (in) directly impaired the immune system. Recently, genetic selection has been used to improve the immune system in dairy cows ( 11 ).

A thousand genes (8-9% of the genome) are responsible for regulation of the immune system in mammals ( 11 ).This number of genes poses many challenges to researchers and breeders before it would be possible to produce an animal with superior immune responses and benefit the farming enterprise.

In this study, it has been hypothesized that BN algorithm based on conditional probabilities could find out gene interactions alike in cellular form. To uncover the regulation scenario, data from bovine transcriptomic leukocytes RAN-Seq were used in a learning Bayesian network approach to find causative stress-related genes underlying the immune system. Constructing the Bayesian network on such data might reveal major genes and biomarkers controlling the immune system, in pathways that would attractive for interventions. By capturing regulator genes on the transcriptome of bovine's immune system, it could be likely possible to single out genes with pivotal effects on general health performance of cattle. They could have a plausible application in constructing a promising breeding program.

2. Objectives

The re-use of previously published data due from transcriptomics technologies is becoming commonplace to derive new hypothetical motivated questions. Newly generated experimental data is often not needed to generate good biological questions. In this study, using previously issued data, we tried to draw out some regulatory genes performed in bovine immune systems. Also, to enrich the results of this study, we have tried to map the set of differentially expressed genes to their orthologous human counterpart genes. The transfer of such animal-based knowledge to the human application can be very beneficial.

3. Materials and Methods

3.1. RNA-Seq Dataset Processing

The transcriptome dataset was used from GEO under accessionnumberGSE37447 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc =GSE37447). This data was respect to RNA-Seq leukocyte of male beef calves in response to weaning stress over different days (0, 1, 2 and 7 d). The whole description of this dataset can be found in ( 4 ). The data quality control, to assess possible biases in the RNA-Seq data, was done by FastQC ( 12 ). Trimming and filtering poor bases quality (quality score < 20) of RNA-Seq data was done by Trimmomatic software ( 13 ). Alignment of reads to bovine reference genome (asia.ensemble.org) was performed by tophat2 (v 2.0.9) ( 14 ). Covering aligned read to counts per gene was accomplished by HT-Seq ( 15 ). Merging HT-Seq was done to make a gene expression matrix with at least 3 counts in every sample. Differentially expressed genes were obtained with the edgeR package ( 16 ) that evaluates the differential expression in read counts of RNA-Seq data by empirical Bayesian approaches. Adjusted P-value was set to 0.05 as a threshold for selecting genes with differential read counts during the time point between samples. Then, normalization was performed using the Limma package ( 17 ). The package was originally designed for the analysis of microarray data, but it has been extended to the analysis of RNA-Seq data in the form of normalized log2-transformed counts by adding a new normalization function termed voom. The voom transformation converts the counts to log-counts per million with the associated precision weights ( 18 ).

3.2. RNA-Seq Bayesian Network

To investigate the gene regulatory network, a Bayesian network on differentially expressed genes was learned using the networkBMA package ( 19 ). To accomplish the analysis, the matrix of differentially expressed genes was transposed, columns and rows represented genes and observations at different time points, respectively. Cytoscape software was used for visualization of the learned regulatory network and extracting graph- theoretic measures ( 20 ).

To further enrich the results, we sought human orthologue genes of the most connected genes found in this study. To do so, the most connected genes to http://www.esyn.org were imported to get the genetic and physical interactions and also the graph theoretic measure of those genes with themselves and with other genes deposited in the H. sapiens (BioGrid) database ( 21 ). To find the gene pathway and ontology of the most connected genes, the InnateDB Pathway Analysis software (www.innatedb.com/batchSearchInit.do) was used ( 22 ). Also, the gene pathway analysis and gene ontology enrichment for whole differentially expressed genes was found by DAVID software ( 23 ). The pipeline of analysis to discover cellular mechanism and major genes associated with the activation of the immune system induced by weaning stress is presented in Figure 1.

Figure 1

Pipeline of RNA-Seq data processing in current study.

4. Results

The details of time series RNA-Seq data, such as the numbers of raw reads, trimmed reads, and mapped reads to bovine genome are provided in Table S1.

In the current study, 17817 genes were detected in the GSE37447 experiment (by tophat2). A summary of differentially expressed genes is given in Table S2. Differentially gene expression analysis detected 220 genes, which have differentially expression level between stress and un-stress calve groups.

According to the RNA-Seq Bayesian gene regulatory network developed in the present study, which was constructed on 220 differentially expressed genes, it was concluded that four genes TERF2IP, PDCD10, DDX10 and CENPE (Fig. 2) were the most connected hub genes, suggesting a significant regulatory effect on other immune-related genes. Statistical parameters of BN showed in Table 1. Table 2 represents the length and number of shortest distances between two connected nodes, and Table 3 represents connected components’ features. It is highly likely that the extracted hub genes in this study would be involved in the immune response and prevent large-scale development of inflammation that may lead to tissue damage (see: Table S3 for further hub genes’ enrichment analysis). Also, the hub genes improve immune function by regulating RNA splicing, DNA repairing, influencing cell cycle and cell proliferation and have an inhibitory effect by apoptosis. The results imply that the genes could have a central role in other types of immune responses because they usually have conserved sequences. By mapping the differentially expressed genes to orthologous genes in H. sapiens (BioGrid) database repository (http://www.esyn.org/ ) ( 24 ), some network measures (Table 4) and their topology (Fig. 3) were obtained. Some parts of the obtained network were shown in Table 4 (the full results can be seen in Table S4).

Figure 2

RNA-Seq bayesian network visualization by Cytoscape. The biological importance of nodes in the network is identified by color (high importance effect to low importance was represented by red to green) and node size (major nodes were represented by larger size).

Table 1

Statistical parameters of Bayesian gene regulatory network

Statistical parameters
	Density	Connectivity	Diameter	Radius	Shortest path	Clustering coefficient
Bayesian network	0	4.38	6	1	1343	0.160024

Table 2

Structure of network paths

Shortest path length	1	2	3	4	5	6
Number of shortest path	424	564	203	83	67	2

Table 3

The network parameters of the most connected genes in GSE37447 experiment

Gene (Ensemble ID)	Gene name	Clustering coefficient	Out- degree	In-degree	Neighbors connectivity	Betweenness centrality	Closeness centrality
ENSBTAG00000015686	TERF2IP	0.003	155	5	3.069	0.016	0.852
ENSBTAG00000031709	PDDCD10	0.003	94	3	3.546	0.009	0.446
ENSBTAG00000002382	DDX10	0.019	68	0	6.224	0.000	0.601
ENSBTAG00000009035	CENPE	0.001	64	0	2.635	0.000	1.000

Table 4

The network parameters of most connected genes inferred from BioGrid repository

Gene	Degree	Radiality	Closeness	Stress	Betweenness	Centroid value	Eccentricity	Collective influence
STK24	1	1.203846	0.025281	0.000000	0.000000	-204	0.250000	0
PDCD10 *	53	1.988462	0.035433	0.476513	0.476513	-100	0.333333	7852
STK25	1	1.203846	0.025281	0.000000	0.000000	-204	0.250000	0
AP2B1	1	1.965385	0.035019	0.000000	0.000000	-204	0.250000	0
TERF2IP *	153	2.750000	0.058065	1.000000	1.000000	98	0.333333	7852
BUB1B	1	0.407692	0.26087	0.000000	0.000000	-34	0.50000	0
CENPE *	36	0.538462	0.514286	0.030546	0.030546	34	1.000000	0
FGFR1OP2	1	1.203846	0.025281	0.000000	0.000000	-204	0.250000	0
TRAF3IP3	1	1.203846	0.025281	0.000000	0.000000	-204	0.250000	0
G3BP2	1	0.211538	0.514286	0.000000	0.000000	-17	0.500000	0
DDX10 *	18	0.276923	1.000000	0.007855	0.007855	17	1.000000	0
SLX4	1	1.965385	0.035019	0.000000	0.000000	-204	0.250000	0
PRC1	1	0.407692	0.26087	0.000000	0.000000	-34	0.500000	0
MAPK1	1	0.407692	0.26087	0.000000	0.000000	-34	0.500000	0

The most connected genes

Figure 3

Topological view of the most connected genes (PDCD10, TERF2IP, CENPE and DDX10) with themselves and other genes in BioGrid.

It was a great surprise that those genes that turned up as hub genes (TERF2IP, PDCD10, DDX10, CENPE) in the GSE37447 experiment, were also the hub genes in current mapping results in H.Sapiens (BioGrid). This culminated in a similar topology (comparing the Fig. 2 and Fig. 3). The results of gene ontology for most connected genes in the innate database can be seen in Table S3. Gene pathway analysis showed that spliceosome and mismatch repair pathway were two over- represented pathways in the weaned group versus control group (P < 0.01) (Table S3).

5. Discussion

It is reported that weaning stress could increase the ability of the immune system to identify and clear pathogens by affecting related pathways such as cytokine signaling, G-protein-coupled receptor (GPRC) signaling, transmembrane transport and homeostasis ( 4 ). Therefore, wisely extracting knowledge from current data is crucial to mining regulator genes, which are involved in the mechanisms of the immune system activation (Fig. 1). In this study, 17817 genes were detected from related dataset (by tophat2), which is more than O’Loughlin et al. ( 4 ), that used Bowtie as aligner and reported 16514 genes which summarize in Table S2.

Network analysis of constructed BN is detailed on Table 1. Network connectivity refers to the number of nodes that are pair-wise connected (connected components). Network diameter and radius show the maximum and minimum distance between two nodes, respectively. In the current network, the diameter was 6 and the radius was calculated as 1. Shortest path parameter obtained was 1343, which shows all the shortest distances between nodes in the network. The network clustering coefficient (0.16) describes the average of all nodes clustering coefficient. A substantial change in network research, currently facilitated by improved computer networks, has recently targeted the extraction of the statistical properties of large-scale networks (over billion vertex). Biological networks are cumbersome to envision and describe without a set of network metric statistics or quantitative measures. For example, Behdani and Bakhtiarizadeh ( 25 ) suggested an integrated gene regulatory network using module inference to construct modules of co-expressed genes with bovine leukocyte RNA-Seq data and assigning transcription factors to these modules using Lemon- Tree algorithms. They identified two transcription factors E2F8 and FOXS1 as novel regulatory candidates in immune response. Even though a through account of possible application of network theoretical exaptation and application in biological sciences can be seen in Pavlopoulos et al. ( 26 ), but still lots of studies need to be performed to pinpoint which network measures are of greatest priority in which biological context. Hub genes network measures are shown in Table 3. It could be seen that DDX10 and CENPE have 0 in-degree value, suggesting they could be the sole regulator in this study. Statistical parameters of BN showed in Table 1. Table 2 represents the length and number of shortest distances between two connected nodes, and Table 3 represents connected components’ features. Node clustering coefficient implies the ratio of the number of edges between the neighbours of each node to the maximum number of edges that could possibly exist between the neighbours of same nodes. Out-degree of a given node reveals the number of edges that coming out of node and the in-degree parameter displays the number of edges that are entering the node. Out-degree of current network ranged between 0 and 155 and in-degree nodes/genes changed between 0 and 4. A neighbour's connectivity for a given node/gene shows the average connectivity of all its neighbours, and in the current study varied from 111.99 to 2.63. The betweenness centrality parameter of a node states the amount of node regulatory effect on other nodes in the network (Yoon et al. 2006), and it changes between 0 and 1. Yoon et al. ( 27 ) reported that betweenness centrality of a given node in a graph describes its influence on other nodes in a network. In the current study, betweenness varied between 0 and 0.016. Closeness centrality depends on average of shortest path between a given node with neighbours, and can be between 0 and 1. Newman ( 28 ) stated that the closeness centrality of each node refers to the magnitude of the influence of neighbours on a given node. The higher values of closeness centrality, the higher regulatory effects of a given node on its neighbours. According to this definition a regulatory gene (TERF2IP, PDCD10, DDX10 and CENP-E) must have a higher value of the parameter. The higher the values of closeness centrality for these genes, the stronger regulatory effects on other nodes.

TERF2IP (Telomeric repeat binding factor 2 interacting protein) could be having an important role in weaning transcriptome network. This gene and several other genes (TERT, POT1, TNKS, TERF1, TINF2 and TERF2) that expressed in telomeres have less nucleotide diversity than other gene families. Reports showed TERF2IP has an alternative effect on the immune system and inflammation by NF-κB pathway, so that inhibition of TERF2IP leads to decreases in pro- inflammatory factors in mesenchymal stem cells (MSCs) ( 29 ). Nuclear factor-kappa B (NF-κB) signaling pathways affect the native and adaptive immune systems, apoptosis, cell cycle, cell differentiation and migration ( 30 ). Therefore, TERF2IP causes an appropriate immune response and controls inflammation by activation of NF-κB ( 31 ). In addition, TERF2IP controls apoptosis through NF-κB led to regulating immune responses and preventing development of immune functions ( 32 ). Gene ontology of this regulator gene approves its proven function related to regulation of NF-κB (Table S3).

The second important gene in this study was PDCD10 (programmed cell death 10). It is also called CCM3 (cerebral cavernous malformation 3). In oxidative conditions, activation of ezrin/radixin/moesin (ERM) protein family can help cell survival. ERM proteins involved in apoptosis, cell adhesion and migration by connection between cAMP signalling pathway and coupled-G protein receptors ( 33 ). It was demonstrated that PDCD10 is necessary to cell viability under stress conditions by activation of ERM protein family. Fidalgo et al. ( 34 ) showed inhibition of PDCD10 led to inefficiency of ERM phosphorylation and made the cell more sensitive to stress condition. It seen that PDCD10 improved immune function by effect on T-cell and leukocyte viability. Gene ontology of this hub gene showed its role in improvement of the immune response by positive affecting on proliferation and negative affecting on apoptotic process (Table S3). According to current study, DDX10 is an important gene with important functions in the immune system under stress conditions. DEAD (Asp-Glu-Ala-Asp) box polypeptide 10 and HRH-J8 are other names for it. DDX10 is an ATP-depended RNA helicase involved in initiation of transcription, RNA splicing, ribosome and spliceosome assembly and mRNA stability ( 35 ). As noted earlier, splicing process is necessary for optimal immune functions. Helicases makes this process easier by facilitating the pre-RNA joining and separating to snRNA ( 36 ). According to results of Bayesian network, increased expression of DDX10 as a regulator in immune response indicate that it can improve immune function by affecting RNA metabolism and RNA splicing.

The last gene that influences immune responses based on current results was CENP-E. The protein was translated from this gene temporarily presents on centromere at special times of cell cycle and leads to proper alignment between homologue chromosomes due to controlling interactions between kinetochores and spindles during mitotic division ( 37 , 38 ). The role is critical for cell survival, as previous studies reported that the inhibition of this gene can lead to apoptosis ( 37 ). During immune response, division of leukocytes occurs at a high rate. Therefore, it is necessary to efficient cell cycle regulation that it does not lead to some immune disorders such as cancer or autoimmunity and establishment of immune cells ( 39 ).

One of the control steps that regulates cell division occurs on metaphase during pairing of homologue chromosomes, therefore creating of a signal to delay the entrance to anaphase occurs in order to verify the accuracy of this alignment. It is well known that CENP-E is one of the proteins that influences the function of check points of the cell cycle and affects the signals for delayed entrance to anaphase ( 40 ). Increased expression of CENP-E by affecting accuracy of leukocyte proliferation can probably improve immune responses. These established roles of this major gene in the Bayesian network is confirmed by significant terms of its gene ontology (Table S3).

Our results indicated spliceosome and mismatch repair pathway were two over-represented pathways in comparison of weaned group versus control group (Table S3). The spliceosome is a large multi-subunit protein and RNA complex that facilitates intron separating, and alternative splicing refers to a process in which several different transcripts from a pre-transcript are produced by spliceosome complex ( 41 ). In addition, some reports have shown that genes involved in different aspects of the biological T-cell functions, have been rich signals that associated with alternative splicing ( 42 ). Activation of the DNA damage pathway may represent a more distinctive feature of oxidative stress in livestock that was induced by many factors such as weaning stress ( 43 ). In stress conditions, neutrophils and lymphocytes increase respiratory burst, which leads to increased reactive oxygen species and the creation of oxidative stress ( 44 ). The mismatch repair system is the main post-replicative pathway for the correction of replication errors that are not corrected by proofreading. Some receptors were shown that simulation of the innate immune system was up- regulated in the DNA damage process. Thus, DNA damage responses such as the mismatch repair mechanism can active the innate immune system. In addition, some evidence has shown some pro- inflammatory factors can accelerate the DNA repair activity in cells during inflammation ( 31 ). Some reports have found that conditions that involve damage to DNA, lead to accelerated immune system response by increasing the expression of ligands for NKG2D receptors ( 45 ). Reports stated some inflammatory modulators such as NF-κB can regulate DNA repair process during the immune responses. Bacterial or viral productions, oxidative stress, and pro-inflammatory cytokines such as IL1 and TNF-α could play a role as signals to activate the NF-κB pathway and to affect the immune responses, apoptosis, inflammation and DNA repair process ( 46 ). Results of gene ontology of differentially expressed genes have shown cellular response to stress, cellular response to DNA damage stimulus, mRNA splicing via spliceosome, RNA splicing, and cell cycle were the most significant terms (Table S5). These results matched with gene pathway analysis and functions of hub genes in the Bayesian regulatory network.

6. Conclusions

According to the Bayesian gene regulatory network developed here, it is concluded that four genes/nodes (TERF2IP, PDCD10, DDX10 and CENP-E) had regulatory effects on other immune related genes. Based on the results, these modulation of these genes could improve immune response and prevent large-scale development of inflammation that may lead to tissue damage. Also, they improve immune function by regulating RNA splicing, DNA repairing, influencing cell cycle and cell proliferation and have an inhibitory effect by apoptosis. The results showed these genes have a central role in other species immune responses, because usually have conserved sequences. Pathway analysis shown weaning stress could damage DNA by creating oxidative conditions and it leads to activating DNA repair mechanisms. DNA damage not only affects DNA repair mechanisms, but also actives immune responses and releases inflammatory mediators that leads to the involvement of spliceosome pathway. In addition, inflammatory mediators can directly affect and enhance the DNA repair mechanisms.