It is intractable to sample myliobatid embryos at specific developmental stages, as there is no published developmental timeline for any myliobatid taxon and no appropriately comparable model taxa have been characterized. Therefore, we conducted a pilot study to delineate the window of pectoral fin development in cownose rays from a wild population in Chesapeake Bay. Then, working with the bycatch of commercial fishermen, we sampled during July and August of 2016 and 2017. Embryos used for RNA-sequencing were preserved in RNALater, incubated at 4°C for 24–72 h, and then stored at −20°C until dissection and RNA extraction. Samples used for staging were fixed in 4% paraformaldehyde for 48 h before being transferred to methanol and stored at −80°C. Embryos were staged in an RNase-free environment, grouped by characters that we observed to be consistently associated with development of the pectoral fins (Table 1), and assigned a numerical index representing these characters (Table 2).

One goal of this study was to evaluate homology between the anterior pectoral fins of skates and the cephalic lobes of myliobatids using comparative transcriptomics. Because developmental-genetic programs can shift in space and time, it is important to examine similar developmental stages in each taxon. Therefore, morphological characters that were shared between developing cownose ray and winter skate (Leucoraja ocellata, following Maxwell et al., 2008) pectoral fins were identified. These were used to determine the most comparable stages of development in the cownose ray, Rhinoptera bonasus (R. bonasus), and little skate, Leucoraja erinacea (L. erinacea) (Table S1) for gene expression analysis.

Of all stages collected, the stage three R. bonasus embryos we sampled for RNA-sequencing represented the most biological replicates (six), and therefore gave us the greatest power to detect differentially expressed genes (DEGs). In addition, the pectoral fins undergo significant developmental processes during this stage, including anterior expansion, similar to stage 31 L. erinacea embryos sequenced by Nakamura et al. (2015). Due to the enhanced power to detect DEGs, the active developmental processes (including anterior expansion of the pectoral fins), and the availability of comparable data in L. erinacea, we focused our molecular analyses on stage three of R. bonasus development.

Nakamura et al. (2015) sequenced RNA from three equal-sized pectoral fin domains of developing embryos of L. erinacea, a representative batoid from the skate family (Rajidae), which diverged from all other batoids approximately 180 million years ago (Aschliman et al., 2012). We used these data for comparison to address the question of cephalic lobe/pectoral fin homology. To replicate the methods of Nakamura et al. (2015), we cut the left pectoral fin of stage three R. bonasus embryos (which includes the cephalic lobe) in thirds. We then added a layer of resolution to our dissections by bisecting each third of the pectoral fin again. The cephalic lobe comprised a distinct domain when the fin was split in this way (Figure 1). Further, when we included RNA-sequencing reads from adjacent dissections, our data represented three evenly-sized pectoral fin domains (following Nakamura et al., 2015; see Figure 1) with which to evaluate homology with the three domains previously examined in L. erinacea.

Following dissection, RNA was extracted from the six R. bonasus pectoral fin domains, as well as two additional fin domains, including pelvic fin. We extracted RNA using Trizol and followed with column-based purification using the PureLink RNA mini kit (Invitrogen) and a subsequent DNase I digestion. Purity was assessed using a NanoDrop with 260/280 ratios registering between 2.0 and 2.5. When there was residual guanidine salt, as indicated by a 260/230 ratio < 1.5, a SPRI bead purification was conducted and the eluate was reanalyzed to confirm that the final 260/230 ratio was greater than 1.5.

Sample aliquots were diluted and RNA integrity numbers (RIN) were evaluated using an RNA Pico chip on an Agilent BioAnalyzer. All samples exhibited RIN scores >8. Total RNA from each sample was quantified using a Qubit and 1 μg of RNA from each sample was used for Poly(A) mRNA enrichment using NEBNext oligo d(T) magnetic beads. After mRNA isolation and fragmentation, cDNA libraries were constructed using the NEBNext Ultra Directional Library Prep Kit for Illumina and each sample was given a unique barcode. Double-sided size selection was performed using SPRI beads targeting 400 bp fragments. Fragment size, concentration, and library purity were verified using a High Sensitivity DNA chip on a BioAnalyzer. In total, 48 libraries representing eight distinct fin domains with six biological replicates were normalized, pooled, and sequenced on an Illumina HiSeq 4000 using 100 bp paired end (PE) reads. Final library quality assessment and sequencing were conducted at the University of California, Berkeley; all other labwork occurred at San Francisco State University.

