Trypanosoma cruzi dysregulates expression profile of piRNAs in primary human cardiac fibroblasts during early infection phase

Trypanosoma cruzi, the etiological agent of Chagas Disease, causes severe morbidity, mortality, and economic burden worldwide. Though originally endemic to Central and South America, globalization has led to increased parasite presence in most industrialized countries. About 40% of infected individuals will develop cardiovascular, neurological, and/or gastrointestinal pathologies. Accumulating evidence suggests that the parasite induces alterations in host gene expression profiles in order to facilitate infection and pathogenesis. The role of regulatory gene expression machinery during T. cruzi infection, particularly small noncoding RNAs, has yet to be elucidated. In this study, we aim to evaluate dysregulation of a class of sncRNAs called piRNAs during early phase of T. cruzi infection in primary human cardiac fibroblasts by RNA-Seq. We subsequently performed in silico analysis to predict piRNA-mRNA interactions. We validated the expression of these selected piRNAs and their targets during early parasite infection phase by stem loop qPCR and qPCR, respectively. We found about 26,496,863 clean reads (92.72%) which mapped to the human reference genome. During parasite challenge, 441 unique piRNAs were differentially expressed. Of these differentially expressed piRNAs, 29 were known and 412 were novel. In silico analysis showed several of these piRNAs were computationally predicted to target and potentially regulate expression of genes including SMAD2, EGR1, ICAM1, CX3CL1, and CXCR2, which have been implicated in parasite infection, pathogenesis, and various cardiomyopathies. Further evaluation of the function of these individual piRNAs in gene regulation and expression will enhance our understanding of early molecular mechanisms contributing to infection and pathogenesis. Our findings here suggest that piRNAs play important roles in infectious disease pathogenesis and can serve as potential biomarkers and therapeutic targets.


