Modeling Transposition of the Great Arteries with Patient-Specific Induced Pluripotent Stem Cells

The dextro-transposition of the great arteries (d-TGA) is one of the most common congenital heart diseases. To identify biological processes that could be related to the development of d-TGA, we established induced pluripotent stem cell (iPSC) lines from two patients with d-TGA and from two healthy subjects (as controls) and differentiated them into endothelial cells (iPSC-ECs). iPSC-EC transcriptome profiling and bioinformatics analysis revealed differences in the expression level of genes involved in circulatory system and animal organ development. iPSC-ECs from patients with d-TGA showed impaired ability to develop tubular structures in an in vitro capillary-like tube formation assay, and interactome studies revealed downregulation of biological processes related to Notch signaling, circulatory system development and angiogenesis, pointing to alterations in vascular structure development. Our study provides an iPSC-based cellular model to investigate the etiology of d-TGA.


Introduction
Congenital heart disease (CHD) represents almost one-third of all congenital diseases [1] and occurs in nearly 1% of live births [2]. The dextro-transposition of the great arteries (d-TGA) is one of the most common and severe types of CHD [1,3] with a prevalence of 20-33 per 100,000 newborns and an incidence of 5-8% of all CHD [4,5]. While d-TGA is most frequent in males (ratio 2:1), there is no correlation with any chromosomopathy [6]. d-TGA is generally reported with Heterotaxy syndrome, but the etiology and morphogenesis of TGA are still largely unknown. Nonetheless, several genes, such as Pitx2 and Tbx2, have been evaluated in animal models [3,7]. During heart development, the cardiac tube expands and begins a morphogenetic looping process that ultimately brings the posterior part of the tube to a rostral position, dorsal to the outflow tract (OFT) [8]. During mid-development, the OFT is remodeled into separate pulmonary and aortic arteries, which involves direct interactions between myocardium, endocardium and neural crest cells [9]. Subendocardial extracellular matrix (ECM) provides physical support and signaling regulation to cells and regulates tissue morphology during embryonic development [10]. The correct rotation of the myocardial wall of the OFT is necessary for the normal position of the great arteries. A high proportion of Perlecan (Hspg2)-null mouse embryos present with complete TGA [11], suggesting the importance of the ECM and basement membranes in cardiac morphogenesis. d-TGA can also be induced by retinoic acid treatments [12], directly involving the Tbx2-Tgfβ2 pathway [7].
Several studies have used animal models of TGA to better understand disease mechanisms, and there have only been a few studies reported in humans (reviewed in [3]). The emergence of human induced pluripotent stem cells (iPSC) [13] has raised the possibility of more robust disease modeling [14]. Due to the importance of endothelial cells (EC) in the vascular system, several studies have focused on designing robust protocols to differentiate iPSC into endothelial cells (iPSC-ECs) [15,16].
In the present study, we sought to generate a d-TGA cellular model based on the use of patient-specific iPSCs. Specifically, we generated two iPSC lines from patients with d-TGA, which were differentiated into iPSC-ECs. Several genes were identified by RNA-seq as potential players in the etiology of this congenital defect. We show that iPSC-ECs from patients with d-TGA have an impaired capacity to develop vascular networks concomitant with changes in the transcriptome, as assessed by bioinformatics analysis. Specifically, gene ontology (GO) biological processes related to embryonic morphogenesis/development, tube development and angiogenesis, together with Notch signaling pathways, were downregulated in EC derived from patients with d-TGA.

