Wild type (CD1) mice were time-mated to obtain tissue for microdissection. Animal experiments followed local legislations regarding housing, husbandry, and welfare (ASPA 1986; United Kingdom). Embryos were collected at E10.5 and E11.5 and accurately staged by counting somites. BA1, BA2 and PBA/OFT tissues were dissected and snap frozen on dry ice and stored at −80°C until RNA extraction. The OFT was harvested with the PBA as it acts as a landmark feature during dissection and maintains physical integrity of the PBA. BAs from three closely staged embryos were pooled for each library, with BA1, BA2 and PBA/OFT from the same embryos being used for each timepoint replicate. Total RNA was extracted using the miRNeasy micro kit (Qiagen, #217084) following the manufacturer’s instructions, eluted into RNase-free water and stored at −20°C until library preparation. Small RNA libraries were generated using the NEBNext Small RNA Library Prep Set (New England BioLabs, #E7330S) following manufacturer’s instructions. For size selection we used gel separation and extracted amplified microRNA cDNA bands corresponding to 140bp. NEBNext Index primer sequences and respective libraries are listed in Supplementary Table S1. Libraries were quality checked using the Agilent 2200 BioAnalyzer TapeStation and sequenced on the Illumina HiSeq 4000 at The University of Manchester Genomics Technologies Core Facility. Small RNA-seq libraries have been deposited under project accession PRJEB64007 available from the European Nucleotide Archive.

NEBNext Small RNA adapter sequences were removed using cutadapt v1.8 (Martin, 2011). Adapter-trimmed reads were filtered to keep those 18–25 nt in length, and mapped against mouse tRNAs [mm 10, GtRNAdb v18.1 (Chan and Lowe, 2016)] and rRNAs [Mus musculus, Silva SSU/LSU r138.1 (Ludwig et al., 2004; Quast et al., 2013)], using Bowtie v1.1.0 (Langmead et al., 2009). Reads that mapped to tRNAs/rRNAs were discarded from further analysis. Remaining reads were mapped to the mm10 primary fasta file (GRCm38, release M23) using bowtie v1.1.0 (Langmead et al., 2009) with the following settings: bowtie -v1 -a -m5 –best–strata. To predict novel microRNAs, we used miRDeep2 v0.1.3 (Friedländer et al., 2011) with combined filtered reads from all our small RNA-seq libraries. Reference microRNAs included stem-loop mouse, mature mouse and mature rat microRNAs, all downloaded from miRBase v22 (Kozomara et al., 2019). Novel pre-microRNAs were filtered using the following criteria: ≥30 0-mismatch reads for the mature arm and ≥10 0-mismatch reads for the star arm, no internal sub-hairpins, ≥50% 5′ arm homogeneity, 0-4 nt overhang at the 3′ arm, hairpin free energy ≤ −0.2 kcal/mol/nt. Custom Python scripts used for filtering are available at github.com/SianGol. Pre-microRNA sequences that met all the above criteria were input into Rfam v14.6 sequence-search (Kalvari et al., 2020) to remove any that overlapped with previously annotated ncRNAs. Remaining novel microRNAs were added to the GTF file of known microRNAs, downloaded from miRBase v22 (Kozomara et al., 2019). Genome-mapped reads were assigned to mature microRNAs and quantified using featureCounts v1.6.0 (Liao et al., 2014) with the following settings: featureCounts -M -g gene_id -s 0.

BA1, BA2 and PBA/OFT E10.5 and E11.5 (Losa et al., 2017) RNA-seq libraries were adapter-trimmed and quality-filtered using Trimmomatic v0.36 (Bolger et al., 2014). We used STAR v2.5.3a (Dobin et al., 2012) to generate the mm 10 genome index (GRCm38, release M23) and map RNA-seq reads to the genome using the mm 10 genome primary fasta file and corresponding GTF file, both downloaded from GENCODE (Frankish et al., 2021). Mapped reads were then assigned to annotated genomic features and quantified using featureCounts v1.6.0 (Liao et al., 2014), with options set as following: featureCounts -g gene_id -s 2.

To remove lowly-expressed genes, a threshold was set for our RNA-seq libraries: 4.2CPM (corresponding to approximately 100 reads) for small RNA-seq, and 0.38 CPM (corresponding to approximately 10 reads) for RNA-seq. Genes with CPM above these thresholds in two or more libraries were kept for further analysis. Differential analysis was performed using DESeq2 (Love et al., 2014) and results were obtained for the following pairwise comparisons: E10.5 BA1 vs. BA2, E10.5 BA1 vs. PBA/OFT, E10.5 BA2 vs. PBA/OFT, E11.5 BA1 vs. BA2, E11.5 BA1 vs. PBA/OFT, E11.5 BA2 vs. PBA/OFT.

3′UTR coordinates of BA transcripts were identified using the mm10 GTF and extract_transcript_regions.py (Floor, 2018). The longest 3′UTR for each gene was used, alongside BA expressed microRNAs, as inputs for in silico microRNA target prediction by seedVicious v1.3 (Marco, 2018). Results were filtered to remove predicted interactions with 6mers, off-6mers, sites with a hybridisation energy of >−7 kcal/mol, and 3′UTRs with <2 predicted microRNAs binding sites.

E10.5 PBA/OFT Gene Ontology (GO) terms were identified using PANTHER (Mi et al., 2018). The background genes used were those differentially expressed (±1.5-fold, adj-p ≤ 0.05) in any BA pairwise comparison, and the input genes were those >1.5-fold (adj-p ≤ 0.05) expressed in the PBA/OFT compared to BA1 and BA2. We defined our PBA/OFT gene set as those annotated under the most significant GO terms returned (≥2-fold, FDR ≤ 0.05) by PANTHER. To perform a hypergeometric test for enrichment, we used the phyper function in R, with PBA/OFT-gene-set-microRNA interactions as the ‘sample’, and all BA-microRNA interactions as the “population”.

NKX2-5^eGFP/w hESCs (Elliott et al., 2011) were seeded in hESC medium (DMEM/F-12 (Gibco, #31765027), 1X None-Essential Amino Acids (Gibco, #11140050), 1X GlutaMAX (Thermo Scientific, #35050038), 0.1 mM 2-ME (Gibco, #21985023), 0.5% penicillin-streptomycin (Sigma-Aldrich, #P0781), 20% KnockOut Serum Replacement (KSR) (Gibco, 10828028), and 10 ng/mL bFGF (Miltenyi, #130-104-924)) at a density of 1.8 × 10⁵ cells/mL on growth factor-reduced Matrigel coated 6-well plates. Twenty-four hours later (day 0) differentiation was induced as described previously (Giacomelli et al., 2020). hESC medium was replaced with BPEL (Ng et al., 2008) supplemented with BMP4 (Bio-techne/R&D, #314-BP-050), 20 ng/mL ACTIVIN A (Miltenyi Biotec, #130-115-009) and 1.75 μM CHIR99021 (Selleckchem, #S1263). On day 3, media was refreshed with BPEL containing 1 μM XAV939 (VWR, #CAYM13596-1). BPEL was refreshed every 3 days thereafter.

Candidate target 3′UTRs were amplified using the following reaction; 1X Q5 Reaction Buffer (New England BioLabs, #B9027S), 20 units/mL Q5 High-fidelity DNA Polymerase (New England BioLabs, #M0491S), 0.6 μM 3′UTR F/R primer (Supplementary Table S1), 0.2 mM dNTP mix (Promega, #U1511), 2 μL DNA, and RNase free water, using the cycling conditions: 98°C 1 min, 35X [98°C 10 s, 58°C 30 s, 72°C 1 min], 72°C 2 min. Amplified 3′UTRs were size selected and purified using the QIAquick Gel Extraction kit (Qiagen, #28704) following manufacturer’s instructions. Purified 3′UTRs were ligated with the pmirGLO dual-luciferase vector (Promega, #E1330) via SacI and SalI restriction sites. Resultant plasmids were sequenced using 0.4 μM pmirGLO F/R custom sequencing primers (Supplementary Table S1).

miR-92b-3p binding sites in each 3′UTR were mutated using the QuikChange II XL site-directed mutagenesis kit (Agilent, #200517) following the manufacturer’s instructions, using custom mutagenesis primers (Supplementary Table S1). Mutated plasmids were sequenced as described above.

NIH/3T3 cells were reverse co-transfected with 100 ng pmirGLO dual luciferase vector (Promega, #E1330) containing appropriate 3′UTRs and 30 nM microRNA mimic in 96-well plates (Invitrogen, miR-92b-3p #4464066, miRNA negative control 1 #4464058) using lipofectamine 2000 (Invitrogen, #11668019) diluted in opti-MEM (Gibco, #31985062). Luciferase activity was measured 48 h later using the Dual-Glo Luciferase Assay System (Promega, #E2920) following the manufacturer’s instructions. Firefly and Renilla luciferase signal were measured using the Promega GloMax-Multi + Detection System. Five technical replicates were performed for each biological replicate. Firefly luciferase values were first normalised to Renilla luciferase, and fold change was calculated relative to a control sample transfected with the dual-luciferase plasmid and no microRNA mimic.

HEK293 and NKX2-5^eGFP/w hESCs were reverse transfected with 30 nM microRNA mimic (Invitrogen, miR-92b-3p #4464066, miRNA negative control 1 #4464058) using lipofectamine 2000 (Invitrogen, #11668019) diluted in opti-MEM (Gibco, #31985062). Samples were incubated for 24 h at 37°C in 5% CO2.

Where hESC cardiomyocyte differentiation timeline samples were used for microRNA and mRNA RT-qPCR, total RNA was extracted using the mirVana miRNA Isolation kit (Invitrogen, #AM1560) following the manufacturer’s instructions. Alternatively, where only mRNA expression was measured, RNA was isolated using TRIzol Reagent (Invitrogen, #155996026) following a standard protocol. For miR-92b-3p quantification, 10 ng total RNA was used as input with the TaqMan Advanced MiRNA cDNA Synthesis Kit (Applied Biosystems, #A28007) following the manufacturer’s instructions. For our internal control, U6, 10 ng total RNA was used as input for the TaqMan MicroRNA Reverse Transcription Kit (Applied Biosystems, #4366596) following manufacturer’s instructions. To measure miR-92b-3p and U6 expression we used the TaqMan Fast Advanced Master Mix Kit (Applied Biosystems, #4444556) according to manufacturer’s instructions, with 1X TaqMan Advanced miRNA Assay (Applied Biosystems, #A25576, assay ID: 477823_mir) or 1X TaqMan small RNA Assay (Applied Biosystems, #4427975, assay ID: 001973) respectively. For mRNA expression we used the QuantiTect SYBR Green RT-PCR kit (Qiagen, #204243) in the following reaction: 1X QuantiTect SYBR Green RT-PCR master mix, 0.3  μM F/R primers (Supplementary Table S1), 1X RT mix, 40 ng RNA, RNase-free water. Cycling conditions: 50°C 30 min, 95°C 15 min, 40X [95°C 20 s, 57°C 30 s, 72°C 30 s], 68°C 7 min, 4°C hold. Fold change was calculated using 2^−ΔΔCt relative to the given control.

HEK293 cells were lysed (20 mM Tris-HCl, 120 mM NaCl, 0.5 mM EDTA, 0.5% NP40 (Thermo Scientific, #J60766-AP), 10% glycerol (Thermo Scientific, #17904), 1X protease cocktail inhibitor (Roche, #4693116001)) and the supernatant was recovered. Proteins were denatured in 1X laemmli buffer. Western blot membranes were incubated with 1:1500 anti-GATA-6 rabbit-mAb DEIE4 (Cellsignal, #5851) or 1:50,000 anti-β-Actin-peroxidase mouse-mAb (Sigma Aldrich, #A3854) in 1% milk. For GATA6 1:10,000 Goat anti-Rabbit IgG HRP (abcam, #ab6271) was used as a secondary antibody in 1% milk.

miR-92 sequences for human, mouse, zebrafish and Drosophila were obtained from miRBase v22 (Kozomara et al., 2019). Gata6 and Tbx20 3′UTR sequences were obtained from the UCSC genome browser versions hg38 and mm 10. To align sequences, we used Clustal Omega (Sievers and Higgins, 2014).

To identify microRNAs that function during development and morphogenesis of mammalian BAs, we generated small-RNA-seq libraries for BA1, BA2, and PBA/OFT tissue from embryonic (E) 10.5 and E11.5 mouse embryos. Both timepoints coincide with mid-gestation and allow us to capture the three separate BA domains prior to their fusion at E11.5-E12.5 (Frazer, 1926). Furthermore, these datasets complement RNA-seq libraries we previously generated from equivalent tissues and timepoints (Losa et al., 2017), and so provide a valuable opportunity to consider expression of both microRNAs and predicted target mRNAs.

BA small-RNA-seq libraries were enriched for 18–25 nt reads (Supplementary Figure S1A), and 73%–86% of the 18–25 nt reads mapped to known microRNAs (Table 1). As shown by principal component analysis (PCA) (Figure 1A), replicate samples clustered closely to one another indicating reproducibility. The sample clustering across the first two principal components coincided with anterior-posterior location and developmental stage respectively. Using Spearman’s rank correlation to determine similarity between datasets, BA1 E11.5 and BA2 E11.5 showed the greatest correlation, followed by PBA/OFT E10.5 and PBA/OFT E11.5 (Supplementary Figure S1B). The PBA/OFT were most distinct in microRNA expression compared to the two anterior BA domains.

Using miRDeep, and a set of post hoc filters for high confidence microRNA annotations, we identified 28 putative novel pre-microRNAs (Supplementary Table S2). The discovery of novel microRNAs is perhaps unexpected in such a well-studied model organism. However, it is well-known that many microRNAs show very specific special and temporal expression patterns (Xu et al., 2024), and therefore under-studied tissues and developmental time points continue to reveal novel microRNA loci. All these predicted loci pass accepted strict criteria, including for submission to the miRBase database (see methods), and as with all microRNAs, further datasets and studies will help to clarify these annotations. Two of the novel microRNAs share seed sequences with mmu-miR-702 and mmu-miR-1839.

After applying an expression cut-off (see methods), we found a total of 550 microRNAs were expressed across the BAs. We explored differential expression in both a spatial and temporal manner. To identify domain-specific microRNAs, we performed BA pairwise comparisons across E10.5 and E11.5 samples. We defined differentially expressed microRNAs as those with ≥1.5-fold change and adjusted p-value≤0.05. In at least one pairwise BA comparison for a given timepoint, 94 microRNAs were differentially expressed at E10.5, and 103 microRNAs were differentially expressed at E11.5. There was considerable overlap between the sets of differentially expressed microRNAs at E10.5 and at E11.5 (Figure 1B), showing that regional microRNA expression is largely maintained across these developmental timepoints. Additionally, at both E10.5 and E11.5, we saw increased microRNA expression progressively across the anterior-posterior axis, whereby most differentially expressed microRNAs demonstrated ≥1.5-fold expression in the PBA/OFT compared to BA1 and BA2 (Supplementary Figure S1C, D).

By considering the expression of microRNAs and their validated targets we can infer their regulatory outcome. For example, 11 mature microRNAs were significantly more highly expressed in BA2 compared to both BA1 and PBA/OFT at E10.5 (Figure 1C). These include miR-743b-3p, miR-741-3p, miR-878-5p, miR-881-3p, miR-871-3p, miR-470-5p, miR-465a/b/c-3p, and miR-465b/c-5p (Figures 1D, E). These microRNAs are transcribed from a large microRNA cluster, Fx-mir, spanning ∼62 kb on Chr X (Figure 1F). Members of this cluster have been described to regulate the neighbouring gene Fmr1 (Ramaiah et al., 2019; Wang et al., 2020), however, each respective study draws contradictory conclusions with regards to whether these microRNAs repress or promote expression of Fmr1. Expression profiles from our BA RNA-seq show that Fmr1 was most lowly expressed in BA2 (Figure 1G), suggesting Fx-mir microRNAs may have a repressive role on Fmr1 in BA2.

The PBA/OFT demonstrated the most distinct microRNA expression across the BAs, at both E10.5 and E11.5 (Supplementary Figure S1C, D). Some of the greatest differentially expressed microRNAs (Figure 1E; Supplementary Figure S1E) are located within the Hox clusters; miR-10b overlaps with both Hoxd3 and Hoxd4, miR-10a is located upstream of Hoxb4, and miR-615 is found within Hoxc5 (Figure 1H). These Hox genes are also more highly expressed in the PBA/OFT E10.5 (Figure 1G), coinciding with spatial collinear Hox gene expression in distinct streams of cranial NC cells that populate the BA’s (Frisdal and Trainor, 2014). Therefore, Hox cluster microRNAs in the PBA/OFT mirror expression of their overlapping or nearby protein-coding genes.

In summary, we have identified differentially expressed microRNAs across the BAs, and BA-specific expression patterns are largely maintained between E10.5 and E11.5. As E10.5 samples demonstrated more defined PCA clustering (Figure 1A) and had a greater number of differentially expressed microRNAs (Supplementary Figure S1C, D), we focused on this timepoint for the remainder of our study.

The microRNA biogenesis pathway has previously been implicated in PBA artery remodelling (Nie et al., 2011) and OFT morphogenesis (Saxena and Tabin, 2010). We combined expression data from both microRNA-seq and RNA-seq datasets to identify individual microRNA mediated regulation in this domain. The approach is shown in Figure 2A. To begin, we generated a comprehensive list of predicted interactions between BA microRNAs and BA protein-coding genes (see methods). In parallel, we curated a subset of protein-coding genes specifically implicated in PBA/OFT development. We identified this subset by extracting genes that were more highly expressed in the PBA/OFT, relative to BA1 and BA2, performing GO analysis and selecting genes annotated under the top GO terms. Many of these terms were related to muscle and cardiac development, consistent with the contribution of this area to the anterior pole of the heart (Figure 2B). We obtained a list of 315 genes associated with developmental processes specific to the PBA/OFT. To identify microRNAs most significantly enriched for predicted interactions with these 315 genes, we performed a hypergeometric p-value test. For each microRNA, this test considers its predicted interactions with all BA genes and its predicted interactions with the 315 genes specific the PBA/OFT. In return, we determine the microRNAs that interact with this set of genes, compared to all BA genes as background, more than would be expected by chance.

The top ranked microRNA by hypergeometric p-value was miR-92b-3p (p-value = 0.0029). This microRNA was 1.7-fold enriched for predicted interactions with the PBA/OFT list of genes compared to interactions with all BA expressed mRNAs. On average, miR-92b-3p was most highly expressed in the PBA/OFT (Figure 2C) and expressed above the 75th percentile of all PBA/OFT microRNAs (Figure 2D). Looking at its predicted interactions, miR-92b-3p was predicted to bind to 30 of the 315 PBA/OFT genes (Supplementary Table S3).

MicroRNAs are understood to stabilise GRNs through broad regulation of multiple targets (Liufu et al., 2017). Therefore, we were interested in identifying interactions between miR-92b-3p predicted targets. To do so we used the STRING database, which returned evidence for multiple functional interaction networks (Figure 2E). One of these protein-protein interaction clusters contained four developmental cardiac TFs, GATA6, TBX20, HAND1 and HAND2, which have previously been shown to have correlated expression (Sharma et al., 2020). In human induced pluripotent stem cells (hiPSCs) undergoing cardiomyocyte differentiation, HAND1, HAND2 and TBX20 are downregulated following GATA6 loss of function (Sharma et al., 2020). This may be a result of direct regulation, as GATA6 binds in the vicinity (100 kb upstream or downstream) of Hand1, Hand2, and Tbx20 in the PBA/OFT (Losa et al., 2017) (Figures 2F–H).

Taken together, we generated in silico microRNA-target predictions using time-matched BA RNA-seq expression datasets. We identified miR-92b-3p as enriched for predicted binding sites for a network of cardiac transcription factors. As one of these transcription factors, Hand2, has previously been validated as a target of miR-92b-3p (Yu et al., 2019), we explored the possibility that miR-92-3p may regulate the other predicted targets Gata6, Tbx20, and Hand1 to broadly modulate a cardiac GRN.

To assess whether miR-92b-3p can interact with its predicted target sites in Gata6, Tbx20 and Hand1 3′UTRs, we used dual-luciferase reporter assays. Whilst we note these assays may not reflect physiological binding, our aim for this experiment was determine if binding can occur at the target sites and if this brings about reduced reporter expression. Portions of each 3′UTR containing the predicted binding sites were cloned into the Promega pmirGLO dual-luciferase plasmid (Figure 3A). Following co-transfection of miR-92b-3p mimic and Gata6 and Tbx20 3′UTR reporters, we saw a significant reduction in luciferase signal, approximately 0.5-fold, compared to the control mimic (Figures 3B, C). In contrast, when miR-92b-3p was co-transfected with the Hand1 3′UTR reporter, a smaller reduction in reporter signal was seen compared to the control mimic (Figure 3D). These results are consistent with the predicted hybridisation energies between miR-92b-3p and the three putative binding sites (Figures 3E–G), whereby Hand1 has the weakest interaction predicted at −7.2 kcal/mol, just meeting our −7 kcal/mol threshold cut-off.

Next, we determined whether luciferase reporter knockdown was specifically caused by the physical interaction of miR-92b-3p with the putative binding sites in Gata6 and Tbx20 3′UTRs. Initial binding occurs between positions 2–4 of a microRNA seed site and the target RNA (Chandradoss et al., 2015), therefore, we mutated nucleotides complementary to the seed region of miR-92b-3p (Figures 3E, F). Following co-transfection of mutated 3′UTRs with miR-92b-3p, we did not observe reduced luciferase reporter expression, confirming luciferase reporter knockdown occurs due to interaction of miR-92b-3p with the identified binding sites located in Gata6 and Tbx20 3′UTRs (Figures 3H, I).

In summary, miR-92b-3p interacts with the predicted target sites in Gata6 and Tbx20, leading to significantly reduced reporter expression. Furthermore, the extent of repression was reflective of the predicted hybridisation energy.