Characterisation of the developing heart in a pressure overloaded model utilising RNA sequencing to direct functional analysis

Abstract Cardiogenesis is influenced by both environmental and genetic factors, with blood flow playing a critical role in cardiac remodelling. Perturbation of any of these factors could lead to abnormal heart development and hence the formation of congenital heart defects. Although abnormal blood flow has been associated with a number of heart defects, the effects of abnormal pressure load on the developing heart gene expression profile have to date not clearly been defined. To determine the heart transcriptional response to haemodynamic alteration during development, outflow tract (OFT) banding was employed in the chick embryo at Hamburger and Hamilton stage (HH) 21. Stereological and expression studies, including the use of global expression analysis by RNA sequencing with an optimised procedure for effective globin depletion, were subsequently performed on HH29 OFT‐banded hearts and compared with sham control hearts, with further targeted expression investigations at HH35. The OFT‐banded hearts were found to have an abnormal morphology with a rounded appearance and left‐sided dilation in comparison with controls. Internal analysis showed they typically had a ventricular septal defect and reductions in the myocardial wall and trabeculae, with an increase in the lumen on the left side of the heart. There was also a significant reduction in apoptosis. The differentially expressed genes were found to be predominately involved in contraction, metabolism, apoptosis and neural development, suggesting a cardioprotective mechanism had been induced. Therefore, altered haemodynamics during development leads to left‐sided dilation and differential expression of genes that may be associated with stress and maintaining cardiac output.


Introduction
According to the British Heart Foundation, congenital heart defects are found in at least one in every 180 births in the UK (www.bhf.org.uk), with cardiomyopathies accounting for 8-11% of those diagnosed with a defect (Pedra et al. 2002). Both environmental and genetic factors affect embryonic heart formation, with blood flow playing a critical role in cardiac remodelling (Midgett & Rugonyi, 2014;Midgett et al. 2017;Courchaine et al. 2019). Animal models of haemodynamic alteration have been found to have phenotypes similar to those seen in human heart disorders, with features of cardiomyopathy (Sedmera et al. 1999;Midgett & Rugonyi, 2014;Midgett et al. 2017). The outflow tract banding (OFT-banding) haemodynamic alteration model has been found to give a wide spectrum of cardiac malformations, such as pharyngeal arch anomalies, ventricular septal defects (VSDs), valve defects and double outlet right ventricle accompanied by increased stress-strain relations and myocardial stiffness (Sedmera et al. 1999;Tomanek et al. 1999;Miller et al. 2003;Buffinton et al. 2013;Midgett et al. 2017;Pang et al. 2017). Further, OFT-banding in the chick leads to heart dilation (Clark et al. 1989;Tobita et al. 2002;Buffinton et al. 2013) and a reduction in compact myocardium (Tomanek et al. 1999). Myofiber alignment was also affected (Tobita et al. 2005). With regard to haemodynamic consequences of banding, chick embryos with banded hearts show an immediate and sustained increase of both peak systolic and end-diastolic ventricular pressure (Clark et al. 1989;Tobita et al. 2002), with an increase in blood flow velocity and wall shear stress (Chivukula et al. 2016). Despite a marked increase of ventricular pressure seen in OFT-banded embryos (Shi et al. 2013), heart rate and cardiac output/ejection fraction were not affected (Clark et al. 1989), suggesting preservation of ventricular function. Therefore, OFT-banding is a model for a heart with a complex congenital heart defect phenotype that is prone to abnormal haemodynamics and pressure overload, and hence heart enlargement (Sanchez-Gomez et al. 2017). Expression analysis in such hearts could pinpoint molecular pathways that may have protective or conversely debilitating long-term effects. Hence, the elucidation of expression mechanisms could help direct potential therapeutic studies or even aid the diagnosis of a heart under stress. However, such gene pathways are currently poorly defined.
In the study described here, the OFT-banded procedure was carried out at Hamburger-Hamilton stage (HH) 21, as the heart undergoes looping and chamber specification (Sedmera et al. 1999;Shi et al. 2013). Analysis was performed at key stages of development: at HH29 [Carnegie stage (CS) 19 for humans and embryonic day (E) 13.5 for mice], as chamber septation has just completed (Martinsen, 2005), and at HH35 (CS22 for humans and E15 for mice), when a fully septated heart with a mature apex to base conduction pattern is seen (Krishnan et al. 2014). Global gene expression analysis by RNA sequencing, together with stereological analysis, was performed on OFT-banded chick hearts and controls at HH29. The RNA sequencing was performed with an optimised procedure described herein for globin depletion in the chick. Furthermore, targeted expression analysis on genes identified by sequencing was performed by qPCR at HH35. Findings from these studies were then used to guide further functional studies regarding apoptosis and glycogen deposition. Expression analysis highlighted genes involved in energy regulation such as PRKAG3 (a regulatory subunit of the key metabolic regulatory enzyme AMPK) and LDHFB (involved in glucose metabolism) to be differentially regulated. Furthermore, genes such as S100A11 and PVALB (associated with improved contractility by increasing Ca 2+ sequestering during repolarisation) were upregulated. Less active AMPK and increased S100A11 expression can potentially lead to reduced apoptosis (Kanamori et al. 2004;Hwang et al. 2010;Law et al. 2017); a decrease in apoptosis was seen in the OFT-banded hearts. Altered AMPK regulation has the potential to increase glycogen storage; however, no increase in glycogen storage by the periodic acid-Schiff with diastase (PAS-D) assay was seen in OFT-banded hearts, suggesting a secondary mechanism. This study has identified left-sided dilation in the OFT-banded heart and showed gene expression that could provide a cardioprotective response to stress by attempting to maintain energy metabolism needs and contractility, together with a decrease in apoptosis.