Generation of Endothelial Cells from Differentiated iPSCs
We generated iPSC lines from two patients with d-TGA (TGA1 and TGA2) and also one control iPSC line from a healthy individual (Ctrl1), which were cultured under feederfree conditions. We used the commercially available iPSC line IMR90-4 as a second control (Ctrl2). The cells were maintained undifferentiated and displayed typical morphology ( Figure 1A) and normal karyotype ( Figure S1). We differentiated the iPSC lines to ECs by the addition of different growth factors ( Figure 1B), following a published protocol [15]. We then isolated cells that were double-positive for CD31 and neuropilin-1 (NRP-1), a coreceptor of vascular endothelial growth factor (VEGF), at day 12, using fluorescence-activated cell sorting (FACS), which ensured an EC population. The results of the sorting analysis revealed a similar but low efficiency of EC differentiation in all four iPSC lines: FACS efficiency was 0.93 ± 0.69 (Ctrl1), 2.8 ± 2.08 (Ctrl2), 0.57 ± 0.23 (TGA1) and 0.6 ± 0.4 (TGA2). Isolated ECs showed a characteristic cobblestone morphology by brightfield microscopy ( Figure 1C) and were able to grow up to four passages. During endothelial differentiation of iPSC, we obtained RNA samples from each cell line at days 0, 5, 12 and 19 for gene characterization by qPCR. As expected, the levels of the pluripotency gene OCT4 decreased over time in all the lines tested. By contrast, the expression of the specific endothelial gene CD31, vascular-endothelial cadherin (CDH5) and receptor tyrosine kinase (TIE2) increased over time ( Figure 1D). All iPSC lines exhibited an increase in CD31 and TIE2 expression over time, indicating the differentiation of iPSC into EC. The expression profile of OCT4 and CD31 was similar in all the iPSC lines; however, some differences were observed when CDH5 and TIE2 were analyzed ( Figure 1C).

RNA-Seq Analysis of Endothelial Cells from Differentiated Control and d-TGA iPSC Lines
We next performed RNA-seq and bioinformatics analysis of sorted double-positive cells at day 20 of differentiation of the four iPSC-EC lines. Despite the known variability that exists between individuals [17,18] and the different origin of the iPSC lines used, the iPSC-EC samples from controls and patients with d-TGA were well separated by the principal component analysis (Figure 2A). These results support the notion that despite the interindividual genotypic differences, there were robust genetic differences between control and TGA lines. (C) Representative brightfield images of iPSC-ECs at day 19. Images were recorded using a 10× objective, scale bar = 10 µm. (D) Expression of genes involved in EC differentiation measured by qPCR in control (blue, • Ctrl1, Ctrl2) and d-TGA lines (red, TGA1, TGA2). Samples were obtained at days 0, 5, 12 and 19 of EC differentiation. Data were normalized to GAPDH expression and are represented as mean ± SD. Three independent differentiations were performed for each iPSC line. ** p < 0.01, by Mann-Whitney U test. Scale bar = 10 µm. A total of 4672 genes were identified by bioinformatics analysis, and 4471 genes had similar expression values in d-TGA and control groups. Identified genes with fewer than 150 counts were removed. We generated a volcano plot to visualize differentially expressed genes (DEG) (p-value ≤ 0.05 and fold change ≥ 2.0) in iPSC-ECs from d-TGA and control lines ( Figure 2B). A total of 201 DEG were identified, of which 62 were upregulated in d-TGA lines and 139 were downregulated. Regarding the upregulated genes, we identified several genes associated with the extracellular matrix (ECM), such as collagen type I alpha 1 (COL1A1), COL1A2 and laminin subunit alpha 4 (LAMA4). We also found an increase in the expression of the transcription factors SIX1 and TBX2, which have important roles in heart morphogenesis [8,19]. A member of the Fox transcription factor family (FOXC1) with an important role in vascular development was also upregulated in iPSC-ECs from patients with d-TGA.
Of the 139 downregulated genes in the d-TGA group, we identified several genes involved in vascular compartment and developmental processes. For example, specific endothelial markers, including PECAM1, CDH5, TIE1 and VEGF receptor 2 (KDR) [20] were all downregulated in iPSC-ECs from patients with d-TGA, as well as genes involved in Notch signaling, including jagged canonical Notch ligand 2 (JAG2) and Delta-like canonical Notch ligand 4 (DLL4). In addition, genes involved in gap junctions, such as gap junction protein alpha 4 (GJA4) and GJA5, which are necessary to establish cell communication, were also downregulated. Likewise, vascular stabilizing genes related to the ECM, such as collagen type IV alpha 1 chain (COL4A1), COL4A2 and HSPG2, were downregulated. A complete list of the upregulated and downregulated DEG is shown in Supplementary Tables S1 and S2, respectively.

