Expression of Extracellular Vesicle PIWI-Interacting RNAs Throughout hiPSC-Cardiomyocyte Differentiation

Extracellular Vesicles (EV) play a critical role in the regulation of regenerative processes in wounded tissues by mediating cell-to-cell communication. Multiple RNA species have been identified in EV, although their function still lacks understanding. We previously characterized the miRNA content of EV secreted over hiPSC-cardiomyocyte differentiation and found a distinct miRNA expression in hiPSC-EV driving its in vitro bioactivity. In this work, we investigated the piRNA profiles of EV derived from key stages of the hiPSC-CM differentiation and maturation, i.e., from hiPSC (hiPSC-EV), cardiac progenitors (CPC-EV), immature (CMi-EV), and mature (CMm-EV) cardiomyocytes, demonstrating that EV-piRNA expression differs greatly from the miRNA profiles we previously identified. Only four piRNA were significantly deregulated in EV, one in hiPSC-EV, and three in CPC-EV, as determined by differential expression analysis on small RNA-seq data. Our results provide a valuable source of information for further studies aiming at defining the role of piRNA in the bioactivity and therapeutic potential of EV.


INTRODUCTION
The secretion of membrane vesicles is an evolutionarily conserved biological phenomenon observed across organisms, that mediates intricate processes of cellular regulation and coordination essential for survival (Deatheragea and Cooksona, 2012). These extracellular vesicles (EV) are nano and micro-sized particles without a functional nucleus, that transport a cargo of complex biomolecules involved in intra-and intercellular communication (Van Niel et al., 2018). Although EV can carry an assortment of genetic material, including both DNA and RNA, they are mostly enriched in small non-coding RNAs (ncRNA) (Nolte-'t Hoen et al., 2012). ncRNAs mediate evolutionarily conserved silencing pathways by regulating mRNA stability or translation, or by targeting epigenetic modifications to specific regions of the genome, and are involved in the control of gene regulatory programs (Statello et al., 2020) in development and disease (Li et al., 2019;Rayford et al., 2021). Multiple ncRNA species have been identified in EV, including microRNA (miRNA), ribosomal RNA (rRNA), transfer RNA (tRNA), long non-coding RNA (lncRNA), PIWI-interacting RNA (piRNA), small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA) (Nolte-'t Hoen et al., 2012;Kaur et al., 2018). Though there is no consensus on their relative percentage, several studies suggest miRNA accounts for 40%-70% of ncRNA in EV (Huang et al., 2013;Yuan et al., 2016). In opposition, other authors argue that most protocols for small RNA library construction use adaptor ligations that favor miRNA, thus introducing a significant bias in biotype identification (Raabe et al., 2014;Wei et al., 2017). Regardless of the reason, miRNAs are the most extensively studied ncRNA, and only a few studies focusing on other EV ncRNA molecules have been performed so far.
Thousands of piRNA have been identified across the mammalian genome. They mostly arise from TEs but can also originate from long intergenic non-coding RNAs (lincRNAs), noncoding, nonrepetitive genomic elements (piRNA clusters), and in rarer cases, from protein-coding genes (Sarkar et al., 2014). PIWI proteins and piRNA were initially thought to be restricted to animal germline (Girard et al., 2006) and stem cells, and their presence in somatic tissues is not well characterized, with controversy remaining around their existence and function. Several studies on piRNA support its existence in the soma (Yan et al., 2011;Rajasethupathy et al., 2012;Ross et al., 2014), with tissue-specific expression patterns, and have demonstrated its functional role in epigenetic reprogramming, regulation of transcription, translation, and development (Wu et al., 2020). Nevertheless, as of yet, few studies have shown that interfering with somatic piRNA or with the genes involved in its biogenesis induces majorly altered phenotypes (Jones et al., 2016), in opposition to what occurs in germ lines in which mutations in PIWI-family genes causes significant defects in germline development, namely infertility, across species (Schmidt et al., 1999;Carmell et al., 2007).
Groundwork on the identification of piRNA cargo of EV and its part in intercellular communication has been emerging but it is rather preliminary in comparison with the extensive portrayal of EV-miRNA.
In a previous study, we characterized the bioactivity and miRNome of human induced pluripotent stem cells (hiPSC)and hiPSC-cardiomyocytes (hiPSC-CM)-derived EV and found distinct EV potency and miRNA expression patterns during differentiation (Louro et al., 2022). Moreover, we corroborated the upregulation of stemness maintenance miRNA clusters in hiPSC-EV and showed that hiPSC-EV targeted the PI3K/AKT pathway, a crucial network in heart physiopathology, by mediating PTEN suppression and increasing AKT activity. Remarkably, it has been reported that piRNAs may also influence the PI3K/AKT signaling pathway (Rajan et al., 2014;Vella et al., 2016), and could thus have an important role in angiogenesis, proliferation and survival (Rajan et al., 2014). In this context, identification of piRNAs in EV secreted by hiPSC, cardiac progenitor cells (CPC) and hiPSC-CM could provide further insights into the mechanisms underlying EV bioactivity and uncover new therapeutic approaches for cardiovascular diseases.