The cownose ray is a non-model organism, and the closest taxon for which a published genome exists is the elephant shark (Venkatesh et al., 2014), which diverged from the cownose ray and other elasmobranchs about 420 million years ago (Inoue et al., 2010). When sequence divergence is >15%, de novo transcriptome assemblies are recommended over reference-based mapping approaches (Vijay et al., 2013). To aid the de novo assembly process, we sequenced additional cephalic lobe, clasper, and pectoral fin domains on an Illumina MiSeq with longer reads (250 and 300 bp PE). The longer reads facilitated de novo transcriptome assembly by helping to bridge low complexity and redundant regions to which the shorter HiSeq reads could map. These samples were not included in our DGE analyses because they were prepared and sequenced using distinct kits and methodologies, which may introduce significant biases in expression levels (van Dijk et al., 2014).

The stages of cownose ray fin development documented here are roughly analogous to stages 28–31 in bamboo sharks and catsharks (Ballard et al., 1993; Onimaru et al., 2018), though we note that sharks lack the anterior expansion of the pectoral fins characteristic of the Batoidea. This anterior expansion appears to occur via similar mechanisms across batoid clades. In Raja eglanteria (Luer et al., 2007), Leucoraja ocellata (Maxwell et al., 2008), Urobatis halleri (Babel, 1966), and Rhinoptera bonasus (this study), pectoral fins begin developing as outgrowths from the body wall, as in all vertebrates. In skates and U. halleri, the anterior pectoral fins adopt a hook-like morphology and continue developing by growing away from the body (anterodistally) as they fuse to the pharyngeal basket in a rostral progression. The process terminates when the anterior pectoral fins fuse to the rostrum anterior to the mouth. Similarly, the anterior pectoral fins of the cownose ray begin developing posterior to the gill arches and assume a hook-like morphology. In the cownose ray, these hooks correspond to the developing cephalic lobes (Figure 2A). As development proceeds, the cephalic lobes unfurl but remain distinguished from the rest of the pectoral fins by a physical “notch,” where tissue outgrowth appears to be inhibited (Figure 2B). Concurrent with the unfurling of the cephalic lobes, the pectoral fins expand anterodistally and fuse to the body, beginning at the posterior pharyngeal arch and progressing rostrally (Figure 2C). Development to this point parallels the process observed in other batoids. However, once the region of the notch fuses at the shoulder near pharyngeal arch I, two processes commence in cownose ray pectoral fins that may be unique to myliobatids.

First, pockets of craniofacial tissue (CFT) begin developing lateral to the mouth (Figure 2D). These CFT pockets have ectodermal components that extend posteriorly to cover the gill arches, and distally to form a triangular sheet which eventually fuses to ectoderm of the pectoral fins. Second, the pattern of fusion near the notch begins to deviate ventrally to position the cephalic lobes within the CFT pockets on the ventral side of the body near the mouth. As development continues, the CFT pockets expand (Figure 2E), envelop (Figure 2F) and fuse to (Figure 2G) the cephalic lobes in a medial-distal progression that mirrors the caudal-rostral progression of pectoral fin fusion. By stage seven, the embryonic cownose ray has assumed the shape of an adult: fusion of pectoral fins, cephalic lobes, and CFT envelopment is complete, and only allometric growth remains. Surprisingly, the entire process of pectoral fin fusion spans just 2–3 weeks during the first month of the eleven-month cownose ray gestation period (J. Swenson, pers. obs.), indicating that the myliobatid body plan is established very early in development.

In cownose rays, cephalic lobes do not develop from independent fin buds emanating from the head; rather, they develop as modifications to the anterior pectoral fins, closely resembling the hooks at the anterior of developing skate and ray pectoral fins (Babel, 1966; Luer et al., 2007; Maxwell et al., 2008). However, unlike the anterior pectoral fins of other batoids, cephalic lobes are distinguished from the rest of the pectoral fin by a small region of reduced tissue outgrowth we call the “notch.” Although most of the transient morphologies and processes observed during cephalic lobe development are also found in developing pectoral fins of other batoids (e.g., the anterior hook, fusion to the pharyngeal arches and rostrum), the notch has only been observed in cownose rays in association with cephalic lobe development, suggesting that the notch is a defining feature of cephalic lobe specification.