iPSC-ECs from Patients with d-TGA Show Alterations in Signaling Pathways
In agreement with the gene expression profiling, GO analysis of DEG demonstrated different biological processes associated with control and d-TGA iPSC-ECs. A list of all significantly overrepresented GO biological processes of upregulated and downregulated genes, as a function of their p-value, is shown in Supplementary Tables S3 and S4, respectively. Many of the GO processes identified in upregulated genes were involved with developmental processes, including animal organ development and morphogenesis, embryonic development and morphogenesis, negative regulation of developmental process and anatomical structure development (Table 1). Table 1. Gene ontology biological processes associated with differentially upregulated genes in iPSC-ECs from patients with transposition of the great arteries versus healthy controls. As a complementary representation, we used a REVIGO treemap to reduce the number of GO terms into 'clusters' and 'superclusters' of related terms, as described [21], where GO terms are visualized with different colors in rectangles representing the p-value ( Figure 2C). The 304 GO biological processes identified were mainly grouped into two superclusters: 'negative regulation of multicellular organismal process', which was the supercluster most represented, and 'animal organ development'. We performed the same enrichment analysis on the downregulated genes; 279 GO terms were identified (p-value < 0.05) and the most representative are listed in Table 2. Of note, many processes were related to blood vessel development and morphogenesis, tube development and morphogenesis, negative regulation of developmental process and regulation of blood vessel endothelial cell proliferation involved in sprouting angiogenesis. These results suggest that the downregulated genes identified in iPSC-ECs from patients with d-TGA are of particular relevance for understanding the mechanisms of d-TGA in our system.

GO_ID
We also created a treemap to group the GO biological processes, which were mostly grouped into the 'regulation of multicellular organismal process' and 'circulatory system development' superclusters ( Figure 2D).
Several GO biological processes were identified when the analysis was performed with up-and downregulated genes. GO analysis of genes downregulated in iPSC-ECs from patients with d-TGA versus controls revealed significant enrichment of genes involved in the circulatory system. We next conducted an interactome analysis to study the interaction network between the targets of DEG in iPSC-ECs from patients with d-TGA versus controls (Figure 3). GO biological processes were identified with gene-set enrichment analysis, and the interaction network was visualized with Cytoscape software v3.7.2.
When targets of upregulated genes were analyzed, biological processes related to anatomical structure morphogenesis and development, regulation of the developmental process, embryonic morphogenesis and development were enriched in iPSC-ECs of patients with d-TGA (Table 3).
When targets of downregulated genes were analyzed, biological processes related to cellular and animal organ development processes were enriched in iPSC-ECs of patients with d-TGA (Table 4). In agreement with the aforementioned results, processes related to circulatory system development, blood vessel development and morphogenesis, and Notch signaling were also identified.

Capillary-like Tube Formation Processes Are Altered in iPSC-ECs from Patients with d-TGA
The results of the bioinformatics analysis suggested that iPSC-ECs from patients with d-TGA might show defective blood vessel development and morphogenesis. To test this, we performed a functional in vitro assay with iPSC-ECs to evaluate their capacity to form a capillary-like tube network, as a readout of angiogenesis. Cells were seeded onto Matrigel in EGM-2 medium, and tube formation was evaluated microscopically after 3 h. Results showed that iPSC-ECs from patients with d-TGA were unable to correctly form tubular structures measured as total loops, total tube length, number of total tubes and branching points ( Figure 4A,B), with significant differences in other parameters, including covered area (%), mean loop area and mean loop perimeter ( Figure 4B). By contrast, the total net abundance was significantly higher in iPSC-ECs from patients with d-TGA, and no differences were found when the mean tube length was measured. These results are in line with the interactome analysis, where biological processes related to the regulation of angiogenesis were downregulated (Figure 3).