RNA Isolation
Total RNA was isolated from EV samples (N = 3 EV preparations from 3 independent cultures of four cell populations, in a total of 12 EV samples) using Norgen Biotek Exosomal RNA Isolation Kit (Cat.58000), according to the manufacturer's instructions.

Library Preparation and Small RNA Sequencing
Next-generation sequencing of small RNAs was performed by Norgen Biotek Corp, as described in (Louro et al., 2022). Briefly, small RNA library preparation was performed by PCR amplification using Norgen Biotek Small RNA Library Prep Kit (Cat. 63600). Library quality control, including estimation of library size and concentration, was achieved using Bioanalyzer. Sequencing was performed using Illumina's NextSeq 500 sequencing platform.

Length Frequency and Small RNA Biotype Distribution Analysis
Read lengths and biotype count files generated by the small RNAseq pipeline on the Genboree Workbench were used for length frequency and biotype distribution analysis.

Analysis of Differentially Expressed piRNAs
Data filtering, normalization and differential expression (DE) analysis were conducted using the EdgeR (Robinson et al., 2009) package (version 3.14) within the R environment (version 4.1). Data was filtered based on a minimum of 5 counts per million and normalized using the trimmed mean of M-values (TMM) method (Robinson and Oshlack, 2010). DE analysis was performed using an exact test for differences between each two groups. False discovery rate was adjusted through the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995). DE was considered significant for |log 2 fold change| ≥ 1, with a p-value ≤ 0.05 and FDR ≤ 0.05. Principal component analysis (PCA) was performed using the DESeq2 (v1.34.0) R package (Love et al., 2014). Soft clustering was implemented with R package Mfuzz (v2.42.0) (Futschik et al., 2005;Kumar and Futschik, 2007). R packages ComplexHeatmap (v2.10.0) and EnhancedVolcano (v1.12.0) were used to plot heatmaps and volcano plots, correspondingly.

RT-qPCR Validation of Differentially Expressed piRNAs
Differentially expressed piRNAs were validated by reverse transcription quantitative real time polymerase chain reaction (RT-qPCR). Briefly, RNA samples were reverse transcribed with TaqMan ™ MicroRNA Reverse Transcription Kit (Applied Biosystems ™ , Thermo Fisher Scientific) and Custom Taqman ™ Small RNA Assays (Applied Biosystems ™ , ThermoFisher Scientific), according to the manufacturer's protocol. RT-qPCR was performed using the TaqMan ™ Fast Advanced Master Mix (Applied Biosystems ™ , ThermoFisher Scientific) on a LightCycler 480 Instrument II 384-well block (Roche) in the following cycles: UNG activation at 50°C for 2 min; polymerase activation at 95°C for 20 s; 40 cycles of amplification with denaturation at 95°C for 3 s, and annealing/extension at 60°C for 30 s. The Cycle threshold (Ct) was determined using the LightCycler 480 software (Roche). Stem-loop primers used for cDNA synthesis and RT-qPCR are available in Supplementary  Table S1. Two endogenous controls (hsa-miR-103a-3p and U6) were selected based on a comprehensive assessment of expression stability in our small RNA-seq data. Relative changes were analyzed using the ΔΔCt method (Livak and Schmittgen, 2001) (CPC-EV vs. hiPSC-EV, CMi-EV vs. CPC-EV).