Another distinction between developing pectoral fins in other batoids relative to cownose rays is manifest in the pattern of tissue fusion. Whereas the anterior pectoral fins in other batoids fuse to the rostrum in the same plane as the pectoral fins, in cownose rays, the pattern of fusion involves a distinct ventral deviation while the pectoral fins are fusing to the shoulder region corresponding to the notch. Ultimately, the unfused cephalic lobes are positioned on the ventral side of the body lateral to the mouth and enveloped within pockets of craniofacial tissue. This process may be unique to myliobatids with ventrally-positioned cephalic lobes.

After observing that myliobatid cephalic lobes develop as modified anterior pectoral fins, we used RNA-Sequencing (Wang et al., 2009) to compare gene expression profiles among pectoral fin domains of two batoids—a rajid skate (Leucoraja erinacea, little skate) and a myliobatid ray (Rhinoptera bonasus, cownose ray)—to determine whether homologous genetic underpinnings pattern the pectoral fins of these species and to reveal patterns that may be shared among divergent batoid lineages.

Using de novo transcriptome assembly and reciprocal-BLAST annotation with multiple vertebrates including Danio rerio, we confirmed expected gene expression patterns associated with the specific tissue domains dissected in R. bonasus (Figure S2). We then evaluated differential expression of genes in each R. bonasus pectoral fin domain relative to the others (Figure 3; Table S4). We found expression of Hoxd10-12 and Grem1 to be confined to the posterior pectoral fin, while Hoxd9 and Hoxa9 were most highly expressed in the mid-posterior pectoral fin, consistent with canonical patterning of paired appendages in vertebrates (Khokha et al., 2003; Zákány et al., 2004; Freitas et al., 2007; Ahn and Ho, 2008; Figure 3). We found genes associated with posterior patterning of pectoral fins in L. erinacea enriched in the posterior of R. bonasus pectoral fins as well (Figure 4, green dots), confirming that our samples were undergoing similar patterning processes and that the stages were comparable with respect to fin development (though note that the posteriorly expressed Ler-Grem1 transcript did not assemble well and was therefore not detected). Further, we confirmed that expression of Tbx4 and Tbx5, genes which define developing pelvic and pectoral fins in vertebrates, respectively (Gibson-Brown et al., 1998; Ahn et al., 2002), exhibited expected expression profiles in R. bonasus (for example, average counts per million in Tbx4 was >20x higher in the pelvic fin (280) when compared to pectoral fin domains (10); whereas the opposite was true with Tbx5). Taken together, these observations indicate that the cownose ray RNA-Seq data reliably show domain-specificity and indicate that similar processes drive development of the posterior pectoral fins of skates and myliobatids.

To evaluate homology between the anterior pectoral fin regions of myliobatid stingrays and rajid skates based on shared expression of anterior patterning genes, we identified genes known to be associated with patterning anterior regions of paired appendages and examined their expression in the anterior pectoral fin region of L. erinacea (anterior third, hereafter APF_Ler) and the comparable region in R. bonasus (hereafter APF_RboC; see Figure 1), which includes the cephalic lobes. By comparing ratios of gene expression in these domains relative to the posterior domains (within each species separately), we identified genes that are similarly enriched in the anterior and posterior pectoral fins of R. bonasus and L. erinacea (Figure 4). These include genes that are known to be associated with anterior patterning in fins and limbs such as Pax9, Alx4, Alx1, and Tbx2 (Onimaru et al., 2015) as well as Foxc2, Runx1, Meis2 and Dmrt3 (Figure 4, orange dots). We also found upregulation of Hoxa2 and Hoxd4 in the anterior pectoral fin regions of both species (Figure 4), which are likely associated with a batoid-specific AER in this region (Nakamura et al., 2015). Taken together, these data support homology of the anterior pectoral fin regions—which includes cephalic lobes in the Myliobatidae—in two disparate batoid species with distinct morphologies.