Discussion
d-TGA is the sixth most prevalent CHD from a total of 27 anatomical subtypes, as reported in a systematic review [22]. The advent of human iPSC technologies in 2007 by Yamanaka et al. [13] has been a major breakthrough for disease modeling and for therapeutic purposes [14], and detailed gene analysis of iPSC lines from patients with CHD is a promising strategy to study the mechanisms involved in the onset of cardiac diseases [23].
We studied the gene expression patterns in iPSC-ECs from patients with d-TGA, the deregulation of which may be related to the onset of disease. Two iPSC lines were successfully generated from patients with d-TGA and were differentiated into ECs. Despite the low efficiency of EC differentiation, the FACS-sorted CD31+/NRP-1+ double-positive cells proliferated and expressed endothelial markers. We used bioinformatics analysis to identify genes differentially expressed in iPSC-ECs from patients with d-TGA. One of these genes was the transcription factor Tbx2, which was overexpressed in ECs from patients with d-TGA. Tbx2 expression is closely related to heart development, as it is expressed in the heart tube, and it has been described to arrest cardiac development at looping [8]. It has also been observed to be highly expressed in the early stages of development, but it decreases over time [24]. Indeed, Tbx2 expression is involved in other aspects of heart looping but is not related to chamber growth [25]. Another gene upregulated in iPSC-ECs from patients with d-TGA was FOXC1, which has been described to promote angiogenic activity in vascular endothelial cells that is determined by the correct balance between pro-angiogenic FOXC1 activity in the endothelium and anti-angiogenic FOXC1 activity in the surrounding cells [26]. Two members of the histone family, HIST1H3C and HIST1H2BB, were also upregulated in iPSC-ECs from d-TGA patients. Altered expression of different histones has been related to vascular diseases, such as abdominal aortic aneurysm, which is characterized by the degradation of the elastic media and remodeling of the aortic ECM [27]. Different genes upregulated in iPSC-ECs from d-TGA (LAMA4, COL1A1, COL1A2 and COL6A1) are involved in ECM organization (GO:0030198). These data suggest that the expression of different genes related to ECM is altered in iPSC-ECs from patients with d-TGA. HSPG2, a glycosylated protein component of the ECM, was downregulated in iPSC-ECs from d-TGA patients. Interestingly, Hspg2-null embryos presented with d-TGA [11], supporting the important role of this gene and genes related to ECM in the development of the disease.
Genetic polymorphisms in GJA4 have been previously related to coronary artery disease [28]. This gene was downregulated in iPSC-ECs from d-TGA patients, suggesting that both Tbx2 and GJA4 might be implicated in the development of d-TGA.
Regarding the downregulated genes expressed in EC from patients with d-TGA, we found several different genes involved in Notch signaling. One of these genes, JAG2, is directly implicated in angiogenesis [29], and DLL4 is required for normal arterial patterning [30]. In addition, DLL4 has an antiangiogenic effect when it competes with JAG1 [31]. Thus, different genes involved in Notch signaling could be directly involved in disease development. CDH5 was also found to be downregulated in iPSC-ECs from patients with d-TGA. CDH5 is a transmembrane cadherin protein required for vascular morphogenesis, as its deficiency impairs vasculogenesis in embryonic stem cells derived from embryoid bodies [32] and in mutant mice embryos [33]. These defects were observed, even though the endothelial cells were well differentiated [33]. The downregulation of KDR, which has an important role in endothelial specification and maintenance, was also observed. The activation of this receptor is critical for enhancing the proliferation and survival potential of iPSC-ECs [15].
At the functional level, we observed that endothelial cells derived from patients with d-TGA were unable to form capillary-like tubular networks. ECs are directly implicated in angiogenesis, which plays a crucial role during embryonic development and organogenesis [34]. These findings might be related to the low expression levels of genes involved in the Notch signaling pathway observed in the bioinformatics analysis.
The major limitation of our study is the lack of a multicellular system to simulate d-TGA in vitro. However, we have generated an in vitro model of d-TGA using iPSCs as a disease modeling approach. The gene expression differences between the ECs from control and d-TGA patients might be related to the etiology of this cardiac disease. The results presented in this work contribute to a better understanding of the mechanisms promoting the onset of d-TGA during embryonic development.
Overall, our results show that ECs from patients with d-TGA display genetic differences versus control counterparts. These changes affect biological processes related to the circulatory system and to vascular development. Accordingly, ECs derived from d-TGA iPSCs could be a good cellular model for the study of this pathology. This work supports the use of patient-specific iPSC-ECs in modeling d-TGA.