RESULTS
piRNA Detection in Extracellular Vesicles piRNA expression was analyzed by small RNA-seq in three biological replicates of EV derived from four distinct cell populations: hiPSC (hiPSC-EV), CPC (CPC-EV), immature cardiomyocytes (CMi-EV), and cardiomyocytes matured by metabolic modulation (CMm-EV), along a continuous differentiation process ( Figure 1A). Analysis of length distribution on reads used for alignment, revealed the presence of three peaks, consistent for all the samples ( Figure 1B). The major peak represented mature miRNA, approximately 22 nt long, while two minor ones represented piRNA, with a size range of 24-34 nt ( Figure 1B). Nearly 60% of the mapped reads matched to miRNA, and the remaining reads to other non-coding RNAs. Filtering of the mapped reads based on piRNA size resulted in the removal of 65%-84% of reads ( Figure 1C).
After an additional filtering step based on mapped biotype for the removal of other non-coding RNAs in the same size range, only 2%-12% of the initial mapped reads remained. The experimental replicates behaved consistently throughout the filtering steps. Overall, piRNA represented the second most abundant small RNA biotype in 3 EV populations, mapping to 7.2% ± 1.1% of sequences in hiPSC-EV, 11.9% ± 1.5% in CPC-EV and 6.8% ± 1.7% in CMi-EV. In CMm-EV piRNA corresponded only to 2.1% ± 0.2% of mappings, being the third most frequent small RNA biotype ( Figure 1D).

Extracellular Vesicles piRNA Expression Profiles Along hiPSC-CM Differentiation
We evaluated the piRNome profile of EV-derived from four cell populations. A total of 293 piRNA were identified in EV samples, and after filtering to a minimum of 5 counts per million in at least 3 samples (the smallest group size), only 37 piRNA were considered reliable for further analysis (Supplementary File S1-TMM normalized counts). The global piRNA expression differed between EV groups, with stage specific patterns emerging ( Figure 2A). Fuzzy clustering was employed to analyze the dynamic expression variations throughout hiPSC-CM differentiation. In general, three broad piRNA clusters were identified ( Figure 2B): one cluster showed low expression in hiPSC-EV and consistently increased throughout cardiomyocyte differentiation and further maturation; a second cluster declined upon hiPSC commitment to the cardiac fate and remained stable from CPC-EV to CMm-EV; and the third cluster was only upregulated in CPC-EV. piRNA representing each of these clusters can be found in Supplementary Table S2.

Differentially Expressed piRNA in hiPSC-CM Derived Extracellular Vesicles
Hierarchical clustering on the sample-to-sample distances from normalized gene counts revealed limited differences among EV populations (Figure 3). hiPSC-EV and CPC-EV clustered together within their corresponding group, while CMi-EV and CMm-EV replicates were assorted in a third cluster, with the exception of CMm-EV replicate 1 that showed the highest difference and did not cluster with any of the other samples. This replicate was maintained, and it did not interfere with DE analysis.
These results were confirmed by PCA. Principal component 1, explaining 28% of the variation, separated hiPSC-EV and CPC-EV from CM-EV, while principal component 2, explaining 14% of the variation, distinguished hiPSC-EV from all other groups ( Figure 3B).
Globally, only 4 piRNA were found to be differentially expressed ( Table 1, Supplemental file 2) at a |log 2 fold change| ≥1 and -log 10 FDR ≥1.3. However, only piR-36770 and piR-30597 were consistently validated by RT-qPCR ( Supplementary  Figures S1A,B). In fact, the magnitude of the log 2 fold change obtained by both techniques differed greatly (Supplementary Figure S1C), although the same tendency was observed for all piRNA probed.

Extracellular Vesicles piRNA Mapping
Differentially expressed piRNA identified both by small RNA-seq and RT-qPCR were mapped to repeat elements contained in transposons using piRNAQuest. While no association between piR-30597 and these elements was found, piR-36770 mapped to numerous LTR/ERV-1 repeats in several chromosomes (Supplementary File S3).