Outflow tract banding, embryo isolation and histological analysis
White fertile chicken eggs (Gallus gallus; Dekalb white strain; Henry Stewart, UK) were incubated at 38°C in a humidified atmosphere under constant rotation for 4 days until HH21 (Hamburger & Hamilton, 1992). After windowing, the inner shell membrane was removed to expose the heart. OFT-banding was performed using 10-0 nylon suture to create a snug double knot around the OFT, which was removed upon harvesting. The suture was passed through and removed immediately on sham controls. Untreated (UT) control embryos were opened and staged but no further procedures performed (Sedmera et al. 1999;Pang et al. 2017). All eggs were sealed and incubated for an additional 3-5 days until HH29-35. Animal work was performed in accordance with national (UK Home Office) and institutional regulations and ethical guidelines. OFT-banded and control (sham and UT) embryos were isolated and externally analysed. Hearts were then fixed in 4% paraformaldehyde, processed and wax-embedded in a transverse orientation. Serial 8-lm sections were taken (DSC1 microtome, Leica, Germany), dewaxed and rehydrated. For histological studies, sections were stained with Alcian blue (Sigma, UK) for 15 min at room temperature (RT) followed by Mayer's haemalum (Raymond Lamb, UK). Images were acquired with a slide scanner (Nanozoomer 2.0-HT, Hamamatsu, Japan).

Stereology and morphometric measurement
Systematic random sampling (Mayhew, 1991) was used to assess tissue proportions throughout HH29 hearts (n = 6 sham, n = 7 OFTbanded). A 96-point grid was placed over every fifth section throughout the heart, and the tissue region and type on each point was identified (4479 points for sham, 7868 points for banded). Tissue regions consisted of right (RA) and left atrium (LA), right (RV) and left ventricle (LV). In addition, the tissue types counted included myocardial wall, extracellular matrix and lumen in the regions of RV and LV, and myocardial wall and lumen for RA and LA. Average tissue proportions of each group were calculated and tested for statistical significance.

Apoptosis
ApopTag Peroxidase In Situ Apoptosis Detection Kit S7100 (Millipore, USA) was used to indicate apoptotic cells in accordance with manufacturer's instructions on 5-lM serial sections. Imaging was performed using Zeiss Axio Scan Z1. Systematic random sampling (Mayhew, 1991) was utilised to count positive cells against total cell count on the left and right ventricular region of the heart, including the ventricular compact myocardium and trabeculae, and the base and the myocardial crest of the IVS to calculate proportions of apoptotic cells for statistical analysis. A total of 303 299 cells were counted. manufacturers' instructions. Quantitative analysis was performed in IMAGEJ using colour deconvoluted images (Ruifork & Johnston 2001). Images were taken using an AxioPlan (Zeiss). Statistics were performed using a two-way analysis of variance (ANOVA).
RNA isolation and cDNA synthesis for RNA sequencing and qPCR RNA isolation on all hearts within direct comparison experiments were performed simultaneously by the same handler. Total RNA was extracted using TRIzol reagent (Sigma) and treated with RNasefree DNase I (Qiagen). RNA purity was checked by NanoDrop 2000c UV/IV spectrophotometer (Thermo Fisher Scientific) at 260 nm absorbance for RNA. A260/280 ratios were 1.90-2.20 and A260/230 ratios were 2.0-2.20. R.I.N values of 9 were consistently seen following RNA purification. The values of 6.9-7.4 used for sequencing was due to sample degradation following multiple freeze thaws and globin depletion treatment. A single-stranded template was used for standard qPCR. cDNA synthesis for RNA sequencing was doublestranded, purified and checked for concentration via Qubit assay (Thermo Fisher Scientific). Reverse transcription reactions were performed using SuperScriptâ II Reverse Transcriptase (Invitrogen, UK) following the manufacturer's instructions. All cDNA synthesis was carried out with 1 µg purified RNA in a 20-µL reaction. cDNA from hearts within the same experiment were always synthesised in the same run, with no reverse transcriptase controls for assessment of gDNA. Hexamer primers were used at 1.25 µM. Following the reaction, samples were instantly placed on ice and then stored at À20°C.

Quantitative PCR
Quantitative PCR (qPCR) was run using the Applied Biosystems 7500 Fast Real-time PCR system using SYBR Select master mix (Bio-Rad). Gene-specific primers were designed to give fragments of 90-210 bp and 62°C T m (AE 1°C) (Supporting Information  Table S2). Primers were designed to exon/exon boundaries or to include introns > 1000 bp, except for ENS-1/ERNI, which contained no introns. For these genes, no template controls showed minimal gDNA expression (≥ 10 Cq later than cDNA samples with Cq values > 38). Primers were optimised with cDNA diluted four times, leading to a dilution series using six points at a ratio of 1 : 5 or 1 : 3, depending on the level of gene expression. Final template dilutions for relative expression analysis were elucidated from a midway point on the linear portion of the range for each gene. All primer standard curve R 2 values were > 0.98 with efficiencies of 96-108%, with reference genes GAPDH, EEF1A1 and TBP having efficiencies of 101-105%. A single heat-denaturing step at 95°C for 25 s, followed by 40 cycles of denaturation (95°C for 20 s), annealing and extension (62°C for 1 min) was used. Each 20-µL qPCR reaction mixture consisted of 10 µL of iTaq TM universal SYBRâ Green Supermix (19), 0.5-0.75 µL (250-375 nM, respectively) primer combinations and optimised template concentrations. All samples were run in triplicate. Following the qPCR reaction, a melt curve was performed. Reference genes were shown to be unaffected by the OFT-banded experimental conditions. Two reference genes were used for analysis: GAPDH and EFF1A1 at HH29 and GAPDH and TBP at HH35. The threshold cycle number for product detection (DT value) was used to calculate the relative gene expression against reference genes using the Pfafll adjusted efficiency 2 ÀΔΔCt method (Pfaffl, 2001).
Statistics were carried out on normalized change in Cq values between OFT-banded and sham.