Introduction
Trypanosoma cruzi is an intracellular protozoan parasite that causes Chagas disease. This disease, which was originally endemic in Central and South America, is now present in all industrially advanced countries due modern globalization (Coura and Vinas, 2010;Gascon et al., 2010;Hotez et al., 2013;Bern et al., 2019;Bonney et al., 2019), qualifying the disease as an emerging global health problem (Tanowitz et al., 2009). Autochthonous T. cruzi transmissions including vertical transmission from mother to the fetus have been reported in both inland states and those sharing a border with Mexico (Dorn et al., 2007;Rassi et al., 2010;Beatty and Klotz, 2020;Lynn et al., 2020). Although the acute phase of the disease is associated with indeterminate symptoms and high parasitemia, it remains an essential prelude to the chronic condition where about 30-40% of infected individuals present with incurable cardiac, neurological, or gastrointestinal tract pathological conditions (Marin-Neto et al., 2007;Matsuda et al., 2009;Pinazo et al., 2010;Py, 2011;Martinez et al., 2019;Ricci et al., 2020). Chagas disease is the world's leading cause of infectious myocarditis and the molecular basis of T. cruzi-induced Chagasic cardiomyopathy remains unknown despite ongoing research (Bonney and Engman, 2008;Rassi et al., 2009). Our laboratory among others have used several in vitro, ex vivo, human and murine cardiomyocyte culture models to show that the parasite induces an increase in the expression of extracellular matrix components during the acute phase of infection, which could lead to the cardiac remodeling observed in some chronically infected individuals (Cruz et al., 2016;Udoko et al., 2016;da Costa et al., 2019). Fibroblasts constitute an essential component of the heart where they play important roles in structure and function of the heart under normal homeostatic conditions (Souders et al., 2009;Ivey and Tallquist, 2016). Fibroblasts also play a role in the regulation of the extracellular and matricellular contents of the heart (Herum et al., 2017;Hortells et al., 2019;Silva et al., 2020). Elucidation of how the parasite dysregulates host gene expression profiles in cardiac myocytes and fibroblasts will advance our understanding of T. cruzi induced cardiac pathogenesis.
Accumulating literature suggests that host small noncoding RNA (sncRNA) molecules such as microRNAs (miRNAs), small interfering RNAs (siRNAs), and P-element induced wimpy testis (PIWI)-interacting RNAs (piRNAs) play important roles in regulating gene expression by forming complexes with Argonaute proteins to recognize specific target sequences (Shimoni et al., 2007;Dana et al., 2017;O'Brien et al., 2018;Rayford et al., 2021). To date, several thousand piRNAs have been identified in somatic and germ stem cells where they are thought to play important epigenetic roles including, but not limited to transposon silencing and de novo methylation to maintain genome integrity (Peng and Lin, 2013;Rayford et al., 2021). Dysregulated expression of piRNAs has been suggested to play important roles in tumorigenesis in a context dependent manner (Rizzo et al., 2016;AboEllail and Hata, 2017;Weng et al., 2018;Liu et al., 2019;Halajzadeh et al., 2020). In cardiac cells, it was suggested that piRNAs interact with PIWIL2 protein to facilitate cardiomyocyte proliferation and regeneration by regulating AKT signaling (Rajan et al., 2014). Patients presenting with myocardial infarction were shown to have a significantly elevated level of piR-2106027, which has been suggested to be an important diagnostic marker for the disease (Xuan et al., 2017;Das et al., 2018). Recently, we showed in infectious disease research that T. cruzi dysregulates the expression profile of known and novel piRNAs in primary human cardiac myocytes during the early phase of infection (Rayford et al., 2020). Our in silico analysis showed that significantly differentially expressed piRNAs could target regions in genes coding for NFATC2, FOS, and TGFB1, which are reported to be critical during the early phase of T. cruzi pathogenesis (Rayford et al., 2020). Despite these important functions of piRNAs, the contextual characteristics and molecular signature of the T. cruziinduced piRNA profile in primary human fibroblasts (PHCF) during the early phase of infection remains unknown. Here, we challenged PHCF with the Tulahuen strain of T. cruzi, clone MMC20A and evaluated the piRNA expression profile (piRNome kinetics) during the early phase of infection. We validated the piRNA expression and target mRNA transcript expression. The significantly differentially expressed piRNAs were computationally predicted to target genes including ICAM1, SMAD2, EGR1, CXCR2, FOS, and CX3CL1. Furthermore, we connected them in biological interaction networks to theoretically predict piRNA pathway-level interactions involved in T. cruzi pathogenesis.
2 Materials and methods 2.1 Primary human cardiac fibroblast culture PHCF were obtained from and cultured following the manufacturer's recommendations (PromoCell, Heidelberg, Germany). Briefly, the PHCF were cultured in fibroblast basal growth medium supplemented with the supplemental mix (PromoCell, Heidelberg, Germany) containing fetal calf serum (0.05 mL/mL), recombinant human epidermal growth factor (0.5 ng/mL), recombinant human basic fibroblast growth factor (2 ng/ mL) and recombinant human insulin (5 ug/mL). The cells were cultured in T75 flasks at 37°C in the presence of 5% CO 2 to approximately 80% confluency (approximately 4 × 10 6 cells) prior to being used in our assays.

Parasite culture and infection assays
Heart myoblast monolayers at 80% confluence, cultured in complete DMEM containing 5% glutamax, 10% fetal bovine serum, 1% each of penicillin/streptomycin, multivitamins and MEM nonessential amino acids (Life Technologies, Carlsbad, CA, USA), were infected with T. cruzi trypomastigotes. Pure cultures of highly invasive T. cruzi trypomastigotes (clone MMC 20A, Tulahuen strain) were harvested from the supernatant of infected heart myoblast monolayers as previously described (Lima and Villalta, 1989;Villalta et al., 1990). The parasites were washed with Hanks Balanced Salt Solution (HBSS) and resuspended in PHCF growth medium without supplement at 1 × 10 7 parasites/mL. For the infection assays, confluent PHCF monolayers (about 80%) were starved in HBSS containing 30 mM HEPES, followed by the addition of T. cruzi trypomastigotes in PHCF growth medium without supplements at a ratio of 10 parasites per cell. Parasitechallenged PHCFs were incubated for 1, 3, and 6 h, respectively, in triplicate for the assays. Total and small RNAs were purified from the samples. Mock-infected (basal media only) PHCFs served as control.

RNA extraction and quality assessment
Control and parasite-challenged PHCF were washed with HBSS. The cells were lysed in QIAzol, an RNA extraction lysis buffer and extracted with chloroform following the manufacturer's instructions (Qiagen, Valencia, CA, USA). The aqueous phase of the extract was mixed with an equal volume of 70% ethanol and passed through the RNeasy Mini spin column. The eluate, which contained small RNA species, was mixed with 0.65 volumes of pure ethanol and passed through an RNeasy MiniElute spin column. The column was washed with 80% ethanol and bound small RNA species were eluted with RNase-free water essentially as described by the manufacturer (Qiagen). Large RNAs bound to the RNeasy mini spin column were washed and eluted with RNase-free water as described by the manufacturer (Qiagen). The quality of the purified RNA was analyzed using the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) to determine the RNA Integrity Number (RIN). Samples with a RIN of at least 8 were considered for further analysis.
2.4 RNA-sequencing of sncRNA, filtering, and expression evaluation Briefly, purified small RNA species ranging from 18 to 30 nt were ligated to the Illumina 3′ adaptor and 5′ adaptor. Ligation products were gel-purified, reverse-transcribed and amplified by Rolling-Circle Replication (RCR). This linear amplification only copies the original DNA template instead of copy-of-a-copy in order to make small, very high-density sequencing compacted templates called DNA NanoBalls (DNBs). The DNBs were compacted on high-density patterned nanoarray and sequenced by combinatorial Probe-Anchor Synthesis (cPAS), a sequencing chemistry technique which is optimized for DNBseq. The combination of linear amplification and DNB technology reduces the error rate while enhancing the signal. Smaller DNB spots with highly concentrated DNA deliver significantly higher dye density than PCR cluster arrays, leading to higher signal-to-noise for optimal imaging signal integrity. High-throughput sequencing (HTS) was conducted using the BGISEQ-500 sequencing platform for small RNA (BGI, Cambridge, MA, USA) utilizing the NGS 2.0 DNBseq technology. The DNBseq sequencing technology (NGS 2.0) combines the power of DNA Nanoballs (DNB), PCR-free Rolling Circle Replication, Patterned Nano Arrays and cPAS to deliver a new level of data clarity and reliability. The BGISEQ-500 sequencing platform has comparable sensitivity and accuracy in terms of quantification of gene expression, and low technical variability as compared to the Illumina HiSeq platform (Fehlmann et al., 2016;Natarajan et al., 2019). High-throughput sequencing (HTS) was conducted using the BGISEQ-500 sequencing platform for small RNA (BGI, Cambridge, MA, USA) (Drmanac et al., 2010;Fu et al., 2018). In total, 29,496,864 raw HTS data were filtered by eliminating low-quality reads (less than 20 in Phred quality score), as well as removing confounders such as adaptors and other contaminants including missing 3′ primers, 5′ primer contamination and incomplete small non-coding RNA (sncRNA) reads of less than 18 bp. Reads that met our analysis criteria were subjected to size filtration for piRNA to select for sncRNA transcripts that were 25-30 bp in length and did not match any known miRNA or siRNA sequences. Bowtie was used to map the reads to reference genomes (Langmead et al., 2009). The piRNA annotation program (Piano) was used to predict known piRNAs via a support vector machine (SVM) algorithm with transposon interaction informatics (Wang et al., 2014). piRNA expression was determined using the standard transcripts per kilobase million mapped (TPM) method (John et al., 2004).

Analysis of differential piRNA expression
The NOISeq method (version 3.34) was used to determine differentially expressed piRNAs. Each time point was screened for significantly different piRNA expression compared to control by calculating the log 2 fold-change (M) and absolute differential value (D) between each pair of time points to build a noise distribution model. A piRNA is counted as differentially expressed if M and D values are likely to be higher than in noise. The significance threshold for differential piRNA expression was set to fold-change ≥2 or ≤2 and a divergence probability ≥0.8. A probability of 0.8 is equivalent to an odds value of 4:1, meaning that a given piRNA is 4 times more likely to be differentially expressed than non-differentially expressed (Tarazona et al., 2011). NOISeq calculates probability of differential expression (qvalue) rather than FDR for multiple comparisons correction. This statistic is considered comparable to the Benjamini and Hochberg FDR as explained in Efron et al. (Efron et al., 2001).

piRNA target prediction
Differentially expressed piRNAs were compared against hg38 RefSeq transcripts in miRanda (Enright et al., 2003) using a highstringency pairing score cutoff of ≥175, an energy cutoff of ≤−30 kcal/mol and a requirement for exact seed region alignment.

Quantitation of mRNA and piRNA expression
Total RNA from T. cruzi challenged PHCF was reverse transcribed into cDNA using iScript cDNA synthesis Kit (Bio-Rad, Hercules, CA, USA). Quantitative real-time PCR (qPCR) for validation of the mRNA targets was done using a customized PrimePCR assay (Bio-Rad, Hercules, CA, USA) containing the genes of interest. In the PrimePCR array experiment, total RNA (1 µg) converted to cDNA using the iScript cDNA synthesis kit was mixed with SsoAdvanced SYBR green 2X master mix and loaded (20 µl/per well) on a customized PrimePCR plate containing primers for selected genes to be validated, including housekeeping genes and controls as described (www.bio-rad.com/PrimePCR). The PCR amplification was carried out on a C1000 Touch Thermal Cycler as described by the manufacturer (Bio-Rad, Hercules, CA, USA). The Ct values of housekeeping genes (beta-2 microglobulin, Ribosomal protein, lateral stalk subunit (RPLO) and, beta-actin) were checked for consistency on all plates across all samples. The data were normalized against the housekeeping genes and control samples using CFX manager analysis software with support from Bio-Rad technical service using the DDCT method. The primer sequences used for quantification of the transcript levels in the PrimePCR were not disclosed by Bio-Rad.
For piRNA quantification, we used the TaqMan Small RNA Assay essentially as described by the manufacturer (Thermofisher, Waltham, MA, USA). Briefly, piRNA-specific stem-loop RT primers were used to create each cDNA template. The TaqMan Small RNA Assay used (assay ID) were as follows: novel_piR_17 (CTNKRW3), hsa_piR_016828 (CTH497C); hsa_piR_017716 (CTGZFME); hsa_piR_016742 (CTFVK2G); novel_piR_167 (CTPRKGZ). For each piRNA, a TaqMan primer-probe set was utilized to amplify each target cDNA by TaqMan qPCR, using U6 as a housekeeping sncRNA. TaqMan real-time PCR assays were carried out on a CFX96 Thermal Cycler (Bio-Rad) as per Thermofisher recommendations. The Ct values of housekeeping gene U6 were checked for consistency on all plates across all samples. The data were normalized against the housekeeping U6 and control samples using CFX manager analysis software using the DDCT method. The primer sequences used to quantify the piRNA levels in the TaqMan Small RNA were not disclosed by Thermofisher.

Target gene selection mapping of biological pathway interactions
To identify piRNA target genes to evaluate with an unbiased approach, we used the Genshot search engine (Lachmann et al., 2019) (https://maayanlab.cloud/geneshot/). Geneshot produces a ranked list of genes most relevant to user-provided search terms. These lists integrate information from publications and gene expression data to identify both common and understudied genes. We used the terms "cruzi" and "fibrosis" to create a prioritized list which was matched with the piRNA target list. The resulting genes were used as input for network construction. Biological interaction network construction was conducted with the GeneMANIA algorithm (Warde-Farley et al., 2010) by querying multiple biological interaction databases including GEO, BioGRID and EMBL-EBI. Predicted piRNA target genes were set as starting seed nodes, and the network was then expanded to one degree of biological interaction using GeneMANIA and visualized with Gephi (Mathieu Bastian and Jacomy., 2009) (https://gephi.org/).

Data availability
All relevant data not presented in the manuscript are located in the SRA database: SRA accession: SRX12127485 3 Results 3.1 piRNAs are differentially expressed in PHCF during early phase of T. cruzi infection To evaluate if T. cruzi can dysregulate piRNA expression in heart cells, we challenged PHCF with invasive T. cruzi trypomastigotes Tulahuen strain clone MMC20A. Small RNAs were purified and sent to BGI (Cambridge, MA, USA) for sequencing, as previously described (Rayford et al., 2020). To evaluate the clustering and variance between biological replicates, we performed principal component analysis. The analytical output showed segregation of 0, 1, 3, and 6 hour time points in the PCA scatter plot. Hierarchical clustering and heatmap show T. cruzi challenged samples clustered independently of the control group. Principal component 1 (33.0%) and principal component 2 (17.4%) showed clustering of each replicate group ( Figure 1A). Furthermore, to evaluate whether the gene expression profiles of sncRNAs at each time point was specific, we generated hierarchical clustering, which showed that biological replicates of each time point clustered together ( Figure 1B).
Others suggested that piRNAs can be generated from transposable element (TE) sequences to target and repress TEs, therefore, piRNAs can be classified into TE families (Weick and Miska, 2014;Nandi et al., 2016;Wang and Lin, 2021). We evaluated the origin of the differentially expressed (DE) piRNAs and mapped them to their various genomic loci such as DNA, long interspersed nuclear elements (LINE), and long terminal repeat (LTR) TE families; we found that during the early phase of T. cruzi infection of PHCF, the majority of DE piRNAs were derived from LTRs, representing 53% at 1 h, 43% at 3 h, and 49% at 6 h, respectively (Figures 1C-E).
3.2 T. cruzi induces differential expression of known and novel piRNAs in PHCF during the early phase of infection To evaluate significant changes in the piRNA expression profile of PHCF during parasite infection, we used NOISeq to analyze normalized DE (Tarazona et al., 2015). Violin plot ( Figure 2A) and volcano plots (Figures 2B-D) were used to visualize and identify DE across the various time points during parasite challenge. We found that 441 piRNAs were significantly DE during our experimental parasite challenge. Of these, 235 piRNAs were upregulated throughout the course of infection, while 3 were downregulated when compared to control ( Figure 2E). Supplemental Tables 1, 2 show statistical significance for each DE piRNA and sequence, respectively. Of the piRNAs we found to be dysregulated, a total of 29 known piRNAs were upregulated. Fifteen (51.7%), were upregulated at all the time points ( Figure 2F). The majority of the dysregulated piRNAs during early T. cruzi infection of PHCF were putative and novel; 220 were upregulated and 3 were downregulated at all time points relative to control ( Figure 2G). Putative piRNAs accounted for all the downregulated piRNAs. The greatest count of DE putative piRNAs was 52 at 1h compared to 28 at 3 h, and 25 at 6 h, respectively.

Parasite induced DE piRNAs were computationally predicted to target fibrogenic and inflammatory genes
Accumulating evidence from our laboratory and others suggest that T. cruzi alters host cell expression profiles during the course of infection. In particular, we and others showed that the parasite dysregulates the expression profile of profibrotic and inflammatory molecules at both the transcript and protein levels through unknown mechanisms (Udoko et al., 2016;Suman et al., 2018).
To determine if any of the DE piRNAs could potentially target any of the reported fibrotic and inflammatory genes, we used the miRANDA algorithm (Enright et al., 2003) to computationally predict potential target binding sites. We then used Geneshot (Lachmann et al., 2019) to obtain a ranked list of genes correlating with search terms "cruzi" and "fibrosis." We used the overlap of the prioritized gene list and piRNA targets to identify genes for further analysis including ICAM1, SMAD2, EGR1, CX3CL1, and CXCR2.
Interestingly, chemokine signaling has been shown to contribute to fibrotic remodeling via fibrogenic/alternative macrophages and through pro-fibrotic actions such as recruitment of activated neutrophils and fibrogenic monocytes. The link between chemokine signaling during early T. cruzi infection in heart cells remains to be investigated. Some of the DE piRNAs were computationally predicted to target CX3CL1 and CXCR2, respectively. Two known piRNAs, piR_016742 and piR_017716 were predicted to target an exon and the 3' UTR of CX3CL1, respectively ( Figure 3C). Similarly, piR_017716 was shown to target the 3' UTR of CXCR2 ( Figure 3D). We previously (B) Heatmap and hierarchical clustering was performed on control and test samples 0, 1, 3, and 6 h, respectively, using Euclidean distance measure and single linkage analysis. Each column represents one sample and each row represents one small non-coding RNA within the data set.
(C-E) Differential expressed putative and known piRNAs induced during parasite challenge of PHCF at 1, 3, and 6 hours respectively, belong to various transposable element subfamilies.
showed in Udoko et al. (Udoko et al., 2016) that T. cruzi dysregulates the gene expression profile of cardiac myocytes during the early phase of infection. In that report we showed that the parasite significantly dysregulates the expression of several profibrotic transcription factors including EGR1. Here, we show that a putative piRNA, npiR_167, which is DE during the early phase of parasite challenge is computationally predicted to target an exon region of EGR1 ( Figure 3E).

Validation of differentially expressed piRNAs and their targets
We validated the expression of the selected dysregulated piRNAs and their target genes (ICAM1, SMAD2, CX3CL1, CXCR2 and EGR1), which are implicated in T. cruzi infection and pathogenesis by qPCR. Total and small RNAs purified from T. cruzi challenged PHCF were used in our validation assays. Stem-loop PCR primers were designed to reverse transcribe each individual specific mature piRNA and U6 snRNA (housekeeping). The generated cDNA was quantified in the second step of the PCR by qPCR; a two-step RT-qPCR. Gene targets were evaluated using RT-qPCR as well. Our qPCR data show that during parasite challenge, the normalized fold change of CX3CL1 transcript showed a significant decrease at 1 h followed by a relative increase in expression at 3 h and 6 h, respectively ( Figure 4A). The targeting piR_016742 was shown to be upregulated at 1 and 3 h relative to 0, with a non-significant decrease at 6 h. This dynamic normalized expression negatively correlates to normalized expression levels of CX3CL1 transcript expression ( Figure 4B). CXCR2 showed a normalized decreased expression at 1 h followed by a normalized increase that was maximum at the 6 h time point ( Figure 4C). We previously reported that piR_017716, computationally predicted to target both CXCR2 and CX3CL1 was significantly decreased at 6 h (Rayford et al., 2022). More studies are needed to understand piRNA mediated regulation of CXCR2. The transcript levels of ICAM1 were significantly decreased at all time points. ICAM1 is uniquely targeted by 2 piRNAs, piR_016828 and npiR_17. Both piRNAs are significantly upregulated at 1 and 3 h relative to 0 h. Putative npiR_17 decreases at 6 h while piR_016828 remains significantly increased (Figures 4E, F). We previously showed that EGR1 expression during T. cruzi infection of cardiac myocytes increases to activate profibrotic molecules. Interestingly, EGR1 transcript expression showed a significant decrease at 1 and 6 h with a significant increase at 3 h relative to 0 h ( Figure 4G). The novel piRNA, npiR_167, predicted to target an exon of EGR1, showed a significant increase at 1 and 3 h, followed by a significant decrease at 6 h ( Figure 4H). Accumulating evidence has shown that AP-1 protein FOS can play a role in early T. cruzi pathogenesis by activating profibrotic genes (Moreira et al., 2017;Choudhuri and Garg, 2020). We predicted in Rayford et al. (Rayford et al., 2020) that several piRNAs can target FOS transcript in T. cruzi challenged PHCM. Here, we evaluated the FOS transcript expression and show that it remained significantly increased from 1 to 6 h compared to control ( Figure 4I), while the npiR_542 computationally predicted to target FOS was also significantly increased at all time points ( Figure 4J).

Biological interaction networks of differentially expressed piRNAs and their targets during early T. cruzi infection
To identify the molecular mechanisms that could be activated in PHCF and regulated by DE piRNAs induced during the early phase of cellular infection by T. cruzi, we used the GeneMANIA algorithm to construct two biological interaction networks based on adhesion molecule and chemokine expression (ICAM1, CX3CL1, CXCR2) and the AP-1 transcription factor family (FOS, SMAD2, EGR1). These two sets of query molecules were expanded to one degree of molecular protein-protein interactions and then connected with validated known and putative piRNAs. The resulting networks were visualized with Gephi. We show that through predicted targeting of CXCR2, piR_017716 may influence the activity of the durotactic mediator vasodilator-stimulated phosphoprotein (VASP), some molecules involved in fibrotic responses including CXCL2, CXCL3, and CXCL8. ICAM1 targeted by piR_016828 and npiR_17, and CX3CL1 targeted by piR_016742 and piR_017716 interacted with various molecules including fibrinogens, integrins, and other regulators of the inflammatory response such as STAT1 and IL2RA ( Figure 5). SMAD2 was predicted to be targeted by piR_016828 and piR_017716 and FOS by piR_542 and piR_018573. Putative piRNA npiR_167 is the sole piRNA predicted to target EGR1, which is directly connected to several AP-1 family proteins. Molecules interacting with FOS, SMAD2, and EGR1 are involved in TGF-b signaling through FOXH, SMAD3, and SMAD4. Epidermal growth factor (EGF) signaling may also be impacted through CAMK2A, ELK1, FOSB, JUN and JUND ( Figure 6).

Discussion
T. cruzi, the etiological agent of Chagas Disease, causes severe morbidity, mortality, and economic burden globally. Recent studies emphasize increased surveillance and screening of T. cruzi infection in the US, following increased detections of autochthonous infections (Dorn et al., 2007;Rassi et al., 2010;Beatty and Klotz, 2020;Lynn et al., 2020). Chronic infections have been associated with the development of cardiomyopathies including but not limited to myocardial hypertrophy, cardiac fibrosis, and cardiomegaly, among others. Though current research aims to understand gene expression regulation in cardiac tissue, few studies have been focused on the contribution of cardiac fibroblasts to infection and pathogenesis. Upon injury, cardiac fibroblasts differentiate into myofibroblasts to produce excessive amounts of extracellular and other profibrogenic molecules. This unregulated response can lead to the onset of T. cruzi-induced cardiomyopathies. However, the molecular mechanisms by which T. cruzi alters the gene expression profile in cardiac fibroblasts Early T. cruzi infection dysregulates expression of piRNAs and their targets. (A-J). Total RNA and small RNA were purified from PHCF challenge with invasive T. cruzi trypomastigotes. Target transcript expression was evaluated using RT-qPCR, while piRNA expression was validated using stem-loop PCR for cDNA synthesis and RT-qPCR for relative expression. Bar graphs show log 10 fold change ± SEM (n=3). Orange bar graphs denote relative mRNA target transcript expression while green bar graphs represent relative piRNA expression. Student's t-test and visualization performed using GraphPad Prism. ns = not significant, *p < .05, **p < .01, ***p < 0.001, ****p < 0.0001.

FIGURE 5
Network of differentially expressed piRNAs targeting adhesion and chemokine signaling molecules operating during T. cruzi challenge of PHCF. Biological interaction networks were created using predicted piRNA targets (ICAM1, CX3CL1, and CXCR2) as primary seed nodes to query pathway and interaction data sources out to one degree of biological interaction. Primary seed notes are denoted as orange circles. Expansion notes (blue) with a predicted binding target of DE piRNAs (green).
remains to be fully understood. Ongoing research seeks to understand alterations in gene regulation during parasite infection. Small noncoding RNAs such as miRNAs and piRNAs are potent modulators of gene expression and regulation, and thus have become important in understanding dysregulation of host gene expression during disease pathogenesis. In recent years, we and others have begun to investigate the role of small noncoding RNAs in T. cruzi pathogenesis. Ferreira et al. linked altered miRNA expression and their targets to various T. cruzi pathologies such as hypertrophy and fibrosis . The traditional sncRNA regulation paradigm proposes miRNAs directly regulate transcript and protein expression via binding to targets when complexed with RISC . In the new world of gene regulation, piRNAs have been suggested to regulate gene expression, transposon silencing, de novo methylation, chromatin remodeling, and translational activation when complexed with PIWI proteins (Ku and Lin, 2014;Rayford et al., 2021). Recently, we showed in infectious disease research that T. cruzi dysregulates the expression profile of host piRNAs computationally predicted to target profibrotic molecules including TGFB1, FOS, and NFATC2 in primary human cardiac myocytes (Rayford et al., 2020). Here, we hypothesize that T. cruzi dysregulates the expression of piRNAs in PHCF during the acute phase of infection. The dysregulated piRNAs can regulate the gene expression profile of profibrotic and inflammatory molecules that we and others reported (Suzin et al., 1987;Waghabi et al., 2005;Kayama et al., 2009;Udoko et al., 2016).
We challenged PHCF with the parasite at different time points and evaluated the piRNome employing an approach that results in lower error accumulation and sequencing bias, thereby improving the quality of our data as we recently described (Rayford et al., 2020). PHCF constitute a good model that closely mimics the natural infection. We anticipate that the data obtained from in vitro assays using PHCF will mimic changes in host piRNome profile induced by the parasite in a natural infection. The findings will guide studies geared towards understanding parasite mediated pathology in complex microenvironment in subsequent studies. In silico analysis showed that replicates clustered independently from each other and controls, indicating substantial alterations in the piRNA expression profile (Figures 1A, B, 2A, D). Though piRNAs are derived from a plethora of origins, we found that the significantly DE piRNAs could be classified into various TE families such as LINEs, LTRs, and class 2 DNA transposons ( Figures 1C-E). The pattern of piRNAs mapped to different transposable element subfamilies agrees with our observations in PHCM (Rayford et al., 2020). Over 85% of the DE piRNAs were derived from LTRs and LINEs, which suggests that these TEs may serve as important precursor sources of piRNA biogenesis in cardiac cells during T. cruzi infection. The high percentage of piRNAs originating from LTRs and LINEs suggests that they play a key role in the regulation of gene expression during the early phase of T. cruzi infection in PHCF compared to PHCM where the large majority of DE piRNAs were from LTRs (Rayford et al., 2020).
Our data show that the parasite induced significant DE of known and novel piRNAs at different time points (Figures 2A-G). The majority of the DE piRNAs that we report here are novel. NCBI BLAST search analysis of the putative piRNA sequences against the T. cruzi sequence database showed no significant alignments, suggesting that the piRNAs are host-derived. These novel putative piRNAs exhibit characteristics of canonical piRNAs that are yet to be annotated in the human genome (Supplementary Tables 1, 2).
We showed in our previous publication that miRANDA and RNA22 algorithms could be adapted for predicting piRNA targets (Rayford et al., 2020). We used miRANDA to computationally predict the piRNA gene targets. We created a list of genes most pertinent to T. cruzi infection and fibrosis by using Geneshot to search aggregate publications and gene expression repositories (Lachmann et al., 2019). From our list of piRNA target genes we identified ICAM1, SMAD2, CX3CL1, CXCR2, and EGR1 for further evaluation and validation in this report (Figures 3, 4). The selected genes have been suggested to play roles in T. cruzi infection and pathology (Blanton et al., 1991;Laucella et al., 1996;Huang et al., 1999;Michailowsky et al., 2004;Shen et al., 2019). We validated the relative expression levels of the selected piRNAs of interest and their gene targets (Figure 4). Several studies have linked ICAM1 as a driver of cardiac remodeling in human, rat, and mouse models (Luc et al., 2003;Davani et al., 2006;Fotis et al., 2012;Lino et al., 2019). Studies examining chronic T. cruzi infection suggest upregulation of various vascular adhesion molecules such as VCAM1 and ICAM1 (Laucella et al., 1996;Huang et al., 1999;Michailowsky et al., 2004). We found that during the acute phase of infection, ICAM1 transcript levels were significantly decreased at all time points following T. cruzi challenge. Conversely, piR_016828 expression was significantly increased at all corresponding time points. Similarly, npiR_17 significantly increased at 1 and 3 hours but significantly decreased at 6 hours. Our data suggest that these two piRNAs predicted to target ICAM1 may in part contribute to the downregulation of ICAM1 transcript observed during early infection. It will be interesting to evaluate the levels of the ICAM and the target piRNAs during the chronic phase of T. cruzi infection as a potential biomarker for cardiac remodeling. Few studies have examined the role of CX3CL1 in T. cruzi pathogenesis, however Zhang et al. observed that the expression of CX3CL1 is increased in liver fibrosis caused by Network of DE piRNAs targeting profibrotic transcription factors during PHCF challenged with T. cruzi. Biological interaction networks were created using predicted piRNA targets (SMAD2, FOS, and EGR1) as primary seed nodes to query pathway and interaction data sources out to one degree of biological interaction. Primary seed notes are displayed as orange circles. Expansion notes (blue) with a predicted binding target of DE piRNAs (green).
Schistosoma infection (Zhang et al., 2020). Here, we observed a significant decrease in CX3CL1 expression at 1 hour following parasite challenge with a significant increase at 3 and 6 hours. The expression of piR_016742 was significantly increased at 1 and 2 hours and significantly decreased at 6 hours. This inverse expression relationship between the piRNA and transcript suggests that piR_016742 may negatively regulate CX3CL1 expression. CXCR2 is a G-protein-coupled receptor that binds to interleukin 8 (IL8) with high affinity. This receptor also binds to CXCL1 and CXCL8 (Li and Frangogiannis, 2021). Previous studies show that CXCR2 activation can be cardioprotective or pathogenic, depending on the cardiovascular disease (Tarzami et al., 2003;Liehn et al., 2013). Schmitz et al. found that CXCR2 antagonism prevented T. cruzi trypomastigote-induced plasma leakage, suppressing the inflammatory response to the parasite (Blanton et al., 1991). During parasite challenge in PHCF, CXCR2 expression was significantly downregulated at 1 h but recovered at 6 h post infection. Interestingly piR_017716 showed nonsignificant increases at 1 and 3 h but significantly decreased at 6 h. This correlation implicates an inverse relationship and potential negative regulation of CXCR2 by piR_017716. EGR1 is a transcription regulator induced by growth factors, cytokines, and stress signals such as radiation, injury, or mechanical stress to transduce signals that activate downstream signaling cascades. Various fibrogenic diseases are associated with elevated levels of EGR1 and TGF-b. We and others show that T. cruzi infection upregulates the expression of EGR1 in various cell types (Mott et al., 2011;Udoko et al., 2016;Shen et al., 2019). In T. cruzi challenged PHCF, EGR1 expression was significantly decreased at the 1 h time point, followed by a significant increase at 3 h, which returned to a significant decrease at 6 h. The expression of npiR_167 does not appear to strongly correlate with EGR1 transcript levels. In another study, we showed that T. cruzi upregulated the transcript and protein levels of EGR1 in PHCM during the early phase of infection (Udoko et al., 2016), an observation which is different in parasite challenged PHCF. EGR1 upregulation in PHCM potentially leads to increased dysregulation of extracellular signals that can alter fibroblast response in cardiac tissue. Within the same period, the parasite induced an increase in the expression of FOS, an AP-1 transcription factor at all time points as we showed before with PHCM (Udoko et al., 2016). The novel piRNA shown to target FOS transcript, npiR_542 was also increased at all time points suggesting that this piRNA may function to stabilize the transcript for increased FOS protein expression. Since increased FOS expression has been associated with fibrogenic responses (Eferl et al., 2008;Nakagawa et al., 2016), its regulatory piRNAs can be developed to serve as fibrogenic biomarkers and/or therapeutic targets. Ongoing studies in the laboratory are geared towards determining the exact molecular mechanisms by which the DE piRNAs regulate gene expression to facilitate T. cruzi pathogenesis.
Selected DE piRNAs and their predicted targets ICAM1, SMAD2, CX3CL1, CXCR2, and EGR1 were used to generate networks that linked the piRNAs to their target genes and expanded them to one degree of biological interaction. Our in silico analysis showed that some of the differentially expressed piRNAs (known and novel) can target genes that we and others reported to be dysregulated during the early phase of T. cruzi infection. Two networks were generated based on chemokine signaling ( Figure 5) and profibrotic transcription factors ( Figure 6). In the chemokine network, CXCR2 had the most proteinprotein interactions and was connected to ICAM1 through CX3CL1 and integrin beta 2 (ITGB2). We observed that the piRNAs dysregulated by the parasite may have an influence on other chemokines and proteins previously implicated in T. cruzi infection and Chagas disease. VASP is a signaling pathway that regulates integrin-extracellular matrix interactions directly modulating the actin ultrastructure. Changes in VASP activity have been linked to cardiomyopathies (Pula and Krause, 2008). T. cruzi infection triggers substantial production of nitric oxide (NO), which has been shown to have protective and toxic effects on the host's immune system (Bergeron and Olivier, 2006). Integrin alpha-L chain (ITGAL) and ITGB2 form lymphocyte function-associated antigen 1 (LFA-1) and binds with ICAMs to facilitate intercellular adhesion of leukocytes. Integrins mediate lymphocyte migration into an infected tissue, and these cells are essential for regulating parasitemia, clearance, and inflammation (Acevedo et al., 2018). Mice treated with anti-LFA-1 displayed increased blood and tissue parasitemia, and quickly succumbed to infection . T. cruzi reduces interleukin 2 receptor alpha (IL2RA) expressed on the surface of immune cells, suppressing the immune response (Majumder and Kierszenbaum, 1996). Levels of pro-platelet basic protein (PPBP) and other peripheral blood mononuclear cells correlate with disease state in Chagas disease (Garg et al., 2016;Zago et al., 2018).
The DE piRNAs were also predicted to target sequences on genes that were within one degree of interaction with the AP-1 transcription factors evaluated in this study (Figures 5, 6). In the profibrotic transcription factor network, FOS had the most interactions. We and others have described alterations in the AP-1 gene and protein expression profiles during parasite infection. In addition to the primary AP-1 family proteins in the seed nodes, other proteins linked to them by one degree of biological interaction could be implicated in piRNA regulation; NAB1 and NAB2 are transcriptional cofactors connected to EGR1. Transgenic mice exhibiting cardiac-specific overexpression of NAB1, revealed that NAB1 inhibits cardiac growth in response to pathological stimuli in vivo. NAB1 overexpression inhibited pressure overload-induced hypertrophy (Buitrago et al., 2005). NAB2 directly binds to EGR1 and regulates EGR1 gene transcription. Additionally, NAB2 expression is induced by TGF-b, which is suggested to play an important role in T. cruzi infection and pathogenesis (Waghabi et al., 2007;Bhattacharyya et al., 2009;Waghabi et al., 2009). In dermal fibroblasts, it was found that expression of NAB2 correlated with persistent fibrotic responses (Bhattacharyya et al., 2009;Bhattacharyya et al., 2011). CAMK2A, connected to SMAD2, is a major mediator of calcium signaling in the heart and it was suggested that T. cruzi uses calcium signaling to facilitate host cell invasion (Benaim et al., 2020). Inhibition of CAMK2A was shown to reduce occurrence of cardiac arrhythmias for in vivo and ex vivo Chagas disease models (Santos-Miranda et al., 2021). Elevated levels of insulin-like growth factor binding protein 7 (IGFBP7) have been observed in patients with heart failure and correlate with worse disease outcomes (Januzzi et al., 2018). In contrast, a study on dermal fibroblasts exposed to T. cruzi soluble factor resulted in reduced levels of IGFBP7 (Mott et al., 2011). SMAD2, a member of the SMAD protein family, is a signal transducer and transcriptional modulator that mediates multiple signaling pathways. The phosphorylated version of SMAD2 is downstream of TGFb signaling, thus it regulates multiple cellular processes. Accumulating evidence establishes that SMAD/TGFb signaling is required for T. cruzi infection and pathogenesis (Udoko et al., 2016;Ferreira et al., 2019;Silva et al., 2019). Accordingly, these selected DE piRNAs highlighted in our manuscript and their gene targets can be suggested to play significant roles in cardiopathogenesis as summarized in Table 1. These piRNAs can be developed as potential therapeutic targets against the onset of T. cruzi induced cardiomyopathies.
Taken together, our data show that during the early phase of T. cruzi challenge of PHCF, the parasite induces DE of piRNAs (both known and novel) that can target and regulate the expression of important profibrotic and inflammatory molecules. This finding is significant considering the extensive regulatory roles of piRNAs. This report showing that a pathogen can dysregulate host piRNA expression, leading to enhanced understanding of parasite-induced gene regulation, which is essential for the identification of biomarkers and the development of molecular intervention strategies during the early phase of T. cruzi infection.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.