DISCUSSION
RNA-carrying vesicles are released by cells into the extracellular space, allowing for the horizontal transfer of genetic material. This mechanism of paracrine communication allows the dissemination of messages that can regulate the function of target cells (Valadi et al., 2007). Several types of short ncRNA were found to be particularly enriched in these vesicles, including miRNA, tRNA, piRNA, and snoRNA. Numerous studies have highlighted the role of ncRNA in the heart (Bian et al., 2021;Martyniak et al., 2021;Bencun et al., 2022), in parallel with the growing evidence on EV as key players in cardiac diseases and heart regeneration (Adamiak and Sahoo, 2018;Santoso et al., 2020;Lima Correa et al., 2021).
Here, we performed an exploratory study into the piRNA content of EV released during the differentiation of cardiomyocytes from hiPSC. In a previous investigation, we established a method for EV isolation and characterization from conditioned culture media and uncovered the EV miRNome from these cell sources (Louro et al., 2022). The pathophysiological role of piRNA in the heart is still poorly understood and has been the focus of recent studies (Rajan et al., 2014(Rajan et al., , 2016Zeng et al., 2021). However, very few have focused on the EV piRNome (Wang et al., 2020).
Herein, we observed that piRNA were the second most abundant small RNA species in EV from undifferentiated and early stages of hiPSC-CM differentiation ( Figure 1D). Our results indicate that piRNA transcripts are present in CM-EV, thus providing evidence of the existence of somatic piRNA. Overall, piRNA were the second most abundant small RNA biotype in hiPSC-EV, reaching a peak in CPC-EV, decreasing in CMi-EV, and considerably declining in CMm-EV ( Figure 1D), suggesting an association with the pluripotent state and early stages of cell differentiation, hypothesis that requires further validation.
Employing a soft clustering algorithm, we identified three different patterns of piRNA expression ( Figure 2B) with similar trends to some of the ones identified by La Greca and colleagues when characterizing the piRNome of hiPSC, early mesoderm progenitor cells, and cardiomyocytes (Greca et al., 2020). The authors identified 447 piRNA transcripts, with 8 distinct patterns, of which 3 consistently decreased over differentiation, and 3 consistently increased until later maturation stages. Compared to cells, EV had a much lower abundance of piRNA (only 37 identified transcripts), a discrepancy that indicates that piRNA accumulation in EV may be limited, in comparison to miRNA.
Even though stage-specific piRNA patterns can be inferred, the four EV populations presented low variation between each other ( Figures 3A,B). In fact, the hierarchical clustering showed that samples were quite similar to each other, the highest difference occurring in hiPSC-EV, which clustered independently from all other EV samples. These results suggest that in opposition to the miRNA cargo of EV, which mimics the content and allows a strict characterization of the differentiation stage of its parent cell, the EV piRNA cargo is, in general, less correlated with cell differentiation profiles. Nevertheless, differences in piRNA profiles between hiPSC-EV and CPC-EV suggest there may be a function for these piRNA in fine-tuning transcript levels during cardiac commitment, piR-36770, overexpressed in hiPSC-EV, mapped to several Human endogenous retroviral H (HERVH) elements, implying a possible role of this piRNA in the regulation of the pluripotent state. In mammals, several ERVs are known to have stage-specific expression during embryonic development, and to be crucial for stem cell self-renewal (Santoni et al., 2012;Ohnuki et al., 2014;Pontis et al., 2019). HERVH is a primate-specific ERV family and one of the most predominant ERVs in pluripotent stem cells. In fact, HERVH repeats contained in lincRNAs exhibit a stem cell-specific expression pattern (Kelley and Rinn, 2012), and are required for pluripotency maintenance (Lu et al., 2014). In hPSC, the HERVH subfamily has numerous binding sites for and recruits several pluripotency transcription factors such as NANOG, OCT4, SOX2, and KLF4 (Kunarso et al., 2010;Ohnuki et al., 2014;Wang et al., 2014;Pontis et al., 2019), crucial for reprogramming towards iPSCs and maintenance of pluripotency and self-renewal features. Moreover, RNA interference knockdown of individual HERVH-derived lncRNAs in iPSC results in loss of pluripotency (Ohnuki et al., 2014). piR-30597 and piR-33415 had already been detected by Vella and colleagues in adult CPC obtained from right atrial biopsies, and cultured as 3D clusters, known as cardiospheres (Vella et al., 2016). piR-30597 was detected by RT-qPCR in EV samples, but no relevant interactions were found for it. To the best of our knowledge, piR-31108 had not been described in the cardiac context. piRNA sequences are stored in databases with distinct nomenclature and accession codes, a lack of consensus may introduce a bias in literature revision, as publications may address the same piRNA by a different alias. To overcome this limitation, we performed a literature search for the differentially expressed piRNA with accession codes for NCBI (https://www. ncbi.nlm.nih.gov), ENA (https://www.ebi.ac.uk/ena), piRNABank (http://pirnabank.ibab.ac.in), piRBase (www. pirnadb.org; version 1.8.0), and piRNAQuest (http:// bicresources.jcbose.ac.in/zhumur/pirnaquest/) databases.
Since we had previously identified the PI3K/AKT signaling as one of the pathways involved in hiPSC-EV bioactivity, we were prompted to search for piRNA in our dataset that targeted this pathway. Lucchinetti and colleagues showed that inhibition of LINE-1 retrotransposon expression in the heart activated the PI3K/AKT signaling (Lucchinetti et al., 2006). Similar results were obtained by Vella, which showed that several piRNA upregulated in cardiospheres targeted LINE retrotransposons throughout the genome, of which LINE-1 was most targeted, and that reduction of L1 expression was accompanied by activation of the AKT pathway (Vella et al., 2016). However, none of the piRNA in our dataset mapped to LINE-1, suggesting another mechanism of action for these piRNA. This reasoning is supported by Du and colleagues who described the role of piR-31115 -also present in our dataset (Supplementary File S2)-in renal cell carcinoma progression by activation of epithelial-mesenchymal transition via the PI3K/AKT signaling, independently of L1 repression (Du et al., 2021), as this piRNA also does not map to LINE-1. Some limitations should be considered when interpreting the results of this study. First, small RNA-seq was performed only on EV and not on their respective cell biofactories. As such we are not able to conclude about any particular enrichment of the identified piRNA in EV versus their parent cells. Second, only an EV subfraction, enriched in small-EV was analyzed. We did not study other EV subpopulations nor compared their small ncRNA content. Third, besides piRNA, total small RNA of the size range of 24-34 nt may contain contaminants of abundant endogenous ncRNAs (such as tRNAs or rRNAs), and degradation products of longer RNAs. Indeed, total small RNA-seq reads provide information of the entire profile of expressed small RNAs, but not on which Argonaute protein a sequence associates with (Iwasaki et al., 2017). Therefore, future studies should combine total small RNA-seq data with small RNA-seq information derived from small-RNAs co-immunoprecipitated with PIWI proteins, technique that as far as we know has not been reported in the literature for EV. Finally, a great number of very low count tags accounted for most of our piRNA dataset and upon filtering based on the number of reads, few piRNA were left to perform DE analysis. In turn, this resulted in few piRNA differentially expressed, which compromised pathway enrichment analysis. As such, though we identified piRNA likely involved in particular cellular transitions, we were unable to disclose the mechanisms underlying their activity.
To conclude, we identified for the first time, global piRNA signatures in EV secreted over hiPSC-CM differentiation. The data reported in this work are a valuable source of information for further studies aiming at understanding the role of piRNA in the bioactivity of EV, in order to clarify their therapeutical relevance.

AUTHOR CONTRIBUTIONS
AL and MS conceived the study. AL and MP performed data acquisition. AL and NV performed data analysis and interpretation. AL wrote the manuscript. All authors revised and approved the manuscript. MS, II, and PG-A coordinated the study.  Supplementary File S1 | TMM Normalized counts for the 37 piRNA transcripts.
Supplementary File S2 | Differential expression analysis between the four EV populations.