Globin depletion oligonucleotide design
Throughout all globin depletion/library preparation procedures, all OFT-banded and sham preparations were performed simultaneously by the same handler. Globin depletion (GD) was performed via enzymatic methods on hybridised RNA/DNA complexes. Before the procedure could be carried out, specific oligonucleotides had to be designed for each globin gene. qPCR confirmed high relative expression of HBAA, HBAD, HBB, HBG1, HBG2 and HBZ genes and oligonucleotides were designed using PRIMER 3 (v. 4.1.0) and BLAST against these globin genes for depletion (Supporting Information  Table S1). Oligonucleotides were designed to the 3 0 end and had a T m 60-70°C.

Globin depletion procedure
The chick-specific designed oligonucleotides were used for GD using a modified Affymetrix protocol (Wu, 2007). A final concentration of 0.75 µM for each oligonucleotide was used for hybridisation to globin transcripts with 3 µg of total RNA. The hybridisation procedure time and temperature were optimised and a hybridisation protocol of 95°C for 2 min, cooled down at 5°C intervals in 50-s steps to 50°C, with longer 1.5-min steps at 65°C and 60°C (the closest temperatures to oligonucleotide T m ) was devised that allowed all oligos to anneal to their targets. Reaction buffer consisted of 500 mM Tris-HCL, 1 M NaCl 59 stock with 19 used in the final hybridisation solution made up to 15 µL with DNase-free water as necessary. Following the 50°C step, samples were placed immediately on ice to help reduce secondary structure formation. Hybridised RNA/DNA was degraded by 4 U RNase H (0.2 U µL À1 in the total 20-µL hybridisation solution) at 37°C in RNase H buffer (NEB) containing SUPERase-In. A 2-U aliquot of DNase (0.07 U µL À1 in 30 µL) was then added to the mix after 10 min, and following a further 10 min, samples were placed on ice and the reaction was instantly stopped with ethylenediaminetetraacetic acid (EDTA) to a final concentration of 30 mM. The additional DNase step helps to degrade the DNA hybridised oligonucleotides and any genomic DNA contamination. Samples were then purified using Purelink RNA columns to allow GD total RNA to be quantified by Qubit Fluorometric assay. Statistical analysis of pre-and post-GD treatment was performed by paired Student's t-test.
Library preparation and RNA sequencing RNA concentration was quantified using the Qubit RNA HS Assay Kit (Thermo Fisher Scientific). RNA integrity was assessed using Agilent 2100 Bioanalyzer (Agilent, USA) with an RNA integrity number (RIN) of 6.9-7.3 used across OFT-banded and sham hearts for sequencing. Following GD, 1 µg RNA from each heart (OFT-banded n = 3, sham n = 3) was used for RNA sequencing library preparation according to manufacturer's instructions of the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs Inc.; NEB UK #E7420S/L). Twelve cycles were used for PCR enrichment. The barcoded strand-specific libraries created were sequenced on the Illumina NextSeq 500 platform. Eight-tenths of a full high output sequencing run was performed according to the manufacturer's instructions. Forty million total reads came from two flow cell reads with each read having 75 bp paired-end reads.