By sequencing six pectoral fin domains in the cownose ray instead of three, we added a layer of spatial resolution to our analyses which also revealed fine-scale expression patterns that were not apparent when considering broad domains for comparison with the little skate. This experimental design revealed high levels of expression of Hoxd8, Ntrk2, and Enpp2 at the anterior of developing cownose ray fins corresponding to the cephalic lobe (Figure 5; Table S5). In addition, the increased degree of spatial resolution allowed us to find enrichment of Wnt3 and Hoxa13 in the anterior region of R. bonasus pectoral fins (Figure 5), genes which are also enriched in the anterior of L. erinacea pectoral fins as indicated by in-situ hybridization (Nakamura et al., 2015; Barry and Crow, 2017). Wnt3 is associated with a novel apical ectodermal ridge (AER) driving outgrowth at the anterior of the L. erinacea pectoral fin (Nakamura et al., 2015). Hoxa13, a gene involved in patterning novel paired fin domains in batoids (O'Shaughnessy et al., 2015; Barry and Crow, 2017), has overlapping expression with Wnt3, suggesting a role for specification of the anterior AER in batoids (Barry and Crow, 2017). Enrichment of Hoxa13 and Wnt3 in the cephalic lobes of R. bonasus suggests that the anterior AER found in skates is correlated with development of the anterior pectoral fin region (including cephalic lobes) in derived batoids as well. We note that the high expression of these genes in the anterior pectoral fin region of R. bonasus was initially masked (see, for example, Hoxa13 in Figure 5), because they are upregulated in the cephalic lobe (red in Figure 5) but downregulated in the anterior (blue in Figure 5), resulting in signal cancellation when the expression scores for these domains were merged. Importantly, we found no evidence for independent posterior patterning by 5′ HoxD or Shh genes in cephalic lobes, supporting the notion that cephalic lobes do not represent a third set of independent paired appendages, but are instead modified anterior pectoral fins.

Unique morphological features arise in association with changes to developmental processes and gene expression patterns. Therefore, genes that exhibit distinct expression profiles in R. bonasus and L. erinacea during development of the pectoral fins may be associated with morphological differences. We have observed that cephalic lobe development is correlated with at least two surprisingly subtle modifications to the process observed in skate pectoral fins: establishment of a “notch” that defines the cephalic lobe-pectoral fin boundary, and a dorsoventral deviation during pectoral fin fusion. Because pectoral fin fusion is common to both skates and cownose rays, it is difficult to identify genes associated with fusion disparities using our comparative approach. However, we were able to identify candidate notch genes by searching our DGE results for genes that have different expression patterns in the APF_RboC and the APF_Ler relative to the posterior pectoral fin (i.e., up or downregulated in the anterior of one species but not the other, see pink genes in Figure 4). We then searched the literature for functional data related to these genes in other taxa and highlighted those known to be associated with processes that likely underlie development of the notch (e.g., cell proliferation, AER signaling) or the compagibus laminam (e.g., calcification, chondrogenic fusion) in Table 3.

Using these criteria, we identified six genes that may be associated with distinct morphologies at the anterior region of the pectoral fins of R. bonasus and L. erinacea (Table 3). Dkk1 is a known inhibitor of the Wnt signaling pathway (Glinka et al., 1998; Fedi et al., 1999; Niida et al., 2004), while Vsnl1 downregulates β-catenin—an important component of canonical Wnt signaling—in mouse kidneys (Ola et al., 2011). These genes are upregulated in R. bonasus but not L. erinacea and may have a role in notch development by inhibiting AER-signaling in this region. Alternatively, Lhx2 and Msx1 promote AER signaling and appendage outgrowth during development (Pizette et al., 2001; Lallemand, 2005; Tzchori et al., 2009). These genes are enriched at the anterior of L. erinacea and not R. bonasus; downregulation of these genes in R. bonasus may preclude AER-driven outgrowth at the notch. Omd is associated with biomineralization (Kawada et al., 2006; Tashima et al., 2015) and may contribute to development of the thicker fin rays in the compagibus laminam of myliobatids. RNA-Sequencing is a powerful tool for identifying genes that may be associated with novel or modified features during development; however, functional studies would be necessary to validate these associations.

While the CFT pockets and cephalic lobes are fusing at the anterior of the pectoral fins in cownose rays, specification of claspers begins at the posterior pelvic fins where they are distinguished by a notch that is reminiscent of the notch that occurs between the cephalic lobes and pectoral fins (see Figures 2B,D). As the claspers differentiate from the pelvic fins, they assume a paddle shaped morphology (Figures 2D,E) before curling into tubes (Figures 2F,G), transiently resembling the curled resting state of cephalic lobes in Mobula. As both claspers and cephalic lobes are modified fin domains supported by fin rays, the morphological similarities we observed during development led us to ask whether similar molecular processes could underlie development of these features.

In addition to the pectoral fin domains described above, we sequenced the transcriptome of pelvic fins of the six stage three cownose rays. At this stage, clasper morphogenesis has not yet begun, but we were able to diagnose sex based on expression of genes known to be exclusively expressed during clasper development in males (i.e. not expressed in pelvic fin development of females) such as Hoxd12, Hoxd13, Hand2 (O'Shaughnessy et al., 2015) and Hoxa13 (Barry and Crow, 2017). We confirmed statistically significant enrichment of Hoxd13, Hoxa13, and Hand2 in the pelvic fins of three of our samples relative to the others (Table 4), indicating that we had sampled three males and three females. We attributed other differentially expressed genes between male and female pelvic fin domains to clasper development, finding enrichment of Sall1, Cdh23, Sall3, Kit, Six2, and Cyp7a1 in male pelvic fins, adding to the repository of genes likely to be associated with clasper development (but note that Sall1 is expressed in posterior pectoral fins as well, see Figures 3, 5).

Interestingly, we found enrichment of five genes in both claspers and cephalic lobes, suggesting a redeployment of outgrowth and patterning pathways during development of these paired fin modifications. In addition to their enrichment in male pelvic fins relative to female, Hand2 and Sall1 are enriched in the cephalic lobes relative to the adjacent pectoral fin domains (Table S5). Ntrk2 is similarly upregulated in claspers (logFC = 1.54, p-value = 0.0002, adj. p-value = 0.13) and highly expressed in cephalic lobes (Table S4). Finally, Hoxa13 is expressed specifically in developing claspers of male pelvic fins, and is also expressed in anterior pectoral fins of the little skate (Barry and Crow, 2017) and cownose ray (Figure 5). Each of these patterns was obscured when examining broad domains in the cownose ray due to low expression in the domain adjacent to cephalic lobe. However, our experimental design allowed us to detect enrichment of these genes specifically in the anterior pectoral fin region corresponding to cephalic lobes (Figure 5). Although we did not detect enrichment of putative notch genes (those outlined in Table 3) in the pelvic fins of stage three males, this may simply reflect temporal differences in notch development, as the cephalic lobe notch is clearly defined at stage three, while the pelvic fin notch does not appear until stage four (Figure 2D). Intriguingly, androgen receptor (AR), another gene associated with clasper specification (O'Shaughnessy et al., 2015), was also upregulated in the cephalic lobe relative to every other pectoral fin domain (p < 0.01) (see Figure 5 and Table S5). AR was not represented in the L. erinacea transcriptome, so the possibility remains that AR expression at the anterior pectoral fin may be common among batoids. Future studies examining expression of genes associated with clasper development in the pectoral fins and cephalic lobes/APF of batoids may reveal that these paired fin modifications rely on similar developmental-genetic pathways for specification of novel domains.

Batoids are distinguished from other cartilaginous fishes by their disc-shaped pectoral fins, which are dorsoventrally flattened and extend anteriorly to fuse to the rostrum (Last et al., 2016). Species of the family Myliobatidae exhibit pectoral fin modifications during development that result in cephalic lobes. Until this study, the process by which cephalic lobes develop, including their physical attachment to developing pectoral fins, had not been documented. Interestingly, the presence of cephalic lobes is not correlated with an increase in the number of fin rays articulating to the propterygium in myliobatids (Hall et al., 2018), as would be expected if cephalic lobes were novel additions to ancestral anterior pectoral fins. Rather, the morphological observations and developmental gene expression patterns presented here suggest that cephalic lobes of myliobatids are homologous to the anterior pectoral fin of the little skate–a representative batoid lacking cephalic lobes. These data indicate that the mechanisms underlying anterior expansion of the pectoral fins were present in the common ancestor of skates and rays. The primary difference between developing cownose ray cephalic lobes and skate pectoral fins is the presence of a notch that defines cephalic lobes by separating the anterior pectoral fin into distinct domains. This subtle developmental modification is associated with specification of a distinct region - the cephalic lobe - that is optimized for feeding. The evolution of cephalic lobes is likely associated with modifications to the rest of the pectoral fins that permitted the evolution of oscillatory swimming (a.k.a. underwater flight) and adaptation to a pelagic or bentho-pelagic lifestyle in the Myliobatidae.

While cephalic lobes are homologous with anterior pectoral fins in skates, there are differences in developmental gene expression patterns associated with unique morphological modifications in myliobatid pectoral fins. Comparative transcriptomics using RNA-sequencing is limited in its capacity for identifying genes that are correlated with distinct phenotypes, as it does not include functional information. Nevertheless, by searching the data for genes that are up or down-regulated at the anterior of one species and not the other, we identified candidate genes that may be associated with myliobatid-specific features such as the compagibus laminam and development of the notch (Table 3; Table S6).

The expression pattern of osteomodulin (omd) is particularly intriguing. In R. bonasus, expression of omd is concentrated at the anterior of the pectoral fin, posterior to the cephalic lobe in the region of the compagibus laminam (Figures 3, 5). The pattern differs from L. erinacea, where omd is slightly upregulated at the posterior (Figure 4). Due to its known role in biomineralization and enriched expression in a pectoral fin domain of R. bonasus that exhibits increased calcification, we hypothesize that omd may be associated with development of the compagibus laminam, a plate of nearly fused and highly calcified fin rays at the anterior of myliobatid pectoral fins (Schaefer and Summers, 2005; Hall et al., 2018) which likely contributes to lift generation during oscillatory swimming.

The notch, a developmental feature that maintains a clear delineation between the cephalic lobe and the rest of the pectoral fin, likely develops as a region that is marked for reduced tissue outgrowth, or interruption of the AER, relative to the ancestral state. Once we identified genes with different expression patterns in the anterior pectoral fins of R. bonasus and L. erinacea, we queried this list for genes that have been shown to be associated with developmental pathways and processes that may be pertinent to notch formation (e.g., Wnt signaling, cell proliferation, apoptosis). Of the identified candidate genes, Dkk1, Msx2, and Lhx2 are compelling candidates for notch induction. Dkk1 is a known Wnt inhibitor. Wnt signaling is an important component of AER induction and maintenance (ten Berge et al., 2008). Upregulation of Wnt inhibiting genes at the notch could provide a mechanism by which AER-facilitated outgrowth could be stunted in the notch while continuing in the adjacent regions. Alternatively, downregulation of genes such as Msx1 and Lhx2 that promote AER maintenance could have a similar effect, resulting in reduced cell proliferation (Table 3). Though we did not identify any candidate notch genes with roles in apoptosis, apoptosis could also be a mechanism by which a notch could form in developing fins. If apoptotic pathways were activated in the notch, AER signaling could occur as a continuous strip throughout the anterior pectoral fin (including cephalic lobe), while an increase in cell death at the notch could prevent tissue outgrowth from keeping pace with adjacent areas. Further research into notch development is needed to confirm the underlying mechanisms, perhaps via functional studies that attempt to induce a notch in a more tractable batoid species such as the little skate (Leucoraja erinacea) via ectopic expression or deletion of the genes highlighted in Table 3.

Development of the anterior region of myliobatid pectoral fins appears to involve a unique combination of Hox gene expression and redeployment of posterior modules. Though most HoxD genes exhibit canonical expression patterns in the pectoral fin of R. bonasus (Figure 5), expression of Hoxd8 is conspicuously enriched in the APF/cephalic lobe of R. bonasus but not L. erinacea (Figures 4, 5), suggesting that there may be a myliobatid-specific Hox code at the anterior of the pectoral fin. In addition, Enpp2 (aka autotaxin), is a known target of Hoxa13 (Williams et al., 2005) and may contribute to AER induction in the anterior pectoral fins of batoids (this gene was not represented in the L. erinacea transcriptome), or specification of cephalic lobes in myliobatids. Mirroring the pattern of Enpp2, Ntrk2 expression is enriched in both the cephalic lobes and the posterior pectoral fins of R. bonasus, as are Wnt3 and Hand2 (Figure 5). By revealing genes with expression patterns that are mirrored at the anterior and posterior of R. bonasus pectoral fins, our domain-specific data suggest that these regions share a handful of common patterning mechanisms, which may contribute to development of boundaries in paired fins.

Modifications to outgrowth and fusion processes may contribute to the morphological variation observed at the anterior pectoral fins among myliobatid taxa. For example, the single fused cephalic lobe of mature bat rays is found in the same plane as the pectoral fin, indicating that the observed deviation during pectoral fin fusion in the region of the notch in cownose rays is lacking in bat rays (Figures 6A,B). This appears to be due to a dorsoventral shift in the pattern of fusion at the notch region in cownose rays and other myliobatid species. While species in the myliobatid genera Rhinoptera, Mobula, and Aetobatus lack skeletal elements in the shoulder (corresponding to the notch), some bat rays of the genus Myliobatis display a radiation of stunted fin rays in this region between the cephalic lobe and pectoral fin (Figures 6C,D). We hypothesize that, relative to other myliobatids, these bat rays exhibit weak or temporally constricted expression of the genes responsible for suppressing tissue outgrowth in the notch. Combined, these observations suggest that the molecular processes driving notch formation and the pattern of fusion may vary among myliobatid taxa, with representation of intermediate phenotypes.

Aberrant phenotypes reflecting incomplete pectoral fin fusion have been observed in various batoid clades, including the Myliobatidae (Ramírez-Amaro et al., 2013), Rajidae (Prado et al., 2008), Torpedinidae (Mnasri et al., 2010), Dasyatidae (Prado et al., 2008; Ramírez-Hernandez et al., 2011), Urotrygonidae (Mejía-Falla et al., 2011), Potamotrygonidae (Oldfield, 2005a,b), Rhinobatidae (Bornatowski and Abilhoa, 2009), and Gymnuridae (Narváez and Osaer, 2016). This prevalent deformity was first described over a century ago (Gill, 1896) and has since been termed the “batman” morphology (Oldfield 2005a; Oldfield 2005b). After observing that a two-step process consisting of distinct phases of outgrowth and fusion underlies development of the anterior pectoral fins in cownose rays, we suggest that the batman morphology may occur when the molecular mechanisms driving outgrowth proceed normally, while the mechanisms underlying fusion are interrupted. It has been suggested that the batman morphology may be associated with exposure to pollutants (Heupel et al., 1999), which could disrupt developmental-genetic processes (Tomanek, 2011; Varotto et al., 2013).

Claspers are male reproductive organs in cartilaginous fishes, which are modifications to the posterior pelvic fins in extant taxa. Morphologically, clasper formation is preceded by development of a notch between the budding claspers and pelvic fins, reminiscent of the notch that develops between the cephalic lobes and the rest of the pectoral fins. This led us to wonder if specification of cephalic lobes in anterior pectoral fins represents a mirror image of clasper specification in posterior pelvic fins, and if we could find evidence for co-option of this ancient developmental pathway in other modified fin domains.

Transcriptomic data revealed that genes known to be involved in clasper specification are indeed enriched in developing cephalic lobes. Barry and Crow (2017) found that Hoxa13 is expressed exclusively in the claspers of the little skate, and not expressed in female pelvic fins. Variation in the pattern of Hoxa13 expression is often associated with the evolution of novel features (Crow et al., 2009; Archambeault et al., 2014; Barry and Crow, 2017), including the autopod (Sordino et al., 1995; Yokouchi et al., 1995; Fromental-Ramain et al., 1996; Metscher et al., 2005), making it a likely candidate for specification of derived batoid features as well. During clasper development, Hand2 is regulated by hormonal input via upregulation of AR (O'Shaughnessy et al., 2015); it was suggested that the evolution of claspers was facilitated by hormonal input to fin developmental networks following the evolution of androgen response elements at the Hand2 locus (O'Shaughnessy et al., 2015). Interestingly, both Hand2 and AR are similarly enriched in developing cephalic lobes. In conjunction with high levels of Ntrk2 and Sall1 expression in both claspers and cephalic lobes (relative to the adjacent fin domain; see Figure 5, Table S5), our data support an intriguing hypothesis: namely, that cephalic lobes, which are anterior pectoral fin modifications, may have evolved by co-opting developmental pathways associated with the evolution of claspers - posterior pelvic fin modifications in male cartilaginous fishes.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.