RNA-Seq Analyses Reveal That Endothelial Activation and Fibrosis Are Induced Early and Progressively by Besnoitia besnoiti Host Cell Invasion and Proliferation

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 . Disease initiates with the acute stage, when the tachyzoites are fastreplicating 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 , 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 postinfection (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.

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 × 10 6 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 R (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 R homogenization according to the manufacturer's instructions. Briefly, samples in RNAlater R were lysed in 900 µL of RLT buffer + β-mercaptoethanol at 10 µl/mL and column purified; an oncolumn 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 R (1.0) Fluorometer and the Caliper GX LabChip R 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://stringdb.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 R VILO TM 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 R 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 highquality 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.

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 (  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). 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  C-X-C motif chemokine ligand 2. It is also known as Growth-regulated protein homolog beta (GROβ). It has chemotactic activity for neutrophils acting after CXCL-1, also implicated in the branching of endothelial cells.

CXCL3
C-X-C motic chemokine ligand 3. It is also known as Growth-regulated protein homolog gamma (GROγ). It has chemotactic activity for neutrophils. It may play a role in inflammation and exert its effects on endothelial cells in an autocrine manner.  NFKBIZ NF-kappa-B inhibitor zeta (IKBZ). It is involved in NF-kappa-B transcription factor complexes regulation. It inhibits NFκβ activity without affecting its nuclear translocation upon stimulation, but it has been also shown to present a transcriptional activation activity, being involved in the induction of inflammatory genes activated through TLR/IL-1 receptor signaling.
*Whilst for some of the comparisons the fold change value would suggest differential expression, those genes were not considered as DE since they failed to meet the fdr requisite (<0.05). For a list of all genes considered as DE in every comparison, see Supplementary Tables 4, 5, 7. 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.   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 . Other ECM related genes were also modulated: Integrin alpha 2, integrin alpha V and integrin V10 were upregulated at 12 hpi whilst fibronectin (FN1) collagenalpha 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 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.
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 , 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.
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 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). 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).
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.

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(Lourido et al., , 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. 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(Besteiro et al., , 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.