RNA sequencing bioinformatic analysis
Following sequencing, compressed FASTQ files were trimmed to filter reads for adapter and low-quality bases. SYTHE was used to remove adaptor sequences (https://github.com/vsbuffalo/scythe), with reads of low sequencing scores then trimmed using SICKLE. Reads were then aligned to the galGal4 reference genome (ICGSC Gallus_gallus-4.0/galGal4) in the UCSC database Nov.2011 (http://ge nome.ucsc.edu) in the context of known gene exon coordinates by TopHat (https://ccb.jhu.edu/software/tophat/index.shtml). Alignment counts were recorded in BAM format with only uniquely mapped read alignment hits (47-90 9 106) having a quality score of MAPQ ≥ 20 used. Htseq count was used for generation of count tables. Read counts were then normalised using FPKM and TPM with resultant figures from both used for PCA using R software. As both normalisation techniques revealed nearly identical PCA plots, FPKM values were used for differential gene regulation analysis by DESeq with corrected P-values of < 0.05 set as the DESeq threshold (Anders & Huber, 2010). Genes with average FPKM values < 1 were filtered from the FPKM data before gene regulation was assessed using DESeq adjusted P-values (FDRs).

Statistical analysis
All data are expressed as AE SEM, but a confidence interval of 95% was also checked. Statistical significance of the differences between OFT-banded and sham hearts was analysed by independent t-test assuming equal or unequal variances; *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. Assumption of normality was tested using a Shapiro-Wilk normality test. For assumption of variance, Levene's test for homogeneity of variance was used. Tests were carried out using R, SPSS and Microsoft EXCEL statistical analysis software.

Stereology reveals enlarged OFT-banded hearts with changes in the chamber components representing left-sided dilation
To define the morphological phenotype of the OFTbanded model, stereological analysis was performed at HH29, confirming a phenotype of left-sided dilation (Fig. 2). Whole OFT-banded hearts were found to be 1.77 times larger than sham controls. A significant reduction of the trabeculae was identified within both ventricular chambers of the OFT-banded hearts. In the right ventricle (RV) a 27.3% decrease was seen (10.66 AE 0.36% sham, 7.75 AE 0.55% OFT-banded; P < 0.001). Decreases were more prominent in the left ventricle (LV) with a 35.9% reduction seen (13.30 AE 1.65% sham, 8.53 AE 0.37% OFTbanded; P < 0.001). Consistent with this, a simultaneous increase of LV lumen was found in OFT-banded hearts (9.94 AE 1.06% and 14.27 AE 1.40% in sham and OFTbanded hearts respectively, an increase of 43.6%; P < 0.05; Fig. 2). This coincided with a thinner myocardium in the LV by 40.7% upon banding (9.41 AE 1.05% sham, 5.58 AE 0.43% OFT-banded; P < 0.01); however, no difference was seen in the RV (P > 0.05; Fig. 2). With regard to the atrial myocardium, a reduction of 37.07% was identified in the myocardium of the left atrium (LA) in OFTbanded hearts (4.99 AE 0.67% sham, 3.14 AE 0.24% OFTbanded; P < 0.01). However, decreases in the right atrium (RA) were found to be insignificant (5.16 AE 0.94% sham, 3.24 AE 0.34% OFT-banded; P > 0.05; Fig. 2). A significant increase of the LA lumen by 47.4% was found in OFTbanded hearts (4.9 AE 0.45% sham, 7.22 AE 0.62% OFTbanded P < 0.01; Fig. 2). In contrast, no significant differences were seen in the lumen of the RA (6.04 AE 0.76% sham, 5.68 AE 0.74% OFT-banded; P < 0.05).

Globin depletion in the chick embryo
Following the identification of morphological phenotypes at HH29, the effect of OFT-banding on gene expression was assessed by global RNA sequencing of whole HH29 hearts. To enhance discovery of possible transcripts influenced by OFT-banding, particularly those of low reads, a library preparation procedure was adapted (Choi et al. 2014). RNase H was used to remove the polyA tail of globin transcripts prior to polyA selection. Expression of globin genes HBAA, HBAD, HBB, HBG1, HBG2 and HBZ in the chick embryo was confirmed by qPCR. Then, to enable their removal, specific oligonucleotides were designed for each globin gene at their 3 0 end (oligonucleotide sequences in Table S1). The RIN was 7.1-7.9 before globin depletion (GD) treatment for all samples, and 6.9-7.3 following GD.
Globin transcripts in GD-treated hearts showed a significant reduction of À11.10 to À17.49 log 2 fold change (FC) by qPCR. This represents values of 0.045-0.0005% that of the original amount in the same hearts before treatment (t = 16.003, df = 5, P < 0.0001; Fig. 3A). Variation in gene expression comes from the variation between the biological replicates, as qPCR readouts show GD-treated quantitation cycles (Cq) were consistent (see Supporting Information Table S3). This demonstrates that GD treatment removed globin genes consistently to a low point but did not remove all transcripts. It also shows that using the specific GD oligonucleotide concentration of 0.75 µM, the initial expression level of the globin gene was not a limiting factor for the number of transcripts remaining post-GD treatment.

Sequencing of chick globin-depleted RNA is highly reproducible
Globin depletion has been shown to lead to reduced RNA concentrations of 33-95% to that of pre-treatment values (Choi et al. 2014). To assess whether gene expression results were affected by the GD procedure, technical verification was required post-RNA sequencing. This necessitated the use of qPCR to confirm that gene regulation identified by RNA sequencing in GD hearts was also seen in the same hearts without GD. A profile of 17 genes with an expression range of À1.35 to 0.81 log 2 FC (OFT-banded n = 3, sham n = 3) was used. No significant differences in expression were seen between GD sequenced hearts and the same hearts prior to treatment by qPCR, thus confirming sequencing expression (t = 0.145, df = 32, P = 0.78; Fig. 3B). Harvesting, RNA extraction, and library preparation were performed simultaneously across all OFT-banded and sham samples to mitigate any batch effects. Housekeeping and qPCR reference genes GAPDH, EEF1A1 and TBP showed consistent expression across all sham and OFT-banded hearts by RNA sequencing at HH29. This demonstrated that library preparation was consistent across samples and confirmed the Genorm, BestKeeper and normFinder data, demonstrating stable expression in HH29 biological repeats and HH35 OFT-banded hearts by qPCR (Fig. S3, Table S3; Supporting Information Fig. S4 also shows efficiency calibration and electrophoresis traces). The consistency of gene expression following GD was further confirmed in studies where independent untreated hearts pre-and post-GD treatment (HH29 hearts, n = 3 per group) were used. Gene FC values remained unchanged when normalised to GAPDH and EEF1A1 by qPCR (Supporting Information Figs S1 and S2). The high consistency across gene regulation shows that the GD library preparation procedure in the chick does not lead to altered gene regulation, and also that the technical repeatability of the RNA sequencing result. The RNA sequencing and qPCR results are highly correlated (r = 0.97, df = 16, P = < 0.01; Fig. 3C) with a high level of sensitivity between small differences in fold change.
OFT-banded heart model shows biological repeatability of gene regulation Biological repeatability of the OFT-banding model was also determined at HH29. Gene regulation in RNA used for sequencing (n = 3 OFT-banded, n = 3 sham) was compared with RNA of newly harvested hearts using qPCR (n = 3 OFTbanded, n = 3 sham). Results showed similar gene regulation log 2 FC results that were significantly correlated (r = 0.91, df = 5, P = < 0.01; Fig. 3D). Expression of selected genes ranged from À1.35 to 0.81 log 2 FC. These data demonstrate the consistency of gene expression within the OFT-banded hearts.

Differential gene expression in OFT-banded hearts
Gene regulation in OFT-banded hearts was assessed by RNA sequencing using HH29 OFT-banded (n = 3) and sham (n = 3) hearts. In all, 40 million paired-end reads were performed on these samples to create a global expression profile. Reads were matched to the galGal4 reference genome, identifying 6920 genes. Principal component analysis (PCA) revealed within-group hearts were clustered together, showing consistency of expression, but that between groups the OFT-banded and sham hearts were distinct (Fig. 4A). PCA was performed with read values normalised to total fragments (FPKM), and for comparison with total transcripts (TPM) with genes with expression levels < 1 FPKM/TPM discounted. However, no proportional differences were found in the PCA between the two normalisation methods. K-means clustering was also performed on FPKM values on a subset of 45 genes linked to cardiac development and stress. This revealed that, without bias, all OFT-banded hearts mapped a distinct centre from that of all sham hearts (Supporting Information Fig. S6).
Differential gene expression was analysed using DESeq with 28 genes being identified as differentially expressed in OFT-banded hearts (using a P-value cut-off of 0.05; Table 1). To ensure differentially regulated genes of biological interest were not missed, significant genes had no minimum fold change cut-off; however, qPCR was used to confirm differential regulation in all genes analysed further. Genes of interest highlighted as significant had corrected fold changes of 1.38-1.87 by RNA-seq and fold changes of 1.46-2.56 by qPCR, confirming their differential regulation (Fig. 3B,C, Supporting Information Table S5), with all genes showing expression levels > 3.5 FPKM before any differential regulation.
DESeq uses corrected false discovery rates (FDR) P-values. Genes identified by this method have been shown to be Fig. 2 Stereological analysis of tissue types contributing to different heart regions upon OFT-banding. Reduction of trabeculae was found in both right and left ventricle of the OFT-banded hearts, with an increase in left ventricular lumen. In addition, the myocardium was thinner in the left ventricle but normal in the right ventricle. Also, banding led to an increase of the lumen and a decrease of myocardium in the left atrium. ECM, extracellular matrix; Myo, myocardium. Sham, n = 6; banded, n = 7. Significant differences are indicated: *P < 0.05; **P < 0.01; ***P < 0.001. Error bars indicate SEM.
highly verified (Anders & Huber, 2010). However, DESeq is a conservative estimator (Soneson & Delorenzi, 2013), particularly with small experiment sizes. To test whether sequencing results could be used to drive further research in genes of near significance, two genes, PRKAG3 and HDAC7, with P-values of 0.06 and 0.07, respectively, were further analysed by qPCR. PRKAG3 returned a significant P ≤ 0.05 by qPCR; however, HDAC7 remained insignificant. The qPCR results show that the sequencing data (full gene regulation results provided in Supporting Information Table S8, and the sequencing data are accessible through GEO Series accession number GSE120498) can be used to target genes for further investigation and that genes just outside of significance could still potentially be targets.

Targeted gene expression study at HH29 and HH35 highlights potential biological mechanisms
To elucidate the effect of OFT-banding beyond HH29, observation of differential expression in genes of interest at HH29 led to their further study at HH35 by qPCR. Genes were selected based on their possible cardiac biological Values are AE1 SEM (n = 6; ****P < 0.0001). (B) Gene expression of globin-depleted and sequenced hearts, compared with the same hearts nonglobin-depleted by qPCR. The lack of statistical differences shows that globin depletion/sequencing preparation does not effect gene fold change (FC) values: t(32) = 0.145, P = 0.78 (n = 17). RNA sequencing non-corrected log 2 FC was used for direct comparison with qPCR log 2 FC. (C) Relationship of corrected (DESeq) RNA sequencing to qPCR. Results are highly correlated and show that qPCR generally gives a small increase in log 2 FC when compared with RNA sequencing: r(16) = 0.97, P ≤ 0.0001) (n = 17). RNA sequencing DESeq corrected log 2 FC compared with qPCR log 2 FC. A closer parallel can be seen between non-corrected RNA sequencing and qPCR log 2 FC as seen in Supporting Information Fig. S5. (D) Relationship of gene expression of RNA-sequenced hearts (n = 3 OFT-banded, n = 3 sham) and new biological replicates (n = 3 OFT-banded, n = 3 sham) at HH29 by qPCR. Statistically significant correlation of biologically independent samples is shown: r(5) = 0.91 (P ≤ 0.01).
significance together with RNA sequencing pathway data (Table S6). Genes had potential roles in metabolism, contraction, apoptosis and neural development. PRKAG3, a gamma regulatory unit of key metabolic regulator AMPK showed a downregulation at HH29 in OFT-banded hearts (À0.68 log 2 FC); this differential regulation increased by 104% at HH35 (À1.18 log 2 FC; Fig. 4B, Table S5). Sequencing also revealed that there is no compensatory response from PRKAG3 isoforms at HH29, with PRKAG2 expressed at very low levels of 0.1-0.4 FPKM across all samples (Supporting Information Table S7). LDHB promotes oxygen metabolism in the heart and at HH29 shows significant upregulation (0.62 log 2 FC); however, differential regulation is eliminated at HH35 (Fig. 4B). Genes which play roles in calcium sequestering and hence are potentially important in contraction, namely S100A11, PVALB, MYL1 and TNNC2, all showed mRNA upregulation at HH29 (Fig. 4B). At HH35, intercellular membrane protein S100A11, cytoplasmic protein PVALB and structural protein MYL1 showed further increases in gene regulation of 43%, 60% and 114%, respectively; however, the structural protein TNNC2 showed a slight reduction in expression (Fig. 4). Neural developmental genes ERNI and ENS-1 show a significant downregulation of À1.08 and À1.35 log 2 FC, respectively, at HH29. ERNI differential expression was maintained at HH35; however, the downregulation of ENS-1 was reduced to insignificant levels (Fig. 4B).