Human-Induced Pluripotent Cell Generation
Foreskin fibroblasts from a healthy individual and a patient with d-TGA (Ctrl1 and TGA1) were obtained through ATCC and the Coriell Institute, respectively. Fibroblasts were reprogramed following the method described by Yu in 2007 [35], using a lentiviral vector expressing the Oct4, Nanog, Lin 28 and Sox2 genes. Different clones were obtained and were validated using the following assays: alkaline phosphatase, surface marker expression, pluripotency gene expression, transgene silencing, proviral integration and in vitro and in vivo differentiation (data not shown).

Endothelial Differentiation
Endothelial cell differentiation was induced following the protocol described by Prasain et al. in 2014 [15]. Briefly, cells were seeded at a density of 6250 cm −2 onto a Matrigel pre-coated 24-well plate (day-2) and were incubated for 2 days in mTeSR1 medium. On day 0, cells were cultured with Stemline ® II Hematopoietic Stem Cell Expansion Medium (Sigma-Aldrich, Darmstadt, Germany) containing 10 ng/mL of Activin A (R&D Systems, Minneapolis, MN, USA), FGF-2 (Tebu-bio, Yvelines, France), BMP-4 (R&D Systems, Minneapolis, MN, USA) and VEGF (Peprotech, London, UK) for 1 day to initiate the differentiation. The following day, the medium was replaced and cells were incubated with the same concentration of FGF-2, BMP-4 and VEGF (differentiation medium) for 11 days. On day 12, FACS was performed to isolate endothelial cells.

Tube Formation Assay
iPSC-ECs were plated onto Matrigel-coated wells to measure tube formation, as described [38]. Briefly, 20 × 10 4 cells were seeded per well into 96-well plates and were incubated for 3 h. Images were acquired using an inverted microscope (Leica DM6000, Leica Microsystems, Wetzlar, Germany) with a 10× magnification and were analyzed with WimTube online software (WimTube: Tube Formation Assay Image Analysis. Release 4.0. https://www.wimasis.com/en/WimTube, accessed on 06 December 2021 ). Five images from independent experiments were analyzed for each iPSC-EC line.

RNA Extraction and Quantitative Real-Time PCR
RNA was isolated from samples on different days of EC differentiation (days 0, 5, 12 and 19) using the RNeasy Plus kit (Qiagen, Dusseldorf, Germany). cDNA was produced using a High-Capacity RNA-to-cDNA Kit (Applied Biosystems, Walthman, MA, USA), and real-time PCR (qPCR) was performed using the LightCycler 480 SYBR Green I Master Kit (Roche Life Science, Switzerland) on a ViiA 7 Real-Time PCR System (Applied Biosystems, Walthman, MA, USA). Primers were provided by Condalab (Madrid, Spain) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as housekeeping control. Reactions were performed in triplicate.

RNA-Seq
Libraries were prepared using the TruSeq Stranded Total RNA Library Prep Kit with Ribo-Zero Human/Mouse/Rat Kit (Illumina, San Diego, CA, USA). Briefly, 300 ng of total RNA was used for ribosomal RNA depletion. Then, ribosomal-depleted RNA was fragmented for 4.5 min at 94 • C. The remaining steps of the library preparation were followed according to the manufacturer's instructions.
Libraries were analyzed on an Agilent Technologies 2100 Bioanalyzer system, using the Agilent DNA 1000 chip to estimate the quantity and validate the size distribution, and were then quantified by qPCR using the KAPA Library Quantification Kit KK4835 (Roche Life Science, Switzerland).
FASTQ files were treated with Cutadapt and aligned with Kallisto to obtain the count table for each gene. Differential expression analysis of RNA-seq expression profiles was performed with the EdgeR package [39]. GO enrichment analyses were performed using g:Profiler [40,41].

Interactome Analysis
Analysis was performed as described [21]. Complete interactome was identified, and GSEA was performed to detect significant GO biological processes, molecular functions and cellular components [40].

Statistical Analysis
Results are represented as the mean ± standard deviation (SD). Comparisons between Ctrl and TGA groups were performed with the non-parametric Mann-Whitney U test, one-way ANOVA and necessary post hoc analysis. Analyses were conducted with Graph-Pad Prism 8 ® software (GraphPad Software Inc., La Jolla, CA, USA). Differences were considered statistically significant at p < 0.05, with a 95% confidence interval.