RNA-Seq Analyses Reveal That Endothelial Activation and Fibrosis Are Induced Early and Progressively by Besnoitia besnoiti Host Cell Invasion and Proliferation
Alejandro Jiménez-Meléndez1 
Chandra Ramakrishnan2 
Adrian B. Hehl2 
Giancarlo Russo3 
Gema Álvarez-García1*- 1SALUVET, Animal Health Department, Faculty of Veterinary Sciences, Complutense University of Madrid, Madrid, Spain
- 2Institute of Parasitology, University of Zurich, Zurich, Switzerland
- 3Functional Genomics Center Zurich, Zurich, Switzerland
The pathogenesis of bovine besnoitiosis and the molecular bases that govern disease progression remain to be elucidated. Thus, we have employed an in vitro model of infection based on primary bovine aortic endothelial cells (BAEC), target cells during the acute infection. Host-parasite interactions were investigated by RNA-Seq at two post-infection (pi) time points: 12 hpi, when tachyzoites have already invaded host cells, and 32 hpi, when tachyzoites have replicated for at least two generations. Additionally, the gene expression profile of B. besnoiti tachyzoites was studied at both pi time points. Up to 446 differentially expressed B. taurus genes (DEGs) were found in BAEC between both pi time points: 249 DEGs were up-regulated and 197 DEGs were down-regulated at 32 hpi. Upregulation of different genes encoding cytokines, chemokines, leukocyte adhesion molecules predominantly at 12 hpi implies an activation of endothelial cells, whilst upregulation of genes involved in angiogenesis and extracellular matrix organization was detected at both time points. NF-κB and TNF-α signaling pathways appeared to be mainly modulated upon infection, coordinating the expression of several effector proteins with proinflammatory and pro-fibrotic phenotypes. These mediators are thought to be responsible for macrophage recruitment setting the basis for chronic inflammation and fibrosis characteristic of chronic besnoitiosis. Angiogenesis regulation also predominated, and this multistep process was evidenced by the upregulation of markers involved in both early (e.g., growth factors and matrix metalloproteinases) and late steps (e.g., integrins and vasohibin). Besnoitia besnoiti ortholog genes present in other Toxoplasmatinae members and involved in the lytic cycle have shown to be differentially expressed among the two time points studied, with a higher expression at 32 hpi (e.g., ROP40, ROP5B, MIC1, MIC10). This study gives molecular clues on B. besnoiti- BAECs interaction and shows the progression of type II endothelial cell activation upon parasite invasion and proliferation.
Introduction
Besnoitia besnoiti is the ethiological agent of bovine besnoitiosis (Besnoit and Robin, 1912), a re-emerging disease in Europe with a progressive dissemination in beef cattle herds and negative impact in cattle welfare and fertility (European Food Safety Authority, 2010; Cortes et al., 2014). This parasitic disease is responsible for both cutaneous and systemic clinical signs, as well as sterility in bulls (Álvarez-García et al., 2014). Disease initiates with the acute stage, when the tachyzoites are fast-replicating in endothelial cells, and evolves with the chronic stage, characterized by the development of bradyzoite-containing tissue cysts located mainly in the subcutaneous tissue and mucous membranes. Acutely infected animals can develop oedemas, orchitis, respiratory distress, and may suddenly die due to a cardio-respiratory failure. Chronically infected cattle develop pathognomonic tissue cysts responsible for skin lesions and the initial orchitis may end up with testis atrophy (reviewed by Álvarez-García et al., 2014).
Endothelial cells (ECs) are one of the main target cells for parasite replication during the acute infection (McCully et al., 1966), and their dysfunction may be determinant for tissue damage, particularly in testes and lungs. Vascular damage has been histologically described mainly in skin, characterized by degenerative and fibrinoid necrotic vascular lesions, vasculitis and thrombosis in medium-sized and small veins, as well as small arteries (Pols, 1960; McCully et al., 1966; Basson et al., 1970; Bigalke et al., 2004). The molecular bases that underlies these host-parasite interactions remain to be elucidated. In vitro studies performed so far in calf umbilical vein endothelial cells (BUVEC) (Maksimov et al., 2016; Taubert et al., 2016) showed that B. besnoiti infection results in an increase in the transcripts of key genes such as P-selectin, intercellular adhesion molecule 1 (ICAM-1), chemokines (CXCL1, CXCL8, CCL5), IL-6 and COX-2 related to endothelial cell activation and leukocyte recruitment.
As obligate intracellular pathogens, protozoan parasites belonging to Toxoplasmatinae subfamily (subphylum Apicomplexa) have developed unique adaptations that allow them to modify the host cellular machinery to alter host gene expression and favor their replication and dissemination known as lytic cycle (Blader et al., 2015; Hakimi and Bougdour, 2015). However, time points of adhesion/invasion, proliferation and egress of B. besnoiti tachyzoites differ from the ones reported in other Toxoplasmatinae parasites, since tachyzoites showed a delayed invasion of host cells, followed by a lag phase of up to 24 hpi, as well as the capacity to survive for up to 24 hpi when still extracellular (Frey et al., 2016). Moreover, species-specific differences in the protein repertoire involved have been reported between Toxoplasma gondii and Neospora caninum (Reid et al., 2012). Proteins contained in specific apicomplexan organelles such as rhoptries (ROP), dense granules (GRA) or micronemes (MIC) (Blader et al., 2015) play a key role at a transcriptional level altering many phenomena (e.g., blockage of host immune pathways, activation of transcription factors, modifications on the chromatin and small non-coding RNA) (Hakimi et al., 2017). In B. besnoiti, the identification of ortholog genes has been hampered by the lack of the whole genome sequence (García-Lunar et al., 2013, 2014), which has only been recently announced (Schares et al., 2017).
Herein, we aimed to dissect the molecular mechanisms that govern the endothelial injury during acute B. besnoiti infection in a recently standardized in vitro model based on primary target bovine aorta endothelial cells (BAEC) by means of RNA-Seq. Remarkably, this cell line was free from bovine viral diarrhea virus (BVDV), which is widely known to alter the transcriptomic profile of endothelial cells in vitro (Neill et al., 2008) and to facilitate early B. besnoiti invasion (Jiménez-Meléndez et al., 2019). Transcriptomics analyses were carried out at two post-infection (pi) time points, representative of early invasion (12 hpi) and intracellular proliferation (32 hpi). In parallel, a repertoire of B. besnoiti tachyzoite differentially expressed genes (DEG) involved in the lytic cycle were also identified.
Materials and Methods
Parasites, Cell Cultures, and Experimental Design
The primary bovine aortic endothelial cells (BAEC) were maintained in M200 medium (Gibco, Life Technologies, Thermo Fisher Scientific, Waltham, MA, USA) containing 20% of FBS, 100 IU/mL of penicillin, 100 mg/mL of streptomycin, 0.25 μg/mL amphotericin B and Large Vessel Endothelial Growth supplement as described by Jiménez-Meléndez et al. (2019). Tachyzoites from the Bb-Spain 1 isolate of B. besnoiti were routinely maintained in a monolayer culture of the monkey kidney MARC-145 cell line as previously described (Jiménez-Meléndez et al., 2017). As quality controls, the absence of Mycoplasma spp. and BVDV were checked according to Jiménez-Meléndez et al. (2019).
Tachyzoites were harvested at 3 days post infection (dpi), when the majority of parasites were still intracellular, by recovering the infected cell monolayer with a cell scraper, followed by repeated passages through a 25-gauge needle at 4°C and separation from cell debris on a PD-10 column (Frey et al., 2016). Tachyzoite viability was confirmed by trypan blue exclusion followed by counting in a Neubauer chamber. Purified tachyzoites were used to infect BAEC confluent cell monolayers as described below.
Infections in BAEC for transcriptomics analyses were performed as previously described (Jiménez-Meléndez et al., 2019). T25 flasks seeded with confluent BAEC monolayers (2 × 106 cells) were infected with purified tachyzoites at a multiplicity of infection (MOI) of 10. Preliminary experiments were performed in order to optimize the infection rate by increasing MOIs (5, 10, 20). The maximum invasion rate for the isolate employed was 20% (unpublished data). Infected flasks were washed at 12 hpi, when most parasites have already invaded but have not yet started replicating, and at 32 hpi, when parasites have replicated at least twice (Jiménez-Meléndez et al., 2019). Cells and parasites were recovered using a cell scraper in 5 ml of cold PBS, pelleted by centrifugation at 1,350 × g for 10 min, resuspended in 300 μl of RNAlater® (Thermo Fisher Scientific, Madrid, Spain) and stored at −80°C until RNA extraction. In parallel, T25 flasks containing non-infected BAECs were similarly collected at 12 hpi. All analyses were performed with four biological replicates and, for each replicate, technical replicates were collected for further validation by quantitative real time PCR (qPCR).
RNA Extraction
Total RNA from the cell culture samples was extracted using QIAGEN RNeasy Mini Kit following a QIAshredder® homogenization according to the manufacturer's instructions. Briefly, samples in RNAlater® were lysed in 900 μL of RLT buffer + β-mercaptoethanol at 10 μl/mL and column purified; an on-column DNase I digestion was included. Total RNA was eluted in RNase-free water. The RNA concentration and purity were assessed spectrophotometrically at 260 nm using NanoDrop 1000 (Implen, Munich, Germany).
Quality Control of Total RNA, Library Preparation, and RNA-Seq
The quality and quantity of the total RNA were further determined in an Agilent TapeStation 4200 to determine the RNA integrity Number (RIN). Only RNA extracts with concentrations between 0.5 and 5 μg/mL, RINs higher than 9.6, and with 260/280 ratios between 1.8 and 2.0 were included in the study. The TruSeq RNA Sample Prep Kit v2 (Illumina, Inc., California, USA) was used in the succeeding steps. Briefly, mRNA was enriched from total RNA samples (100–1,000 ng) using poly-A selection before transcription into double-stranded cDNA. The cDNA samples were fragmented, end-repaired and polyadenylated before ligation of TruSeq adapters containing the index for multiplexing Fragments containing TruSeq adapters on both ends were selectively enriched with PCR. The quality and quantity of the enriched libraries were validated using the Qubit® (1.0) Fluorometer and the Caliper GX LabChip® GX (Caliper Life Sciences, Inc., USA). The product is a smear with an average fragment size of ~260 bp. The libraries were normalized to 10 nM in 10 mM Tris-Cl/0.1% Tween 20, pH8.5.
The TruSeq PE Cluster Kit v3-cBot-HS or TruSeq SR Cluster Kit v3-cBot-HS (Illumina Inc.) was used for cluster generation using 10 pM of pooled normalized libraries on the cBOT. Stranded sequencing was performed on the Illumina HiSeq 4000 with single end sequencing at 125 bp using the TruSeq SBS Kit HS4000 (Illumina Inc.).
Computational Analysis of RNA-Seq Data
Reads were quality-controlled with FastQC. Sequencing adapters were removed with Trimmomatic (Bolger et al., 2014), and reads were hard-trimmed by 5 bases at the 3' end. Successively, reads at least 20 bases long, and with an overall average phred quality score >10 were aligned to the reference genome and transcriptome of Bos Taurus (UCSC/bosTau7) and B. besnoiti (BbLIS14 in-house assembled genome, unpublished) with STAR v2.5.1 (Dobin et al., 2013), with default settings for paired end reads.
Distribution of the reads across genomic isoform expression was quantified using the R package GenomicRanges (Lawrence et al., 2013) from Bioconductor Version 3.0. Differentially expressed genes were identified using the R package DESeq2 (Love et al., 2014) from Bioconductor Version 3.0. DESeq2 performs an internal normalization, where the geometric mean is calculated for each gene across all samples. The counts for a gene in each sample is then divided by this mean. The median of these ratios in a sample is the size factor for that sample. We used a threshold of ≥1.5fold difference (fold change) in mRNA levels (measured as normalized averaged mapped sequence reads per gene) and a false discovery rate (FDR) ≤ 0.05 in order to categorize genes as differentially expressed (DE).
Functional Enrichment and Network Analysis
Functional enrichment was applied to all comparisons performed using the genome from Bos taurus available at the Panther website (http://www.pantherdb.org/) and the statistical overrepresentation test, using an FDR of 0.05 as threshold. The outcome of the gene ontology (GO) analysis was used for semantic clustering using REVIGO in order to identify non-redundant GO terms (Supek et al., 2011). The analyses were performed with an allowed similarity of 0.7 both for up-regulated and down-regulated genes at both infection time points, and with default settings in advanced options.
Furthermore, STRING 10 version 11.0 (Szklarczyk et al., 2015) was used to identify interesting associations between the significant genes identified in our study and to construct Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa, 2019). Using the STRING database (http://string-db.org/), multiple proteins were chosen from the website interface setting the maximum relevance to 0.9. The differentially expressed genes (DEGs) names were inserted as the input in the list of names, and Bos taurus was chosen as the organism.
Further annotation and GO assignment in Besnoitia besnoiti were performed using the Blast2GO platform and NCBI's nr protein database. Blast2GO assigns descriptions using a language processing algorithm to hit descriptions, which extracts informative names and avoids low-content terms (Conesa and Gotz, 2008). These descriptions were manually curated where necessary.
Quantitative Real-Time PCR (qPCR) for Transcriptome Validation
Reverse transcription was performed using the master mix SuperScript® VILO™ cDNA Synthesis Kit (Invitrogen, Paisley, UK) in a 20 μL reaction using up to 2.5 μg of total RNA. cDNA was sequentially diluted to 1:20, 1:80, 1:320 and 1:1,280, and all dilutions were analyzed by quantitative real-time PCR (qPCR). Quantitative real-time PCRs were performed in 25 μL volumes using 12.5 μL of Power SYBR®PCR Master Mix (Applied Biosystems, Foster City, CA, USA), 10 pmol of each primer and 5 μL of the diluted cDNA samples. Reactions were performed in an ABI 7500 FAST Real Time PCR System (Applied Biosystems, Foster City, CA, USA). Relative expression was calculated using the comparative method 2−ΔΔCt (Livak and Schmittgen, 2001) after normalization with housekeeping genes actin β (ACTB) (ENSBTAG00000026199) and separately GAPDH (ENSBTAG00000014731) for bovine genes (Puech et al., 2015; Horcajo et al., 2017) and BbACT1 and GAPDH1 for the B. besnoiti genes (BbLIS14 in-house assembled genome, unpublished). The primers used to amplify the target genes are listed in Supplementary Tables 1, 2.
Results and Discussion
Sequence Mapping and Quality of RNA-Seq Data
More than 240 million reads were produced in total with 12 different samples. After alignment, an average of 90% of the high-quality reads mapped to the reference B. taurus genome, which is a global indicator of sequencing accuracy. Approximately 0.174 % of the total RNA population was attributable to B. besnoiti, with a higher percentage of mapped reads to B. besnoiti in samples collected at 32 hpi (Supplementary Table 3). To study reproducibility and the experimental variation between replicates, normalized RNA-Seq data were subjected to principal component analysis (PCA), showing that all samples involved in the present work clustered into three biologically distinct groups (Supplementary Figure 1).
Genes expressed commonly associated with ECs, such as thrombospondin, PECAM1 (CD31), Vimentin or Endothelial Specific Molecule (ESM-1), were found, with mean signals above 46.6k normalized counts in all samples confirming the validity and quality of the uninfected control sample.
Host Cell Transcriptome Modulation by B. besnoiti Tachyzoites
Our results showed an early BAEC modulation at 12 hpi characterized by a proinflammatory and procoagulant state upon B. besnoiti invasion and a progression of endothelial activation with upregulation of fibrosis markers along the lytic cycle up to 32 hpi. The gene description and fold change values of a selection of Bos taurus differentially expressed genes (DEGs) at each time point studied are shown in Table 1.
Table 1. Gene description and fold change values of a selection of Bos taurus differentially expressed genes (DEGs) at each time point studied.
Early Parasite Adaptation to Intracellular Lifestyle Seems to Modulate Host TCA Cycle Metabolic Machinery and Down Regulates Mechanisms Involved in Endothelium Protection
When B. besnoiti infected BAEC at 12 hpiwere compared to uninfected BAEC, only 57 DEGs were identified, of which 15 were up-regulated and 42 down-regulated (Supplementary Table 4). This low number of DEGs may be due to the low infection rate achieved with the reference strain employed (BbSpain1 is a low invader and low prolific isolate; Frey et al., 2016), since only up to 20% of the cells were infected (data not shown) biasing the results. Moreover, it has been reported that parasites modulate host cell gene expression only when they have already adapted to the intracellular lifestyle inside the parasitophorous vacuole (Hakimi et al., 2017). Thus, enriched GO terms under the three domains assayed (biological process, molecular function or cellular component) (PANTHER overrepresentation test, Complete GO) were not found. However, two Panther Pathways were significantly enriched: Succinate—proprionate conversion and methylmalonyl pathway. KEGG pathways analysis showed an enrichment of “Carbon metabolism” pathway (bta01200). These results indicate that genes in the tricarboxylic acid cycle (TCA) pathway are regulated in infected BAECThis glucose deprivation from the parasitized host cell might arrest cell cycle progression in conditions of starvation (Kalucka et al., 2015). It is known that the closely related apicomplexan Toxoplasma gondii scavenges carbon sources (glucose and glutamine) from the host cell to fuel its own TCA cycle (Blader and Koshy, 2014; Jacot et al., 2016), and is able to induce a G2/M host cell cycle arrest in primary bovine endothelial cells (BUVEC) (Velásquez et al., 2019). Previous results in B. besnoiti have already shown that the TCA pathway is a key energy pathway for tachyzoites according to Fernández-García et al. (2013), who compared the tachyzoite and bradyzoite stage using proteomic approaches and showed that most of the identified proteins were glycolytic enzymes. Taubert et al. (2016) also showed that the parasite triggered an up-regulation of glycolysis and glutaminolysis in vitro using pathway-specific inhibitors.
Genes involved in the protection of endothelial cells (ECs) were shown to be downregulated at 12 hpi (Table 1, Supplementary Table 4, Figure 1A). Among these, we found genes associated to normal EC function involved in the response to oxidative stress, such as NOX5 and SOD3, protease inhibitors (SERPIN5) and endothelial specific molecules (EPAS1, ECM2). NOX5 is responsible for the generation of reactive oxygen species (ROS) in response to several stimuli and its expression is accompanied by an increase in EC proliferation and promotes their organization into capillary networks (Fulton, 2009). SOD3 is a secreted enzyme responsible for the redox balance in specific tissues, including ECs, preventing oxidative damage and preserving nitric oxide (NO) availability (Fukai and Ushio-Fukai, 2011). SERPIN5 is a glycoprotein that can inhibit serine proteases, including plasminogen activators (Azhar et al., 2013). Our results are thus indicative of endothelial dysfunction in B. besnoiti infected BAEC. Additionally, the highest fold change observed for DEG in infected BAEC at 12 hpi corresponded to ADAMTS1 (Supplementary Table 4), a gene belonging to a family named as a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS). The regulation of ADAMTS1 may suggest an effect on extracellular matrix organization (ECM) as a response to tachyzoite invasion. This family of protease enzymes has been shown to be capable of degrading essentially all components of the ECM to facilitate tissue remodeling (Dancevic et al., 2016). In particular, ADAMTS1 may contribute to the dissemination of Toxoplasma infected leukocytes into immune-privileged sites (Seipel et al., 2010) and its gene expression is also modulated upon N. caninum infection of trophoblast cells (Horcajo et al., 2017).
Figure 1. Heatmaps of a selection of Bos taurus and Besnoitia besnoiti differentially expressed genes. (A) Heatmap showing patterns of expression (normalized reads) in a selection of B. taurus differentially expressed genes implicated in relevant pathways mentioned in the manuscript (e.g., endothelial activation, leukocyte recruitment, proinflammatory phenotypes, extracellular matrix organization). A complete list of B. taurus genes that are included in the figure is available in Supplementary Table 11. (B) Heatmap displaying the transcriptional abundance of a selection of B. besnoiti genes. A complete list of B. besnoiti genes that are included in this figure is available in Supplementary Table 12. The heatmaps were generated using Heatmapper (http://www2.heatmapper.ca/). The genes were clustered using the Pearson computing distance method.
396 genes were DEG when infected BAEC at 32 hpi were compared to uninfected cells, of which 161 were upregulated at 32 hpi. Among those genes upregulated at 32 hpi, the highest fold change observed was for MT1A, and ADAMTS1 (Supplementary Table 5). When this subset of DEG were Analyzed (Supplementary Table 5), terms related to both innate and adaptive immune responses, as well as responses to wounding, were found (Supplementary Table 6, Figure 2A). Besides, KEGG pathways such as “TNF-α signaling pathway (bta 04510)” or “TGF-β signaling pathway (bta04350)” were present (Supplementary Table 6). These results correlate with the activation of ECs described at 12 hpi, but it seems that the proinflammatory response is in a later stage as a consequence of the initial injury of ECs. A remodeling of ECM and angiogenesis regulation are present, evidenced by the expression of genes coding for thrombospondin precursors, fibronectin, ADAMTS, integrin β7 (ITGB7) and vascular endothelial growth factor A (VEGFA), among others.
Figure 2. (A) The scatterplot shows the cluster representative genes of enriched biological process Gene Ontology (GO) terms among differentially expressed Bos taurus genes in Besnoitia besnoiti-infected BAEC cells at 32 hpi vs. non-infected BAEC summarized using REVIGO (Supek et al., 2011). (B) The scatterplot shows the cluster representatives of enriched biological process GO terms among differentially expressed Bos taurus genes in Besnoitia besnoiti-infected BAEC cells at 32 hpi vs. Besnoitia besnoiti-infected BAEC cells at 12 hpi summarized using REVIGO. The axes, semantic spaces x and y, mean that more semantically similar GO terms are closer in the plot.
Type II Endothelial Activation Is Progressively Induced by B. besnoiti Tachyzoites in BAEC Along the Lytic Cycle
Endothelial cell dysfunction was corroborated when infected BAEC at 32 hpi and infected BAEC at 12 hpi were compared, since the highest number of DEG (n = 445) were found, compared to the pair wise comparisons that included infected BAEC at 12 hpi or 32 hpivs. non-infected BAEC (Supplementary Table 7, Figure 1A), together with relevant enriched pathways. Of this subset of DEGs, 249 DEGs were up-regulated and 196 DEGs were down-regulated in infected BAEC at 32 hpi (Supplementary Table 7).
Upon B. besnoiti infection, BAEC showed features compatible with a type II activation (e.g., in response to pathogen stimulation or to cytokines such as IL-1, TNFα and IFNγ) (Mai et al., 2013; Gimbrone and García-Cardeña, 2016) (Figure 2B). Endothelial cell activation comprises a sequence of events involved in many pathological processes: (i) loss of vascular integrity; (ii) upregulation of leukocyte adhesion molecules that leads to leukocyte recruitment, (iii) coagulation including platelet aggregation, coagula formation and reduced fibrinolytic potential (Gimbrone and García-Cardeña, 2016). DEGs related with EC activation were found overexpressed at both pi time points. Genes coding for proteins involved in early activation steps were overexpressed at 12 hpi, whereas DEGs related with a later stage of activation were found overexpressed at 32 hpi.
Loss of vascular integrity was reflected by several genes and GO terms (Supplementary Table 8, Figure 2B) associated with ECM remodeling found to be regulated upon infection. We identified matrix metalloproteinases (such as MMP14) upregulated at 12 hpi, which are able to degrade ECM for tissue remodeling (Parks et al., 2004; Nagase et al., 2006). Previous studies in the closely related parasites, T. gondii and N. caninum, have shown a regulation of other MMPs, which may be important for crossing biological barriers and intra-organic dissemination (Horcajo et al., 2017). Other relevant proteins for ECM remodeling included the upregulation of ADAMTS as previously mentioned. This loss of integrity suggested by these in vitro studies may be responsible for changes in vascular permeability that are thought to be responsible for oedemas and hemorrhages during the acute phase of bovine besnoitiosis (Álvarez-García et al., 2014). Other ECM related genes were also modulated: Integrin alpha 2, integrin alpha V and integrin V10 were upregulated at 12 hpi whilst fibronectin (FN1) collagen-alpha and integrin beta 8 were found to be upregulated at 32 hpi. In fact, integrin-mediated signaling was found to be enriched using Panther tools (Supplementary Table 8). In ECs, integrins mediate the blood flow shear stress response (Muller et al., 1997; Shyy and Chien, 2002). Finally, another DEG involved in the integrity of the endothelial barrier was claudin-1, a gene coding for a protein present in the host cell tight-junctions (Morita et al., 1999), which was also upregulated at 32 hpi. The overexpression of claudin-1 may reflect a wound-healing response initiated after the initial injury and activation of endothelial cells upon parasite invasion, as it has been described for claudin-1 in other systems (Shi et al., 2018).
Another key step in the endothelium activation is the expression of surface markers, which may lead to an increase in adhesion of leukocytes. We found as DEGs canonical markers of endothelial activation such as selections (SELE and SELP), as well as VCAM and ICAM1, upregulated in infected BAEC at 12 hpi when compared to 32 hpi (Supplementary Table 7). E-selectin expression depends on P-selectin upon pro-inflammatory modulators such as IL-1 (Leeuwenberg et al., 1992), which was also upregulated at 12 hpi. VCAM-1 has been shown to present a sustained expression pattern in cultured endothelial cells upon cytokine stimulation (Choi et al., 2018), and shows selective adhesive properties for mononuclear leukocytes and lymphocytes (Osborn et al., 1989). Activated ECs have an intrinsic capacity to secrete chemokines and cytokines such as IL-1, IL-6, CCL2, or CXCL-1, as they form part of the innate immune system and work as a first danger signal sensor (Mai et al., 2013). The secretion of these signaling molecules is responsible for the generation of localized signaling loops which may contribute to the increase in vascular permeability that has been hypothesized to occur in acute bovine besnoitiosis and is induced by other bovine pathogens responsible for endothelial dysfunction such as bluetongue virus (BTV) (DeMaula et al., 2002). The different chemokines secreted by ECs are able to selectively attract different subpopulations of leukocytes, such as monocytes (CCL2 also known as MCP1, Deshmane et al., 2009), neutrophils (CXCL1, Scapini et al., 2004) or eosinophils (CCL24, also known as eotaxin2, Tsai et al., 2016). These findings are supported by TNF-α- KEGG pathway, genes of which were found to be overrepresented in the set of DEGs upon B. besnoiti infection in the 32 vs. 12 hpi comparison, being overexpressed at 12 hpi (Supplementary Table 8, Figure 3). TNF-α is a key proinflammatory cytokine under the control of the NF-kβ transcription factor together with other potent upregulated inflammatory cytokines, such as IL-6 and IL-1A. In healthy barrier cells (such as epithelial and endothelial cells), a steady expression of IL-1A has been demonstrated (Garlanda and Mantovani, 2013; Garlanda et al., 2013), but its expression can be induced by canonical proinflammatory mediators (Weber et al., 2010; Rider et al., 2012). This cytokine is responsible for the establishment of an inflammatory loop, where stressed or damaged cells produce an IL-1A-dependent activation of chemokines that recruit leukocytes. Also, IL-1A is able to modulate the expression of other proinflammatory cytokines, such as IL-1B, at a transcriptional level, which may amplify the inflammatory loop (Di Paolo and Shayakhmetov, 2016). Besides, in vivo studies performed in mice have shown that interactions between TNF-α and IL-6 exacerbate oxidative stress and reduce the phosphorylation of endothelial nitric oxide synthase (eNOS), leading to an endothelial dysfunction (Lee et al., 2017). In this context, the pleiotropic transcription factor NF-kB seems to play a pivotal role since it coordinates the expression of additional effector proteins that appeared as DEGs, such as E-Selectin (SELE), VCAM1, IL-6 and CCL-2, conferring a locally coordinated proinflammatory phenotype. Moreover, several inhibitors and NFkB2 genes were upregulated at 12 hpi indicating that this pathway might be modulated upon parasite-infection.
Figure 3. KEGG pathway for TNFα (bta04668) with the annotated DEGs in our results between Besnoitia besnoiti-infected BAEC cells at 32 hpi vs. Besnoitia besnoiti-infected BAEC cells at 12 hpi (Kanehisa, 2019). Red stars represent DEGs genes in the pathway. Leukocyte recruitment: CCL2, CXCL1, CXCL2, CXCL3; Inflammatory cytokines: IL-6; Intracellular signaling (negative): NFKBIA; IκBa; Remodeling of extracellular matrix: MMP14; Cell adhesion: ICAM1, VCAM1, SELE.
A key gene with an antioxidant effect is metallothionein 1A (MT1A), which was found to be upregulated in infected BAEC at 32 hpi and may reflect an attempt to counteract the inflammatory phenotype induced by the parasite proliferation. Also, genes in the KEGG PPAR peroxisome pathway, which provides vascular protection upon oxidative stress (Mukohda et al., 2015) were found to be modulated upon B. besnoiti infection (Supplementary Table 8).
Our results partially agree with the findings reported by Maksimov et al. (2016), where the transcription of a small panel of adhesion molecules, chemokines, and regulatory molecules (n = 12) was studied upon B. besnoiti infection by qPCR up to 48 hpi. However, in their study a few statistically significant differences were observed, i.e., up regulation of ICAM at 24 hpi, P-SELE at 6 and 12 hpi, chemokines (CXCL-1, CXCL-8, and CCL5) between 6 and 48 hpi and of IL-6 and COX-2 at 48 hpi. In contrast, the relevance of other genes such as VCAM was not proved. These differences in their results compared with our observations might be explained by the different in vitro model employed, including the host cell (BUVEC), B. besnoiti isolate (Bb1Evora04) and the MOI.
Coagulation, exemplified by fibrinolysis pathways, has been also found to be modulated by the parasite in the present work, with the upregulation of several molecules with fibrinolytic properties such as plasminogen activator, urokinase receptor (PLAUR) and plasminogen activator, tissue type (PLAT) at 12 hpi. In addition, fibronectin and thrombospondin are adhesins which favor the platelet adhesion that were also upregulated at 32 hpi.
These pathways induced are in accordance with the microscopic lesions found in vivo, such as vasculitis and thrombosis, in skin and circulatory system (Pols, 1960; McCully et al., 1966; Basson et al., 1970).
BAEC Infected With B. besnoiti Show a Profibrotic Phenotype
Fibrosis plays a crucial role in the pathogenesis of scleroderma, which is characteristic of chronic besnoitiosis. Fibrosis includes the following sequential steps: inflammation, myofibroblast differentiation, ECM accumulation and angiogenesis (Sahin and Wasmuth, 2013). Several fibrosis markers were upregulated at 12 hpi, including chemokines responsible for the recruitment of monocytes and macrophages, such as CCL2, IL-6, and IL-1A and mediators of macrophage differentiation such as macrophage colony stimulating factor 1 (MCSF-1). This macrophage recruitment could explain the predominance of activated macrophages among the leukocyte populations in the inflammatory foci around tissue cysts in chronic besnoitiosis (Frey et al., 2013). The interaction of CCL2 with its receptor CCR2 has been implicated in many fibrotic processes (Sahin and Wasmuth, 2013). Besides its role as a monocyte chemoattractant, this molecule is able to mediate the activation of fibroblasts to produce TGF-β and thus, to stimulate collagen synthesis leading to the deposit of collagen and fibroid tissue (Gharaee-Kermani et al., 1996). Indeed, macrophages are major sources of growth factors that induce ECM deposition. Besides, the KEGG pathway for TGF-β (Supplementary Table 6, Figure 4) was enriched in our subset of DEGs. Since scleroderma is commonly associated with the chronic stage of bovine besnoitiosis, CCL2 may play a key role in the pathogenesis of bovine besnoitiosis. Other relevant DEG were proheparin-binding EGF-like growth factor, which is involved in macrophage mediated cellular proliferation and has mitogenic properties for fibroblasts (Jin et al., 2002), and integrins that allow the anchoring of ECs to components of ECM, such as fibronectin or collagen, and regulating the TGF-β pathway (Henderson et al., 2013).
Figure 4. KEGG pathway for TGF-β (bta04350) with the annotated DEGs in our results between Besnoitia besnoiti-infected BAEC cells at 32 hpi vs. non-infected BAEC cells (Kanehisa, 2019). Red stars represent DEGs in the pathway: Mothers against decapentaplegic homolog 9 (SMAD9); Activin A receptor type 1C (ACVR1C); Inhibin beta E subunit (INHBE); Inhibitor of DNA binding 1 (ID1); ID2; ID3; latent transforming growth factor beta binding protein 1 (LTBP1); thrombospondin 1 (THBS1); v-myc avian myelocytomatosis viral oncogene homolog (MYC).
Angiogenesis was represented by several DEGs. Thus, DEGs involved in both early (growth factors and matrix metalloproteinases) (Brigstock, 2002) and late steps (integrins and vasohibin) evidenced the regulation of this multistep process. The activation of ECs in response to angiogenic factors can lead to a degradation of the endothelial barrier by the secretion of extracellular proteinases (such as MMP or ADAMTS). This degradation lead to a branch point in the vessel wall, with the synthesis of integrins that allow the endothelial cells to migrate toward the angiogenic stimulus. Afterwards, ECs are re-organized to form tubules which are cleaved and interconnected to form a network. Our results show a modulation of several angiogenic molecules (e.g., integrins, growth factors, ADAMTS1, MMP14), including an upregulation of vasohibin-1 at 32 hpi, that is specifically expressed by ECs (Kosaka et al., 2013). Figure 5A shows a graphical abstract of the regulation of host pathways/genes during infection.
Figure 5. Graphical abstract with the main findings of up- and down-regulated host (A) and parasite (B) pathways/genes during an in vitro Besnoitia besnoiti infection in primary bovine endothelial cells. (A) A progressive endothelial dysfunction along the parasite lytic cycle is evidenced. Bos taurus differentially expressed genes are associated with reduced endothelial protection (SOD3, NOX5, ECM2, EPAS1, SERPIN5); a proinflammatory, fibrinolytic and profibrotic phenotype represented by cytokines (IL-1A, IL-6), chemokines (CXCL-2, CXCL-3, CCL-2, CCL24), surface adhesion markers (ICAM-1, VCAM-1, SELE, SELP), genes involved in the coagulation cascade (PLAT, PLAUR) and growth factors (HBEGF, CSF-1; VEGFA) for both post-infection times assayed. (B) Genes coding for surface proteins (SRS), microneme proteins (MIC), AP-2 transcription factors, rhoptry kinases (ROP), rhoptry neck proteins (RON), and genes coding for proteins such as apical membrane antigen-1 (AMA-1) and gliding associated proteins (GAP) were regulated along B. besnoiti lytic cycle. Mo, macrophages; cT, cytotoxic T lymphocytes.
qPCR Validation of Bos Taurus Genes
A total of 17 B. Taurus DEGs (CCL2, CCL24, CXCL2, CXCL3, IL-6, IL-1A, ICAM-1, VCAM-1, SELE, PLAT, PLAUR, NFKB2, thrombospondin1 (THBS1), MT1A, ADAMTS1, FN1, ADAMTS4) were selected for qPCR validation due to their roles in DEG-enriched pathways and/or relevant signaling routes (endothelial activation, leukocyte recruitment, fibrosis and ECM remodeling). Quantitative real-time PCR results showed a similar profile compared with the RNA-Seq results, with a similar significance and direction of fold change in >85% of the validations (Supplementary Figure 2). Even though these genes (SELE, IL-6, IL-1A, CCL2, CCL24, CXCL2, CXCL3) did not appear as DEGs in infected BAEC at 12 hpi vs. non infected BAEC.
Besnoitia besnoiti Transcriptome
Besnoitia besnoiti Tachyzoites Display a Similar Gene Expression Profile of Genes Involved in the Parasite Lytic Cycle at Early and Late Stage Infection
Our results showed that several B. besnoiti genes exclusively present in apicomplexan parasites were expressed at both infection time points, such as dense granule proteins (GRAs) or calcium-dependent kinases (CDPKs), without being differentially expressed when both post-infection times were compared (Figure 1B).
Similarly, in N. caninum, GRAs, such as NcGRA7, NcNTPase (Pastor-Fernández et al., 2016) showed higher expression after early invasion and at egress. This was also described for several T. gondii GRA proteins (Radke et al., 2005) that were shown to be relevant in invasion of the host cell, maturation of the parasitophorous vacuole and egress (Blader et al., 2015).
The present work identified at least eight putative CDPKs in the B. besnoiti genome (Supplementary Table 9) and expressed on the tachyzoite stage. CDPKs are kinases absent in the mammalian host and have therefore garnered significant interest as potential drug targets, allowing the development of specific bumped kinase inhibitors (BKIs), which have been shown to be effective both in vitro and in vivo against a wide range of apicomplexan parasites (reviewed by Van Voorhis et al., 2017) including B. besnoiti (Jiménez-Meléndez et al., 2017). CDPK1 is required for microneme secretion in T. gondii and is thus essential for parasite invasion (Lourido et al., 2010, 2012), while CDPK3 is involved in the K+-mediated egress pathway (McCoy et al., 2012). Besides, BKI treatment of B. besnoiti tachyzoites impaired parasite invasion and proliferation (Jiménez-Meléndez et al., 2017). Also, genes coding for aspartyl proteases (ASP), enzymes important in the cleavage of proteins, have become of increasing interest in the drug design against T. gondii (Dogga et al., 2017) and have been detected in the present work.
However, despite the similar expression profile, we successfully identified 105 DEGs when the two infection time points were compared, and only eight genes were upregulated at 12 hpi compared to 32 hpi (Supplementary Table 10; Figure 1B). We were not able to identify significantly enriched GO terms or pathways due to the high proportion of hypothetical proteins (~40%) in the DEGs.
Among those DEGs, many of which are involved in the lytic cycle in other Toxoplasmatinae parasites, we found genes coding for surface proteins, microneme proteins, AP-2 transcription factors, rhoptry kinases, rhoptry neck proteins, and genes coding for proteins such as apical membrane antigen-1 (AMA-1) and gliding associated proteins (GAP) (Supplementary Table 10), that were regulated along the lytic cycle. Those proteins had not been previously identified by proteomics approaches due to the absence of the whole genome sequence (Fernández-García et al., 2013; García-Lunar et al., 2013, 2014).
SAG-related sequence (SRS) proteins mediate low-affinity parasite interactions with the host-cell membrane for parasite adhesion. Our results have shown that 22 genes encoding SRS proteins were upregulated in tachyzoites at 32 hpi, whilst only one was upregulated at 12 hpi (Supplementary Table 10). One gene encoding a glideosome proteins -GAP80-, which is essential for gliding motility, was upregulated at 32 hpi (Supplementary Table 10).
After the initial recognition of the host cell mediated by SRS proteins, a tighter interaction mediated by microneme proteins leads to the formation of the moving junction (MJ). This MJ is established between the microneme protein AMA1 and proteins from the neck of the rhoptries (RON proteins), which are secreted before kinases from the bulb of the rhoptries (ROPs) (Blader et al., 2015). Our results have shown a differential expression profile for RON3, RON4, RON5, RON6, RON8 and RON10, all of them upregulated at 32 hpi (Supplementary Table 10). This finding may indicate that at 32 hpi, tachyzoites are already preparing the invasion machinery for the entry in new host cells after egress. In T. gondii, it has been shown that RON8 associates to the MJ complex for proper invasion (Straub et al., 2009). On the other hand, RON4 and RON5 have a structural function in the architecture of the complex (Besteiro et al., 2009, 2011). Additionally, RON6 forms a complex with other RON proteins, such as RON2, RON4 and RON8, and is associated with AMA1 (Muniz-Feliciano et al., 2013). RON10 forms a complex with RON9 in T. gondii but it is not detected in the MJ during invasion and it is not associated with the core complex formed by RON2/4/5/8 and AMA1. Finally, it was reported that RON3 was also DE, but in T. gondii it has an unknown function; in P. falciparum, RON3 has been shown to play an important role in merozoite invasion intro erythrocytes (Zhao et al., 2014).
Among the rhoptry bulb proteins, several ROPs have been shown to be DEG, such as ROP5B, ROP17, ROP40, which were overexpressed at 32 hpi (Supplementary Table 10). In T. gondii, several ROPs have been described as effectors that hijack host immune responses (Hakimi et al., 2017). For instance, ROP5, ROP17, and ROP18 forming ROP5/17/18 complex function as virulence factors to prevent the accumulation of immune related GTPases (iRG) at the parasitophorous vacuole (PV) (Saeij et al., 2006). So far, the phosphorylation of iRGs has not been described in B. besnoiti as a strategy to subvert the host-immune responses, but since mice are not a natural host for B. besnoiti and iRGs, as well as TLR11 and TLR12, are absent from the Bos taurus genome this may not be as relevant as described for T. gondii infections in mice. Further studies should address whether these ROPs represent putative virulence factors in B. besnoiti.
Of the genes that encode for inner-membrane complex (IMC) proteins, that play a crucial role in parasite replication, we found 3 DEG: IMC18, IMC22 and IMC24, all of them upregulated at 32 hpi (Supplementary Table 10).
Among transcription factors, the AP2 family have members which coordinates many developmental processes in T. gondii (Walker et al., 2013). However, our results have shown only AP2X6 transcription factor as DEG (Supplementary Table 10), whose function in other coccidian's has not been fully elucidated yet. Figure 5B shows a graphical abstract of the regulation of parasite genes during infection.
qPCR Validation of Besnoitia besnoiti Genes
Eleven B. besnoiti genes were selected for qPCR validation according to their implications in key processes in other apicomplexan parasites (MIC2, MIC11, GAP80, ROP5B, ROP17, ROP40, GRA7, GRA10, SRS22A, PEP, SRS). The qPCR results showed a similar profile to the RNA-Seq results, with a similar significance and direction of fold change in >90% of the validations (Supplementary Figure 3).
Concluding Remarks
This work represents the first transcriptomics approach by means of RNA-Seq to study in vitro host-pathogen interaction between B. besnoiti and target cells during the acute infection. The different steps of ECs activation were progressively induced, a process associated to vascular injury and correlated with microscopic lesions described in tissues from naturally and experimentally infected cattle. Rapid progression of a proinflammatory, procoagulant and profibrotic state was triggered upon parasite invasion.
Additionally, we have identified for the first time many orthologous genes from other Toxoplasmatinae parasites, some of which are upregulated during proliferation (e.g., MIC2, MIC11, ROP40, RON2, GAP80, SRS). This work may help to search for potential drug and vaccine targets and promising markers of disease and prognosis (such as IL-6, IL-1A, CCL2, that showed some of the highest fold change values), that should be further studied in vivo.
Data Availability Statement
The data that support the findings of this study have been deposited in Gene Expression Omnibus (GEO) repository at https://www.ncbi.nlm.nih.gov/geo, with reference number GSE139306.
Author Contributions
GA-G and AH conceived the study and participated in its design. AJ-M performed in vitro infection experiments, functional enrichment and network analysis and wrote the manuscript, with interpretation of results and discussion inputs from CR, GR, AH, and GA-G. GR performed computational analysis of RNA-Seq data. AJ-M designed and performed RT-PCR analyses. All authors read and approved the final manuscript.
Funding
This work was financially supported through a research project from the Spanish Ministry of Economy and Competitiveness (Ref. AGL2016- 75202-R) and Community of Madrid (Ref. P2018/BAA-4370, PLATESA2-CM). AJ-M was supported by a grant from the Spanish Ministry of Education, Culture and Sports (grant n° FPU13/05481).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We would like to acknowledge members from SALUVET research group for their excellent support, specially to Pilar Horcajo for her expert reviewing of the manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcimb.2020.00218/full#supplementary-material
Supplementary Figure 1. 2D mdsPlot of top genes from the samples analyzed in the present work.
Supplementary Figure 2. RT-PCR validation results for Bos taurus differentially expressed genes in every comparison performed. (A) RT-PCR validation results for 12 hpi vs. C-. Positive fold-change values indicate gene upregulation in infected BAEC at 12 hpi. (B) RT-PCR validation results for 32 hpi vs. C-. Positive fold-change values indicate gene upregulation in infected BAEC at 32 hpi. (C) RT-PCR validation results for 32 vs. 12 hpi comparisonPositive fold-change values indicate gene upregulation in infected BAEC at 32 hpi.
Supplementary Figure 3. RT-PCR validation results for Besnoitia besnoiti differentially expressed genes. Positive fold-change values indicate gene upregulation in infected BAEC at 32 hpi.
Supplementary Table 1. Primers used for RT-PCR validation of Bos taurus genes.
Supplementary Table 2. Primers used for RT-PCR validation of Besnoitia besnoiti genes.
Supplementary Table 3. Mapped reads against Bos taurus and Besnoitia besnoiti in all samples analyzed in the present work.
Supplementary Table 4. Differentially expressed genes (DEGs) mapped to Bos taurus at 12 hpi vs. C- comparison. Positive fold-change values indicate gene upregulation in infected BAEC at 12 hpi.
Supplementary Table 5. Differentially expressed genes (DEGs) mapped to Bos taurus at 32hpi vs. C- comparison. Positive fold-change values indicate gene upregulation in infected BAEC at 32 hpi.
Supplementary Table 6. Results from the Gene Ontology (GO) analysis in the category for biological process complete (bp complete) at 32 hpi vs. C- comparison.
Supplementary Table 7. Differentially expressed genes (DEGs) mapped to Bos taurus at 32 hpi vs. 12 hpi comparison. Positive fold-change values indicate gene upregulation in infected BAEC at 32 hpi.
Supplementary Table 8. Results from the Gene Ontology (GO) analysis in the category for biological process complete (bp complete) at 32 vs. 12 hpi.
Supplementary Table 9. Expression profile of selected families of B. besnoiti genes (dense granules, GRA; calcium dependent protein kinases, CDPKs; AP2 transcription factors and Aspartyl proteases, ASP).
Supplementary Table 10. Differentially Expressed Genes (DEGs) mapped to Besnoitia besnoiti at 32 vs. 12 hpi comparison. Positive fold-change values indicate gene upregulation in infected BAEC at 32 hpi.
Supplementary Table 11. Complete list of the abundance of the expression data from B. taurus genes included in Figure 1A.
Supplementary Table 12. Complete list of the abundance of the expression data from B. besnoiti genes included in Figure 1B.
References
Álvarez-García, G., García-Lunar, P., Gutiérrez-Expósito, D., Shkap, V., and Ortega-Mora, L. M. (2014). Dynamics of Besnoitia besnoiti infection in cattle. Parasitology 141, 1419–1435. doi: 10.1017/S0031182014000729
Azhar, A., Singh, P., Rashid, Q., Naseem, A., Sazzad Khan, M., and Aman Jairajpuri, M. (2013). Antiangiogenic function of antithrombin is dependent on its conformational variation: Implication for other serpins. Protein Pept. Lett. 20, 403–411. doi: 10.2174/0929866511320040004
Basson, P., McCully, R., and Bigalke, R. (1970). Observations on the pathogenesis of bovine and antelope strains of Besnoitia besnoiti (Marotel, 1912) infection in cattle and rabbits. Onderstepoort J. Vet. Res. 37, 105–26.
Besteiro, S., Dubremetz, J., and Lebrun, M. (2011). The moving junction of apicomplexan parasites: a key structure for invasion. Cell. Microbiol. 13, 797–805. doi: 10.1111/j.1462-5822.2011.01597.x
Besteiro, S., Michelin, A., Poncet, J., Dubremetz, J. F., and Lebrun, M. (2009). Export of a Toxoplasma gondii rhoptry neck protein complex at the host cell membrane to form the moving junction during invasion. PLoS Pathogens 5:e1000309. doi: 10.1371/journal.ppat.1000309
Bigalke, R., Prozesky, L., Coetzer, J., and Tustin, R. (2004). “Besnoitiosis,” in Infectious Diseases of Livestock, Vol. 1, eds J. A. W. Coetzer and R. C. Tustin (Cape Town: Oxford University Press), 2, 351–359.
Blader, I. J., Coleman, B. I., Chen, C., and Gubbels, M. (2015). Lytic cycle of Toxoplasma gondii: 15 years later. Annu. Rev. Microbiol. 69, 463–485. doi: 10.1146/annurev-micro-091014-104100
Blader, I. J., and Koshy, A. A. (2014). Toxoplasma gondii development of its replicative niche: in its host cell and beyond. Eukaryot. Cell. 13, 965–976. doi: 10.1128/EC.00081-14
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatic 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Brigstock, D. R. (2002). Regulation of angiogenesis and endothelial cell function by connective tissue growth factor (CTGF) and cysteine-rich 61 (CYR61). Angiogenesis 5, 153–165. doi: 10.1023/a:1023823803510
Choi, H., Kim, N., Kim, B., Seo, M., and Heo, J. (2018). TNF-α-induced YAP/TAZ activity mediates leukocyte-endothelial adhesion by regulating VCAM1 expression in endothelial cells. Int. J. Mol. Sci. 19:3428. doi: 10.3390/ijms19113428
Conesa, A., and Gotz, S. (2008). Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int. J. Plant Genomics 2008:619832. doi: 10.1155/2008/619832
Cortes, H., Leitao, A., Gottstein, B., and Hemphill, A. (2014). A review on bovine besnoitiosis: a disease with economic impact in herd health management, caused by Besnoitia besnoiti (Franco and Borges, 1915). Parasitology 141, 1406–1417. doi: 10.1017/S0031182014000262
Dancevic, C. M., McCulloch, D. R., and Ward, A. C. (2016). The ADAMTS hyalectanase family: biological insights from diverse species. Biochem. J. 473, 2011–2022. doi: 10.1042/BCJ20160148
DeMaula, C. D., Leutenegger, C. M., Bonneau, K. R., and MacLachlan, N. J. (2002). The role of endothelial cell-derived inflammatory and vasoactive mediators in the pathogenesis of bluetongue. Virology, 296, 330–337. doi: 10.1006/viro.2002.1476
Deshmane, S. L., Kremlev, S., Amini, S., and Sawaya, B. E. (2009). Monocyte chemoattractant protein-1 (MCP-1): an overview. J. Interferon Cytokine Res. 29, 313–326. doi: 10.1089/jir.2008.0027
Di Paolo, N. C., and Shayakhmetov, D. M. (2016). Interleukin 1α and the inflammatory process. Nat. Immunol. 17:906–913. doi: 10.1038/ni.3503
Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635
Dogga, S. K., Mukherjee, B., Jacot, D., Kockmann, T., Molino, L., Hammoudi, P., et al. (2017). A druggable secretory protein maturase of Toxoplasma essential for invasion and egress. Elife 6:e27480. doi: 10.7554/eLife.27480
European Food Safety Authority (2010). Scientific Statement on Bovine Besnoitiosis. Available online at: http://www.efsa.europa.eu
Fernández-García, A., Álvarez-García, G., Marugán-Hernández, V., García-Lunar, P., Aguado-Martínez, A., Risco-Castillo, V., et al. (2013). Identification of Besnoitia besnoiti proteins that showed differences in abundance between tachyzoite and bradyzoite stages by difference gel electrophoresis. Parasitology 140, 999–1008. doi: 10.1017/S003118201300036X
Frey, C. F., Gutiérrez-Expósito, D., Ortega-Mora, L. M., Benavides, J., Marcen, J. M., Castillo, J. A., et al. (2013). Chronic bovine besnoitiosis: intra-organ parasite distribution, parasite loads and parasite-associated lesions in subclinical cases. Vet. Parasitol. 197, 95–103. doi: 10.1016/j.vetpar.2013.04.023
Frey, C. F., Regidor-Cerrillo, J., Marreros, N., García-Lunar, P., Gutiérrez-Expósito, D., Schares, G., et al. (2016). Besnoitia besnoiti lytic cycle in vitro and differences in invasion and intracellular proliferation among isolates. Parasit Vectors 9:115. doi: 10.1186/s13071-016-1405-9
Fukai, T., and Ushio-Fukai, M. (2011). Superoxide dismutases: role in redox signaling, vascular function, and diseases. Antiox. Redox Signal. 15, 1583–1606. doi: 10.1089/ars.2011.3999
Fulton, D. J. (2009). Nox5 and the regulation of cellular function. Antiox. Redox Signal. 11, 2443–2452. doi: 10.1089/ARS.2009.2587
García-Lunar, P., Regidor-Cerrillo, J., Gutiérrez-Expósito, D., Ortega-Mora, L., and Álvarez-García, G. (2013). First 2-DE approach towards characterising the proteome and immunome of Besnoitia besnoiti in the tachyzoite stage. Vet. Parasitol. 195, 24–34. doi: 10.1016/j.vetpar.2012.12.040
García-Lunar, P., Regidor-Cerrillo, J., Ortega-Mora, L. M., Gutiérrez-Expósito, D., and Álvarez-García, G. (2014). Proteomics reveals differences in protein abundance and highly similar antigenic profiles between Besnoitia besnoiti and Besnoitia tarandi. Vet. Parasitol. 205, 434–443. doi: 10.1016/j.vetpar.2014.09.003
Garlanda, C., Dinarello, C., and Mantovani, A. (2013). The interleukin-1 family: back to the future. Immunity 39, 1003–1018. doi: 10.1016/j.immuni.2013.11.010
Garlanda, C., and Mantovani, A. (2013). Ligands and receptors of the interleukin-1 family in immunity and disease. Front. Immunol. 4:396. doi: 10.3389/fimmu.2013.00396
Gharaee-Kermani, M., Denholm, E. M., and Phan, S. H. (1996). Co-stimulation of fibroblast collagen and transforming growth factor beta1 gene expression by monocyte chemoattractant protein-1 via specific receptors. J. Biol. Chem. 271, 17779–17784.
Gimbrone, M. A. Jr., and García-Cardeña, G. (2016). Endothelial cell dysfunction and the pathobiology of atherosclerosis. Circ. Res. 118, 620–636. doi: 10.1161/CIRCRESAHA.115.306301
Hakimi, M., and Bougdour, A. (2015). Toxoplasma's ways of manipulating the host transcriptome via secreted effectors. Curr. Opin. Microbiol. 26, 24–31. doi: 10.1016/j.mib.2015.04.003
Hakimi, M. A., Olias, P., and Sibley, L. D. (2017). Toxoplasma effectors targeting host signaling and transcription. Clin. Microbiol. Rev. 30, 615–645. doi: 10.1128/CMR.00005-17
Henderson, N. C., Arnold, T. D., Katamura, Y., Giacomini, M. M., Rodriguez, J. D., McCarty, J. H., et al. (2013). Targeting of α v integrin identifies a core molecular pathway that regulates fibrosis in several organs. Nat. Med. 19, 1617–1624. doi: 10.1038/nm.3282
Horcajo, P., Jiménez-Pelayo, L., García-Sánchez, M., Regidor-Cerrillo, J., Collantes-Fernández, E., Rozas, D., et al. (2017). Transcriptome modulation of bovine trophoblast cells in vitro by Neospora caninum. Int. J. Parasitol. 47, 791–799. doi: 10.1016/j.ijpara.2017.08.007
Jacot, D., Waller, R. F., Soldati-Favre, D., MacPherson, D. A., and MacRae, J. I. (2016). Apicomplexan energy metabolism: carbon source promiscuity and the quiescence hyperbole. Trends Parasitol. 32, 56–70. doi: 10.1016/j.pt.2015.09.001
Jiménez-Meléndez, A., Fernández-Álvarez, A., Calle, A., Ramírez, M. A., Diezma-Díaz, C., Vázquez-Arbaizar, P., et al. (2019). Lytic cycle of Besnoitia besnoiti tachyzoites displays similar features in primary bovine endothelial cells and fibroblasts. Parasit. Vectors 12:517. doi: 10.1186/s13071-019-3777-0
Jiménez-Meléndez, A., Ojo, K. K., Wallace, A. M., Smith, T. R., Hemphill, A., Balmer, V., et al. (2017). In vitro efficacy of bumped kinase inhibitors against Besnoitia besnoiti tachyzoites. Int. J. Parasitol. 47, 811–821. doi: 10.1016/j.ijpara.2017.08.005
Jin, K., Mao, X. O., Sun, Y., Xie, L., Jin, L., Nishi, E., et al. (2002). Heparin-binding epidermal growth factor-like growth factor: hypoxia-inducible expression in vitro and stimulation of neurogenesis in vitro and in vivo. J. Neurosci. 22, 5365–5373. doi: 10.1523/JNEUROSCI.22-13-05365.2002
Kalucka, J., Missiaen, R., Georgiadou, M., Schoors, S., Lange, C., De Bock, K., et al. (2015). Metabolic control of the cell cycle. Cell. Cycle 14, 3379–3388. doi: 10.1080/15384101.2015.1090068
Kanehisa, M. (2019). Toward understandind the origin and evolution of cellular organisms. Protein Sci. 28, 1947–1951. doi: 10.1002/pro.3715
Kosaka, T., Miyazaki, Y., Miyajima, A., Mikami, S., Hayashi, Y., Tanaka, N., et al. (2013). The prognostic significance of vasohibin-1 expression in patients with prostate cancer. Br. J. Cancer. 108, 2123–2129. doi: 10.1038/bjc.2013.169
Lawrence, M., Huber, W., Pages, H., Aboyoun, P., Carlson, M., Gentleman, R., et al. (2013). Software for computing and annotating genomic ranges. PLoS Comput. Biol. 9:e1003118. doi: 10.1371/journal.pcbi.1003118
Lee, J., Lee, S., Zhang, H., Hill, M. A., Zhang, C., and Park, Y. (2017). Interaction of IL-6 and TNF-α contributes to endothelial dysfunction in type 2 diabetic mouse hearts. PLoS ONE 12:e0187189. doi: 10.1371/journal.pone.0187189
Leeuwenberg, J. F., Smeets, E. F., Neefjes, J. J., Shaffer, M. A., Cinek, T., Jeunhomme, T. M., et al. (1992). E-selectin and intercellular adhesion molecule-1 are released by activated human endothelial cells in vitro. Immunology 77, 543–549.
Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-delta delta C(T)) method. Methods. 25, 402–408. doi: 10.1006/meth.2001.1262
Lourido, S., Shuman, J., Zhang, C., Shokat, K. M., Hui, R., and Sibley, L. D. (2010). Calcium-dependent protein kinase 1 is an essential regulator of exocytosis in Toxoplasma. Nature 465, 359–362. doi: 10.1038/nature09022
Lourido, S., Tang, K., and Sibley, L. D. (2012). Distinct signalling pathways control Toxoplasma egress and host-cell invasion. EMBO J. 31, 4524–4534. doi: 10.1038/emboj.2012.299
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8
Mai, J., Virtue, A., Shen, J., Wang, H., and Yang, X. (2013). An evolving new paradigm: endothelial cells–conditional innate immune cells. J. Hematol. Oncol. 6:61. doi: 10.1186/1756-8722-6-61
Maksimov, P., Hermosilla, C., Kleinertz, S., Hirzmann, J., and Taubert, A. (2016). Besnoitia besnoiti infections activate primary bovine endothelial cells and promote PMN adhesion and NET formation under physiological flow condition. Parasitol. Res. 115, 1991–2001. doi: 10.1007/s00436-016-4941-5
McCoy, J. M., Whitehead, L., van Dooren, G. G., and Tonkin, C. J. (2012). TgCDPK3 regulates calcium-dependent egress of Toxoplasma gondii from host cells. PLoS Pathogens 8:e1003066. doi: 10.1371/journal.ppat.1003066
McCully, R. M., Basson, P. A., Van Niekerk, J. W., and Bigalke, R. D. (1966). Observations on Besnoitia cysts in the cardiovascular system of some wild antelopes and domestic cattle. Onderstepoort J. Vet. Res. 33, 245–276.
Morita, K., Sasaki, H., Furuse, M., and Tsukita, S. (1999). Endothelial claudin: claudin-5/TMVCF constitutes tight junction strands in endothelial cells. J. Cell. Biol. 147, 185–194.
Mukohda, M., Stump, M., Ketsawatsomkron, P., Hu, C., Quelle, F. W., and Sigmund, C. D. (2015). Endothelial PPAR-γ provides vascular protection from IL-1β-induced oxidative stress. Am. J. Physiol. Heart. Circ. Physiol. 310, H39–H48. doi: 10.1152/ajpheart.00490.2015
Muller, J. M., Chilian, W. M., and Davis, M. J. (1997). Integrin signaling transduces shear stress-dependent vasodilation of coronary arterioles. Circ. Res. 80, 320–326.
Muniz-Feliciano, L., Van Grol, J., Portillo, J. C., Liew, L., Liu, B., Carlin, C. R., et al. (2013). Toxoplasma gondii-induced activation of EGFR prevents autophagy protein-mediated killing of the parasite. PLoS Pathogens 9:e1003809. doi: 10.1371/journal.ppat.1003809
Nagase, H., Visse, R., and Murphy, G. (2006). Structure and function of matrix metalloproteinases and TIMPs. Cardiovasc. Res. 69, 562–573. doi: 10.1016/j.cardiores.2005.12.002
Neill, J. D., Ridpath, J. F., Lange, A., and Zuerner, R. L. (2008). Bovine viral diarrhoea virus infection alters global transcription profiles in bovine endothelial cells. Dev. Biol (Basel). 132, 93–98. doi: 10.1159/000317148
Osborn, L., Hession, C., Tizard, R., Vassallo, C., Luhowskyj, S., Chi-Rosso, G., et al. (1989). Direct expression cloning of vascular cell adhesion molecule 1, a cytokine-induced endothelial protein that binds to lymphocytes. Cell 59, 1203–1211.
Parks, W. C., Wilson, C. L., and Lopez-Boado, Y. S. (2004). Matrix metalloproteinases as modulators of inflammation and innate immunity. Nat. Rev. Immunol. 4, 617–629. doi: 10.1038/nri1418
Pastor-Fernández, I., Regidor-Cerrillo, J., Álvarez-García, G., Marugán-Hernández, V., García-Lunar, P., Hemphill, A., et al. (2016). The tandemly repeated NTPase (NTPDase) from Neospora caninum is a canonical dense granule protein whose RNA expression, protein secretion and phosphorylation coincides with the tachyzoite egress. Parasit. Vectors 9:352. doi: 10.1186/s13071-016-1620-4
Pols, J. W. (1960). Studies on bovine besnoitiosis with special reference to the aetiology. Onderstepoort J. Vet. Res. 28, 265–356.
Puech, C., Dedieu, L., Chantal, I., and Rodrigues, V. (2015). Design and evaluation of a unique SYBR green real-time RT-PCR assay for quantification of five major cytokines in cattle, sheep and goats. BMC Vet. Res. 11:65. doi: 10.1186/s12917-015-0382-0
Radke, J. R., Behnke, M. S., Mackey, A. J., Radke, J. B., Roos, D. S., and White, M. W. (2005). The transcriptome of Toxoplasma gondii. BMC Biol. 3:26. doi: 10.1186/1741-7007-3-26
Reid, A. J., Vermont, S. J., Cotton, J. A., Harris, D., Hill-Cawthorne, G. A., Konen-Waisman, S., et al. (2012). Comparative genomics of the apicomplexan parasites Toxoplasma gondii and Neospora caninum: coccidia differing in host range and transmission strategy. PLoS Pathogens 8:e1002567. doi: 10.1371/journal.ppat.1002567
Rider, P., Kaplanov, I., Romzova, M., Bernardis, L., Braiman, A., Voronov, E., et al. (2012). The transcription of the alarmin cytokine interleukin-1 alpha is controlled by hypoxia inducible factors 1 and 2 alpha in hypoxic cells. Front. Immunol. 3:290. doi: 10.3389/fimmu.2012.00290
Saeij, J. P., Boyle, J. P., Coller, S., Taylor, S., Sibley, L. D., Brooke-Powell, E. T., et al. (2006). Polymorphic secreted kinases are key virulence factors in toxoplasmosis. Science (New York, NY). 314, 1780–1783. doi: 10.1126/science.1133690
Sahin, H., and Wasmuth, H. E. (2013). Chemokines in tissue fibrosis. Biochim. Biophys. Acta 1832, 1041–1048. doi: 10.1016/j.bbadis.2012.11.004
Scapini, P., Morini, M., Tecchio, C., Minghelli, S., Di Carlo, E., Tanghetti, E., et al. (2004). CXCL1/macrophage inflammatory protein-2-induced angiogenesis in vivo is mediated by neutrophil-derived vascular endothelial growth factor-A. J. Immunol. 172, 5034–5040. doi: 10.4049/jimmunol.172.8.5034
Schares, G., Venepally, P., and Lorenzi, H. A. (2017). Draft genome sequence and annotation of the apicomplexan parasite Besnoitia besnoiti. Genome Announc. 5, e01200–17. doi: 10.1128/genomeA.01200-17
Seipel, D., Oliveira, B. C., Resende, T. L., Schuindt, S. H., Pimentel, P. M., Kanashiro, M. M., et al. (2010). Toxoplasma gondii infection positively modulates the macrophages migratory molecular complex by increasing matrix metalloproteinases, CD44 and alpha v beta 3 integrin. Vet. Parasitol, 169, 312–319. doi: 10.1016/j.vetpar.2009
Shi, J., Barakat, M., Chen, D., and Chen, L. (2018). Bicellular tight junctions and wound healing. Int. J. Mol. Sci. 19:E3862. doi: 10.3390/ijms19123862
Shyy, J. Y., and Chien, S. (2002). Role of integrins in endothelial mechanosensing of shear stress. Circ. Res. 91, 769–775. doi: 10.1161/01.res.0000038487.19924.18
Straub, K. W., Cheng, S. J., Sohn, C. S., and Bradley, P. J. (2009). Novel components of the apicomplexan moving junction reveal conserved and coccidia-restricted elements. Cell. Microbiol. 11, 590–603. doi: 10.1111/j.1462-5822.2008.01276.x
Supek, F., Bosnjak, M., Skunca, N., and Smuc, T. (2011). REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE 6:e21800. doi: 10.1371/journal.pone.0021800
Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–52. doi: 10.1093/nar/gku1003
Taubert, A., Hermosilla, C., Silva, L., Wieck, A., Failing, K., and Mazurek, S. (2016). Metabolic signatures of Besnoitia besnoiti-infected endothelial host cells and blockage of key metabolic pathways indicate high glycolytic and glutaminolytic needs of the parasite. Parasitol. Res. 115, 2023–2034. doi: 10.1007/s00436-016-4946-0
Tsai, C. S., Huang, C. Y., Chen, C. H., Lin, Y. W., Shih, C. M., Tsao, N. W., et al. (2016). Eotaxin-2 increased toll-like receptor 4 expression in endothelial cells in vitro and exacerbates high-cholesterol diet-induced atherogenesis in vivo. Am. J. Transl. Res. 8, 5338–5353.
Van Voorhis, W. C., Doggett, J. S., Parsons, M., Hulverson, M. A., Choi, R., Arnold, S., et al. (2017). Extended-spectrum antiprotozoal bumped kinase inhibitors: a review. Exp. Parasitol. 180, 71–83. doi: 10.1016/j.exppara.2017.01.001
Velásquez, Z. D., Conejeros, I., Larrazabal, C., Kerner, K., Hermosilla, C., and Taubert, A. (2019). Toxoplasma gondii-induced host cellular cell cycle dysregulation is linked to chromosome missegregation and cytokinesis failure in primary endothelial host cells. Sci. Rep. 9:12496. doi: 10.1038/s41598-019-48961-0
Walker, R., Gissot, M., Huot, L., Alayi, T. D., Hot, D., Marot, G., et al. (2013). Toxoplasma transcription factor TgAP2XI-5 regulates the expression of genes involved in parasite virulence and host invasion. J. Biol. Chem. 288, 31127–31138. doi: 10.1074/jbc.M113.486589
Weber, A., Wasiliew, P., and Kracht, M. (2010). Interleukin-1 (IL-1) pathway. Sci. Signal. 3:cm1. doi: 10.1126/scisignal.3105cm1
Keywords: Besnoitia besnoiti, tachyzoite, primary bovine aorta endothelial cells (BAEC), transcriptome, RNA-Seq
Citation: Jiménez-Meléndez A, Ramakrishnan C, Hehl AB, Russo G and Álvarez-García G (2020) RNA-Seq Analyses Reveal That Endothelial Activation and Fibrosis Are Induced Early and Progressively by Besnoitia besnoiti Host Cell Invasion and Proliferation. Front. Cell. Infect. Microbiol. 10:218. doi: 10.3389/fcimb.2020.00218
Received: 09 January 2020; Accepted: 20 April 2020;
Published: 15 May 2020.
Edited by:
Tiago W. P. Mineo, Federal University of Uberlandia, BrazilReviewed by:
Carsten Lüder, Universitätsmedizin Göttingen, GermanyVictoria Jeffers, University of New Hampshire, United States
Copyright © 2020 Jiménez-Meléndez, Ramakrishnan, Hehl, Russo and Álvarez-García. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Gema Álvarez-García, gemaga@ucm.es