OFT-banded hearts show no changes in glycogen deposition
Expression data showing downregulation of AMPK regulator PRKAG3 could be indicative of altered glycogen storage, with glycogen storage disease and PRKAG2 (isoform of PRKAG3) mutation seen in some cases of cardiomyopathy (Porto et al. 2016;Banankhah et al. 2018). To see whether PCA of RNA sequencing FPKM shows consistency of expression within OFT-banded (OFT) and sham groups, with a distinction seen between groups. (B) Expression analysis of targeted genes at HH29 and HH35 by qPCR shows a general increase in differential expression in OFTbanded hearts, with the exception of LDHB and ENS-1, where decreases are seen (*P < 0.05; significantly differentially regulated in OFT-banded hearts by qPCR). HH29 data are from hearts used for RNA sequencing compared with data from isolated HH35 hearts (n = 3 OFT-banded, n = 3 sham).
any phenotype of glycogen storage was present, glycogen deposition was analysed by periodic acid-Schiff stain with diastase (PAS-D) at HH35 in OFT-banded and sham hearts (n = 3 per group). HH35 was chosen because this stage represented the highest level of differential regulation of PRAKG3. Glycogen storage disease was not present; no differences were seen in glycogen deposition between OFTbanded and control hearts (Fig. 6).

Discussion
Environmental and genetic factors affect embryonic heart formation, with blood flow playing a critical role in cardiac remodelling (Midgett & Rugonyi, 2014;Courchaine et al. 2019). Though many haemodynamically altered heart models exist in the literature (Sedmera et al. 1999;Tomanek et al. 1999;Miller et al. 2003;Buffinton et al. 2013;Chivukula et al. 2016;Midgett et al. 2017), a global molecular characterisation of whole hearts during cardiogenesis has not been performed. Additionally, RNA sequencing was carried out on library preparations with globin transcripts removed. The globin depletion technique has been shown to increase the detection of non-globin transcripts in expression studies, allowing greater detection of genes with low read counts (Wu, 2007;Choi et al. 2014). However, it has not been performed previously in the chick embryo. Not all species have high levels of globin mRNA, with a recent study showing that equine and bovine blood have low levels (Correia et al. 2018). However, the data presented here show six globin genes with very high expression relative to reference genes such as GAPDH, suggesting the need for globin depletion in the chick.
The OFT-banded heart is a typical pressure overloaded embryonic heart model that gives rise to a range of structural malformations such as pharyngeal arch anomalies, persistent truncus arteriosus, ventricular septal defects and valve defects (Clark & Rosenquist, 1978;Sedmera et al. 1999;Tobita et al. 2002;Hall et al. 2004;Midgett et al. 2017;Pang et al. 2017), while also showing a left-sided dilation. In addition, previous data from our group show changes in the ECM of the epicardium at a molecular and morphological level with a progressive aberrant phenotype from HH29 to HH35 (Perdios et al. 2019). However, overall embryonic growth was not affected upon banding, as the cardiac output and overall embryo size are the similar in control and banded embryos (Clark et al. 1989;Keller et al. 1997;Pang et al. 2017). In this study, evidence is provided of molecular mechanisms involved in OFT-banded hearts, highlighting potential cardioprotective responses, with further analyses characterising the levels of apoptosis and glycogen accumulation. Though some dilation of the ventricles has been reported previously (Buffinton et al. 2013), this study found dilation of both the left ventricle and atria, with an increased lumen and a decrease in myocardium and trabeculae. Though OFT-banded models have been shown generally to give rise to a thicker myocardium (Clark et al. 1989;Sedmera et al. 1999;Tobita et al. 2002), analysis at later stages (HH36) has shown a reversal of this, with a significant reduction seen (Tomanek et al. 1999). Left atrial dilation has been postulated as a consequence of more extensive left ventricular damage, and dilated cardiomyopathy patients with dilation of the ventricle accompanied by atrial dilation are associated with a higher risk of morbidity (Modena et al. 1997). Therefore, the left-sided dilation and reduced myocardium seen herein could be indicative of a more advanced OFTbanded phenotype (Sedmera et al. 1999;Tobita et al. 2002). Heart dilation has been seen in humans subsequent to a previous diagnosis of a hypertrophic heart, suggesting the heart is failing (Hamada et al. 2010). Conversely, the developing or postnatal heart with a complex congenital defect (having two or more heart defects) is prone to abnormal haemodynamics, with a dilated heart seen as a secondary effect in some cases (Sanchez-Gomez et al. 2017;Courchaine et al. 2018). In addition, a recent study has indicated that as heart development progresses (HH22 to HH36), the proximal part of the OFT (conus) becomes part of the right myocardial wall (Lazzarini et al. 2018); this process could rescue the right ventricular wall from a dilating phenotype. Therefore, to assess the molecular mechanisms involved with this phenotype during cardiogenesis, we performed global RNA sequencing on OFT-banded and control hearts.
Sequencing found differential regulation of genes predicted to be involved in heart metabolism, apoptosis, neural development and contractility (see pathway analysis Table S6). qPCR confirmed that the most differentially expressed genes identified can be used as a guide for further functional studies. This differential gene expression as a consequence of the OFT-banding procedure could be expected as a response to stress in order to maintain cardiac output.
One of the more unexpected sequencing results was that two of the three top differentially regulated genes were the neural developmental genes ENS-1 and ERNI. These genes are expressed during proliferation and are downregulated upon differentiation (Jean et al. 2013). They are known to be expressed in pluripotent cells of the epiblast and later in the prospective neural plate (Mey et al. 2012). Another gene, Gbx2, known to be involved in cardiac neural crest activation (Calmont et al. 2009;Simoes-Costa & Bronner, 2015) showed an upregulation in OFT-banded hearts. Multipotent cells derived from the neural crest have been detected in the conduction system of the heart into late development (Keyte & Hutson, 2012), and ablation of neural crest cells has led to a retardation of the apex to base conduction pattern (Gurjarpadhye et al. 2007). This would suggest that these cells have a role in cell differentiation in the heart. The downregulation of ENS-1/ERNI, along with the upregulation of Gbx2, suggests that OFT-banded hearts could suffer from abnormalities caused by aberrant neural crest differentiation. Further investigations to analyse neural crest differentiation would be of interest in future studies.
Expression of PRKAG3, a gamma regulatory unit of key metabolic regulator AMPK, was confirmed in the embryonic chick heart at HH29 and HH35. PRKAG3 and PRKAG1 are the dominant isoforms, with PRKAG2 showing near-zero levels of expression at HH29. In the normal human adult heart, PRKAG2 is highly expressed, whereas PRKAG3 shows very low expression (Pinter et al. 2012). Gamma subunits are co-ordinately regulated in response to metabolic requirements (Pinter et al. 2013), and it is likely that in the embryonic chick heart, PRKAG3 is expressed to compensate for the low levels of PRKAG2. Further, PRKAG2 has been associated with heart disorders such as cardiac hypertrophy, left ventricular dilation, supraventricular tachyarrhythmias, ventricular pre-excitation and glycogen storage cardiomyopathies (Porto et al. 2016). The gamma subunit in response to AMP/ADP alters AMPK conformation, promoting phosphorylation and subsequently stopping dephosphorylation and AMPK inactivation (Zaha & Young, 2012). Thus, the reduction in PRKAG3 upon OFT-banding would be expected to lead to reduced sensitivity of AMPK to AMP/ADP/ATP levels. AMPK is a cell master metabolic regulator with many functions involved in protein kinase and transcriptional regulatory activity (Zaha & Young, 2012;Bairwa et al. 2016). Therefore, affecting AMPK activity could have many knock-on effects in OFT-banded hearts and is an area for further study.  6 Glycogen storage is not affected in HH35-banded hearts. Sham and OFT-banded heart sections were treated with periodic acid-Schiff (PAS). PAS without the addition of diastase (À) (a,b) or with a diastase reaction taking place (+) (c,d). There was no apparent difference in the treatment groups or between the right (a',b') and left side (a'',b'') of the hearts in non-diastase-treated and diastase-treated right (c',d') and left sides (c'',d''). LA, left atrium; LV, left ventricle; RA, right atrium; RV, right ventricle. Scale bars: 1000 lm (a-d), 100 lm (a',a' ',b',b'',c',c'',d',d''). n = 3 for each group.
Differential expression was also found for another metabolic gene that plays a key role in glucose metabolism, lactate dehydrogenase B (LDHB). LDHB is the predominant isoform in the adult heart (Valvona et al. 2016), and our sequencing confirms this is also the case in the chick embryo. It favours the reaction of lactate to pyruvate, thus promoting the TCA cycle and oxidative metabolism (Schisler et al. 2015). Under conditions of pressure overload and anaerobic stress, the heart can alter metabolism dependence in order to maintain output in response to substrate and/or oxygen availability and load (Das et al. 1987;An & Rodrigues, 2006;Lai et al. 2014). The data described here show significant upregulation at HH29; however, differential regulation is eliminated at HH35. This would suggest that at HH29 the OFT-banded heart is trying to promote oxidative metabolism in an effort to maintain cardiac output. However, by HH35, metabolism may be becoming increasingly anaerobic, perhaps as oxygen levels are not high enough to maintain output, and the heart switches again to try to preserve cardiac output. This study analysed hearts during a period of active remodelling during cardiogenesis; metabolism is known to be involved in heart remodelling (Doenst et al. 2013). Therefore, these results would suggest a disruption to energy metabolism in OFTbanded hearts that could have a bearing on heart morphology.
The PAS-D assay was performed to examine any differences in glycogen storage between the OFT-banded and control hearts. Glycogen stains a purple-magenta colour using PAS and the amount of glycogen can be measured semi-quantitatively by the intensity of the stain (McManus, 1948). The addition of diastase acts as a control for the diastase negative sections, as glycogen is not the only PAS-positive element. PAS-D can also identify other PAS-positive elements that can be a sign of a disorder (e.g. increase in mucins) because the glycogen is digested and washed out (McManus, 1948). Although the data described here found differential regulation in genes associated with glycogen storage, the PAS-D assay showed that the glycogen storage of the experimental model was not significantly affected in comparison with shams, suggesting glycogen storage disease was not present.
This study has also shown significant upregulation in multiple Ca 2+ sequestering genes, which play roles in the heart's ability to contract. Such increases could be expected as part of the heart's response to stress in an attempt to increase contractility. A recent study has shown that treatment with a drug expected to increase contractility (phosphatase-inhibitor-1) in a pressure overload model actually led to heart deformities, but non-pressure overload hearts showed the expected improvement with treatment (Schwab et al. 2018). The differential expression of PVALB and S100A11 are of particular interest. PVALB expression was increased upon OFT-banding in HH29 hearts, with differential expression seen to increase at later stages. PVALB, expressed in fast-conducting skeletal muscle, has been used in preclinical rat trials where its expression accelerated relaxation by the action of sequestering calcium from TNNC2 to the sarcoplasmic reticulum (Szatkowski et al. 2001;Wolfram & Donahue, 2013). Its action is independent of ATP (Wolfram & Donahue, 2013), and with our results showing disruption to genes involved in energy metabolism, a protein that is not energy dependent could be an ideal therapeutic target. S100A11 has been found to be co-expressed with S100A1 on internal cell membranes and at Z-discs in skeletal muscle (Arcuri et al. 2002). S100A11 has not been linked directly with pressure overload (Thuny et al. 2012) but has been found to be upregulated following b-adrenergic stimulation in the rat heart (Inamoto et al. 2000). S100A1-deficient mice have been shown to have impaired cardiac contractility in response to pressure overload (Inamoto et al. 2000). Our sequencing at HH29 reveals S100A11 to be dominantly expressed over S100A1, and hence a potential cardioprotective effect seen to be increasing from HH29 to HH35. Overexpression of S100A11 is associated with reduced apoptosis, whereas suppression has been seen to lead to apoptosis (Kanamori et al. 2004). The decrease in apoptosis seen in OFT-banded hearts suggests the increase in S100A11 expression may be involved.
Upregulation of PDCD4 is associated with apoptosis and is a target of microRNA-21 and microRNA-499-5p, both of which have been linked to cardioprotection (Cheng et al. 2010;Liu et al. 2012;Jia et al. 2016;Li et al. 2016). Therefore, downregulation of PDCD4 expression with a concomitant decrease of apoptosis in this study suggests the embryonic hearts were undergoing an intrinsic cardioprotective mechanism in response to pressure overloading. Further, AMPK has been found to be a key regulator of cell death in response to stress, with AMPK activation coinciding with cell cycle arrest and induction of apoptosis (Hwang et al. 2010;Law et al. 2017). Apoptosis is not a characteristic of the normal human adult heart; however, increased apoptosis has been linked to a heart that is failing (Narula et al. 2000;Zhao et al. 2016). A decrease in apoptosis has been shown to provide a protective/defensive mechanism in the fetal heart (Schaffer et al. 2000;Tao et al. 2015). Therefore, these studies support the notion that the heart is undergoing a protective response, and are further evidence that the heart is not yet failing.
With a key period of chamber development being haemodynamically altered, from banding at HH21 to RNA sequencing and structural analysis at HH29, it could be that the differential gene regulation seen in metabolic genes, though attempting to respond to changes in load, had a direct bearing on the anatomical phenotype seen in OFTbanded hearts. The pressure overloaded developing chick heart was found to have left-sided dilation and a decrease in apoptosis. Upon analysis of the global gene expression profile using RNA sequencing on globin -depleted RNA, changes in gene expression related to metabolism, apoptosis, neural development and contractility, with further targeted expression studies at a later stage of development performed. These data suggest the heart is responding in a cardioprotective manner. Together, this study provides insights into the effect that altered haemodynamics has on heart structure and function, which will help increase the understanding of heart development while subjected to pressure overload. This may focus future functional and therapeutic studies to elucidate the mechanisms involved.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Fig. S1. Relative mRNA quantity comparison in RNase-treated (+RNase, globin-depleted) and non-RNase treated (-RNase, nonglobin-treated) samples. An approximately one-fold decrease in relative quantity was seen in globin-depleted samples. Three untreated control HH29 hearts where used, with the same hearts used in both groups. Fig. S2. Normalised expression of GJA5 in non-experimental control hearts for RNase-treated (globin-depleted) and non-RNasetreated (non-globin-depleted). Fig. S3. Analysis of reference gene stability by qPCR in HH35 OFT-banded hearts by Genorm, BestKeeper and normFinder (n = 4 OFT-banded, 4 sham, 4 untreated hearts). Fig. S4. qPCR reference gene primer optimisation. Fig. S5. RNA sequencing correlation in OFT-banded and sham controls, using non-corrected fold change (FC) values compared with qPCR, demonstrate a close relationship: r(16) = 0.98, P ≤ 0.0001 (n = 17). Fig. S6. OFT-banded hearts show unbiased mapping to distinct clusters based on a profile of 45 genes linked to cardiac development and stress (all genes are shown in Tables 1 and S7)  Table S1. Oligonucleotide design for globin gene hybridisation. Table S2. qPCR primer sequences Table S3. qPCR on globin-depleted and non-globin-depleted mRNA show consistency of mean Cq Table S4. RNA sequencing of housekeeping and qPCR reference genes at HH29 shows consistent expression. Table S5. Comparison of targeted genes by RNA sequencing and qPCR at HH29 and HH35. Table S6. Molecular pathway analysis of significant differentially regulated genes in OFT-banded hearts. Table S7. Fragment per kilobase per million (FPKM) of genes of biological interest from RNA sequenced OFT-banded and sham control hearts Table S8. RNA sequencing gene regulation (OFT-banded vs sham) by DEseq