Prenatal Developmental Trajectories of Fluctuating Asymmetry in Bat Humeri

Fluctuating asymmetry (random fluctuations between the left and right sides of the body) has been interpreted as an index to quantify both the developmental instabilities and homeostatic capabilities of organisms, linking the phenotypic and genotypic aspects of morphogenesis. However, studying the ontogenesis of fluctuating asymmetry has been limited to mostly model organisms in postnatal stages, missing prenatal trajectories of asymmetry that could better elucidate decoupled developmental pathways controlling symmetric bone elongation and thickening. In this study, we quantified the presence and magnitude of asymmetry during the prenatal development of bats, focusing on the humerus, a highly specialized bone adapted in bats to perform under multiple functional demands. We deconstructed levels of asymmetry by measuring the longitudinal and cross-sectional asymmetry of the humerus using a combination of linear measurements and geometric morphometrics. We tested the presence of different types of asymmetry and calculated the magnitude of size-controlled fluctuating asymmetry to assess developmental instability. Statistical support for the presence of fluctuating asymmetry was found for both longitudinal and cross-sectional asymmetry, explaining on average 16% of asymmetric variation. Significant directional asymmetry accounted for less than 6.6% of asymmetric variation. Both measures of fluctuating asymmetry remained relatively stable throughout ontogeny, but cross-sectional asymmetry was significantly different across developmental stages. Finally, we did not find a correspondence between developmental patterns of longitudinal and cross-sectional asymmetry, indicating that processes promoting symmetrical bone elongation and thickening work independently. We suggest various functional pressures linked to newborn bats’ ecology associated with longitudinal (altricial flight capabilities) and cross-sectional (precocial clinging ability) developmental asymmetry differentially. We hypothesize that stable magnitudes of fluctuating asymmetry across development could indicate the presence of developmental mechanisms buffering developmental instability.

Fluctuating asymmetry (random fluctuations between the left and right sides of the body) has been interpreted as an index to quantify both the developmental instabilities and homeostatic capabilities of organisms, linking the phenotypic and genotypic aspects of morphogenesis. However, studying the ontogenesis of fluctuating asymmetry has been limited to mostly model organisms in postnatal stages, missing prenatal trajectories of asymmetry that could better elucidate decoupled developmental pathways controlling symmetric bone elongation and thickening. In this study, we quantified the presence and magnitude of asymmetry during the prenatal development of bats, focusing on the humerus, a highly specialized bone adapted in bats to perform under multiple functional demands. We deconstructed levels of asymmetry by measuring the longitudinal and cross-sectional asymmetry of the humerus using a combination of linear measurements and geometric morphometrics. We tested the presence of different types of asymmetry and calculated the magnitude of size-controlled fluctuating asymmetry to assess developmental instability. Statistical support for the presence of fluctuating asymmetry was found for both longitudinal and cross-sectional asymmetry, explaining on average 16% of asymmetric variation. Significant directional asymmetry accounted for less than 6.6% of asymmetric variation. Both measures of fluctuating asymmetry remained relatively stable throughout ontogeny, but cross-sectional asymmetry was significantly different across developmental stages. Finally, we did not find a correspondence between developmental patterns of longitudinal and cross-sectional asymmetry, indicating that processes promoting symmetrical bone elongation and thickening work independently. We suggest various functional pressures linked to newborn bats' ecology associated with longitudinal (altricial flight capabilities) and cross-sectional (precocial clinging ability) developmental asymmetry differentially. We hypothesize that stable magnitudes of fluctuating asymmetry across development could indicate the presence of developmental mechanisms buffering developmental instability.

INTRODUCTION
Bats (Chiroptera) have a highly specialized postcranial morphology adapted to perform self-powered flight, a trait unique among mammals. Compared to terrestrial mammals, the forelimb skeleton of bats is characterized by aerial-based functional adaptations that provide structural support to the wing (bone elongation of distal bones) (Panyutina et al., 2015), withstand biomechanical loading (increased bone density) (Swartz, 1997;Dumont, 2010), enhance flight maneuverability (proximo-distal demineralization) (Papadimitriou et al., 1996), and facilitate the generation of flight power (increase of muscle mass in the shoulder girdle) (Panyutina et al., 2015). Moreover, morphological adaptations related to different functional pressures can occur within a single bone (e.g., increased muscle attachment area and biomechanical load resistance in the humerus; Panyutina et al., 2015). Such a variety of long bone adaptations arise ontogenically, mainly through two bone remodeling processes: bone elongation and bone thickening (Enlow, 1963). The former is linked to axial bone deposition via endochondral ossification at growth plates (Webber et al., 2015), whereas the latter results from endosteal and periosteal bone deposition via intramembranous ossification in the periosteum/endosteum (Kronenberg, 2003). These two processes derive from independent morphogenetic pathways that respond to different factors, indicating decoupled developmental trajectories controlling the optimal development of individual forelimb bones (Atchley and Hall, 1991;Montoya-Sanhueza et al., 2019). However, our understanding of the differences in developmental trajectories of bone thickening and elongation in the bat wing remains largely limited.
Uncovering the historical trajectories that led to the morphological diversification and specialization of bats has been greatly limited by a markedly incomplete fossil record (Brown et al., 2019). Evolutionary developmental biology has emerged as a promising approach to study the evolution of bats while circumventing limitations in the fossil record (Adams, 2008;Cooper et al., 2012;Camacho et al., 2019;López-Aguirre et al., 2019,a,b). Studies have provided evidence for the ontogenetic mechanisms behind forelimb specialization in bats and the evolution of vertebrate flight (Sears et al., 2006;Adams, 2008;Cretekos et al., 2008;Farnum et al., 2008;Cooper et al., 2012;Adams and Shaw, 2013), ecology-driven deviations in chiropteran development from general mammalian patterns, and the phylogenetic signal in postcranial development (Adams, 1992;Koyabu and Son, 2014;López-Aguirre et al., 2019,a,b).
For body plans that are naturally symmetrical, deviations from an "ideal" state have been interpreted as a signal of reduced fitness (Dongen, 2006). Accordingly, quantifying the deviation from perfect symmetry can indicate the amount of stress an organism undergoes and its homeostatic capacity (i.e., buffering of instabilities to maintain fitness) (Gummer and Brigham, 1995;Aparicio and Bonal, 2002). The main regulatory mechanisms that influence phenotypic symmetry occur ontogenically, when the genotypic and phenotypic mechanisms involved in morphogenesis can be destabilized by genetic or environmental stressors (i.e., developmental noise) (Hallgrímsson, 1998(Hallgrímsson, , 1999Hallgrímsson et al., 2003;Kellner and Alford, 2003). Evolutionary studies have provided evidence for the heritability of an organism's capacity to buffer developmental noise (developmental stability/instability, DI), suggesting that natural selection can act as a regulator of phenotypic asymmetry (Tocts et al., 2016).
Phenotypic asymmetry in animals with bilateral asymmetry can be quantified by computing the morphometric differences between the right and the left side of the body (rightleft, R-L) (Palmer, 1994). Advances in theoretical framework to quantify body asymmetry has enabled the recognition of three types of asymmetry, each with a different biological interpretation: fluctuating asymmetry (FA) is characterized as random deviations from "ideal" perfectly symmetrical phenotypes, directional asymmetry (DA) is described as a natural tendency to have consistently asymmetrical phenotypes (i.e., one side always larger than the other), and antisymmetry (AS) represents a pattern where symmetrical phenotypes are least favored and asymmetry is equally distributed across both sides (Klingenberg, 2015). All three types of asymmetry have also been described in mathematical terms: FA is characterized by a normal distribution of asymmetry values (R-L) along a value mean of zero, DA is described by a normal distribution of asymmetry values along a mean different to zero, and AS is identified where asymmetry values have a bimodal distribution and most values are different from zero (Palmer, 1994;Klingenberg, 2015). Combining the biological and mathematical interpretations of phenotypic asymmetry, FA has been regularly used as a possible indicator of DI, although some studies argue that DA and AS can also be indicators of DI (Palmer, 1994;Leamy and Klingenberg, 2005). Despite the utility of studying FA and DI, the efficacy of the theoretical framework traditionally applied to detect a real FA-DI link has been debated (Palmer, 1994;Dongen, 2006;Klingenberg, 2015).
The FA-DI link can vary in response to ecological, genetic, environmental, and developmental factors, stressing the need to study it at multiple scales to test a variety of hypotheses (Kellner and Alford, 2003). Swaddle and Witter (1997) summarized and Kellner and Alford (2003) postulated and tested predictions for a list of hypotheses on the ontogeny of FA describing developmental patterns of asymmetry and possible evolutionary mechanisms shaping them. These include the following: small fluctuations during early growth are magnified during later morphogenesis (i.e., magnification of asymmetry hypothesis), side-biased environmental influences can skew growth toward asymmetrical phenotypes (i.e., directional external cues hypothesis), accumulative growth of independent subunits will homogenize morphogenesis, reducing asymmetry throughout development (i.e., coin toss hypothesis), and developmental feedback mechanisms will stabilize asymmetric growth between structures by either promoting or constraining growth (i.e., compensatory growth hypothesis). Studies on the developmental basis of FA have also been restricted in scope (mostly focused on postnatal development) and study groups (invertebrates and captive populations) (Hallgrímsson, 1999;Hallgrímsson et al., 2003;Leamy and Klingenberg, 2005;Blackburn, 2011;Perchalski et al., 2018), limiting our understanding of variation in wild non-model species and potential insights into the mechanisms controlling FA in early development. Based on the limited studies available, decreasing magnitudes of FA across prenatal development have been reported as a result of the interaction between variations in timing of growth and growth rates (Hallgrímsson, 1998;Hallgrímsson et al., 2003;Kellner and Alford, 2003). In contrast, magnitudes of FA have been reported to increase during the development of the mammalian skeleton due to cumulative variability in growth regulation and/or bone remodeling (i.e., morphogenetic drift model; Hallgrímsson, 1998).
Asymmetry in bats has been studied at the cranial and postcranial level (Juste et al., 2001a,b;Lüpold et al., 2004;Voigt et al., 2005;López-Aguirre and Pérez-Torres, 2015), all studies being based on the analysis of adult specimens only. Forelimb FA has been associated with differential reproductive success (Lüpold et al., 2004;Voigt et al., 2005), suggesting that sexual selection favors symmetric individuals in Saccopteryx bilineata (Voigt et al., 2005) and a significant link between asymmetry and reproductive potential in Carollia perspicillata (Monteiro et al., 2019). A correlation between wing FA and resistance to environmental stress and resilience to anthropogenic habitat change has also been assessed, suggesting high resilience in Neotropical bat species (de Figueiredo et al., 2015;Castillo-Figueroa, 2018). Compensatory growth has been reported in the wing of the vampire bat Desmodus rotundus as a way to maintain wingspan symmetry (Ueti et al., 2015), while sex-based differences in the magnitudes of FA have been reported in the wing of D. rotundus and the cranium of Artibeus lituratus (López-Aguirre and Pérez-Torres, 2015; Ueti et al., 2015). Research also indicates that variation in the levels of FA across morphological (both cranial and postcranial) traits could depend on functional importance, with FA decreasing in traits under higher functional demands for feeding and locomotion (Gummer and Brigham, 1995;López-Aguirre and Pérez-Torres, 2015;Robaina et al., 2017). Despite the repeated study of bat forelimb FA, its ontogenetic basis and how it varies across different bone growth dimensions (i.e., bone elongation and thickening) remain unknown.
The objective of this study was to assess whether the ontogenetic trajectories of phenotypic asymmetry during bone elongation and thickening in bats are decoupled. We analyzed the presence and magnitude of FA in the prenatal morphogenesis of the humerus in bats, representing the first developmental study of FA in Chiroptera. We focused on the humerus because it represents a clear example of multiple functional demands acting on a single bone (i.e., withstanding torsional and bending stress, increasing muscle insertion area, and controlling the maneuverability of the wing) (Swartz et al., 1992;Panyutina et al., 2015). Humeral cross-sectional shape has been found to reflect foraging differences across bat taxa (López-Aguirre et al., 2019). Furthermore, prenatal limb FA has been described as an accurate indicator of DI in human fetuses (Klingenberg and Nijhout, 1999;Broek et al., 2017). We quantify asymmetry based on bone elongation and cross-sectional cortical bone deposition as a way to exemplify the multipatterned process of bone growth (Klingenberg and McIntyre, 1998). Humeral length asymmetry is commonly used in bat studies (Gummer and Brigham, 1995;Voigt et al., 2005;de Figueiredo et al., 2015;Ueti et al., 2015;Robaina et al., 2017;Castillo-Figueroa, 2018), whereas crosssectional asymmetry is commonly measured in other mammals (Macintosh et al., 2013;Wilson and Humphrey, 2015;Perchalski et al., 2018). Based on previous studies inferring compensatory growth in the wing of bats and previous evidence of decreasing FA across development in terrestrial mammals (Hallgrímsson et al., 2003;Ueti et al., 2015), we hypothesize that asymmetry of the humerus will decrease throughout ontogeny. Because bone thickening and elongation in mammals represent independent morphogenetic trajectories, we predict that cross-sectional and longitudinal asymmetry will not be correlated.

Sampling
A total of 66 prenatal specimens from 11 bat species (Table 1, see Figures 1A-C) were collected through taxonomic fieldwork in Vietnam by VTT and DK under collection permit no. 972/UBND-TH issued by Tuyen Quang Provincial People's Committee and research and ethics permit no. 322/STTNSV of the Institute of Ecology and Biological Resources, Vietnam Academy of Sciences. All specimens were fixed in Serra's fixative (ethanol, formalin, and glacial acetic acid mixed 6:3:1 by volume) for 48 h, then transferred, and preserved in 70% ethanol. 3D scanning of the embryos and fetuses was performed using a microfocal X−ray computed tomography system at the University Museum, University of Tokyo (TXS225−ACTIS; TESCO; Tokyo, Japan), with 70-kV source voltage and 114 µA source currents at a resolution of 36 µm. All osseous skeletal elements were segmented by the first author using the thresholding tool and the predetermined bone setting in MIMICS v. 20 software (Materialise NV, Leuven, Belgium). To standardize finer manual segmentation in early embryos, thresholding of Hounsfield unit values of osseous tissue was performed using the half-maximum  height method (i.e., gradual change in computerized tomography values at the boundary of a structure) (Spoor et al., 1993). This method homogenously retrieved diaphyseal osseous tissue only, removing segmenting inconsistencies related to the imaging of diaphyseal cartilaginous tissue.
Smoothing techniques were not used to prevent artificially changing the dimensions of the periosteal surface of the models. After selecting specimens with at least partially ossified and unbroken humeri, three individuals were excluded from the final sample. The left and right humeri were segmented and exported as STL files for further processing and analysis.

Data Collection
Humeral development was described based on bone elongation (i.e., humeral length, HL) and cortical bone deposition (i.e., periosteal diameter as a proxy for cross-sectional cortical thickness, CT) ( Figure 1D). To estimate cross-sectional and HL measurements, 3D humeri models were imported into Rhinoceros 5.0 (Robert McNeel & Associates, Seattle, WA).
To remove the non-shape effects of translation, rotation, and scale, all humeri models were aligned to a standard position in 3D space, following a protocol for long-bone cross-sectional geometry (see Wilson and Humphrey, 2015). HL was automatically measured as the length of the long axis of a rectilinear box (i.e., bounding box) enclosing the model created using the BoundingBox command in Rhinoceros 5.0, preventing human measurement error (ME). Crosssections of the left and right humeri were extracted at the midshaft (i.e., 50% of HL) for a total of 126 crosssections [63 × 2 (left and right sides)]. The midshaft was extracted by placing a perpendicular axis intersecting the 3D model at the midpoint of the length of the rectilinear bounding box, dividing HL in half. Small deviations in the placing of the cross-sections could occur due to asymmetrical development of cartilaginous epiphyseal tissue. Humeral CT was quantified using a geometric morphometrics-based approach.
Following the method described in Wilson and Humphrey (2015), we used a set of 16 equiangular landmarks semiautomatically placed along the periosteal surface of each cross-section ( Figure 1D). The cross-sections were aligned with a set of 16 equiangular radii along the centroid of each cross-section, aligning the radii with the anatomical axes of the bone (e.g., radii 1 and 9 represent the anteroposterior axis). The landmarks were automatically placed on the intersection of individual radii and the periosteal surface of each cross-section. CT was quantified as the average of interlandmark distances between pairs of landmarks that formed linear axes (e.g., landmarks 1 and 9; see Figure 1D) using the interlmkdist function in R package Geomorph 3.2 . Our automated cross-sectional geometric morphometrics-based approach enables the quantification of CT while circumventing the lack of identifiable homologous landmarks early in prenatal development. To control for the effect of matching bilateral symmetry (i.e., the left and right sides of the body are mirror images), the landmark coordinates of the cross-sections of the right side were reflected along the antero-posterior axis by multiplying the coordinates of that axis by −1. Given the lack of Carnegie staging systems for many non-model taxa (nine of the 11 species in our sample), staging of developmental series was based on crownto-rump length and bone ossification sequence as described in López-Aguirre et al. (2019a) and following the general patterns described in bat development (Cretekos et al., 2005). All developmental stages were represented by at least three individuals (Supplementary Table 1).

Estimation of Asymmetry
Asymmetry in bilateral organisms can be described as the difference between both sides of the body (e.g., right-left) (Palmer, 1994), which we estimated based on HL (longitudinal asymmetry) and CT (cross-sectional asymmetry). Individual longitudinal and cross-sectional asymmetries were quantified as the signed difference between right and left HL, negative and positive values indicating the directionality of asymmetry (Palmer, 1994). Individual longitudinal FA was described as the size-corrected signed difference between sides, and per stage FA was calculated as the variance of individual FA measurements across a population (i.e., FA6 index) (Palmer, 1994). FA6 is a single-trait index that expresses FA as the variance in asymmetry proportional to trait size (i.e., HL) in an individual so as to truly represent DI and not developmental bone growth, ].

Data Analysis
The presence of FA and DA in longitudinal and cross-sectional asymmetry was statistically tested with full-factorial ANOVAs using side (left or right), individual, and duplicate as factors (FA ∼ side + individual + side/individual; see Table 2; Monteiro et al., 2019;Rivera and Neely, 2020). The side factor provides a statistical test for DA, whereas the side-individual interaction provides statistical tests for FA. Measurement error was not computed because the automated protocol implemented to obtain HL and CT ensures that no human error could affect the measuring process. Developmental trajectories of individual longitudinal and cross-sectional FA were explored using box plots, while per stage FA6 was explored using bar plots. Statistical differences across developmental stages were tested using two ANOVAs, using the stage of each individual as a factor and unsigned longitudinal and cross-sectional FA (longitudinal FA∼stage and cross-sectional FA∼stage).
We tested the association of longitudinal and cross-sectional FA across development using a linear regression model based on ordinary least squares. We used linear regression models to test the effect of the number of specimens and species in values of average longitudinal and cross-sectional FA per developmental stage to assess whether uneven sample composition affected our results. Using a subsample of the best sampled genus, we explored the association between peaks in humeral growth and peaks in

Individual Longitudinal and Cross-Sectional Asymmetry
Distribution of individuals' values of longitudinal and crosssectional FA demonstrate a normal distribution with a mean near zero (Figure 2), supporting the presence of FA in longitudinal and cross-sectional humeral asymmetry throughout development. Moreover, 42.86% of individuals had negative values of longitudinal FA (Figure 2A), indicating that a narrow majority of individuals had larger right humeri (57.14%). Furthermore, 39.68% of individuals showed negative values of cross-sectional asymmetry ( Figure 2B). Across datasets (crosssectional and HL), only three individuals were found to have perfectly symmetrical humeri (i.e., R-L = 0) for HL in developmental stages 5, 7, and 8. Box plots of individuals' longitudinal and cross-sectional FA values illustrated independent trajectories of cross-sectional and longitudinal humeral FA across development (Figure 3). Unsigned longitudinal FA indicated two peaks of high values of longitudinal FA in individuals early in prenatal development (stages 2, 3, and 6), separated by a sharp decrease in longitudinal FA values between stages 4 and 5. Longitudinal FA values steadily decreased from stage 7 onwards, with a noticeable increase in dispersion in stages 2 and 6 ( Figure 3A). A similar pattern was found in our subsample of the best sampled species, with relatively stable longitudinal FA after stage 1 and an increase in variability in stages 8 and 9 (Supplementary Figure 1). Crosssectional FA showed low magnitudes of FA in individuals in the early stages of prenatal development (stages 1 and 2), followed by an increase between stages 3 and 5 ( Figure 3B). Individuals in stage 5 showed the highest dispersion of FA values. Crosssectional FA consistently decreased from stage 6 onwards. The cross-sectional FA in our subsample remained relatively stable after stage 3, followed by a decrease from stage 7, similar to our general results (Supplementary Figure 1).

Presence and Magnitude of Fluctuating Asymmetry
The ANOVA of asymmetry on HL supported the presence of both DA and FA in longitudinal asymmetry ( Table 2). Nevertheless, FA accounted for only 0.007% of variation, DA for less than 0.0001% of the variation, and individual variation for 99.25% of variability in raw HL. The markedly low values of longitudinal FA suggest that the extraction of the cross-sections did not vary due to asymmetrical growth rates between epiphyseal plates. The ANOVA of cross-sectional asymmetry found statistical support for the presence of FA and DA ( Table 2). DA accounted for only 13.20% of cross-sectional variation and FA for 31.69% of variation, respectively. Individual variation explained the highest  proportion of variation, accounting for 55.11% of cross-sectional asymmetrical variation.

Developmental Trajectories of Fluctuating Asymmetry
The per-stage longitudinal and cross-sectional FA (FA6) indicated disparate trajectories across development (Figure 4). The longitudinal FA values showed two peaks ( Figure 4A): one early in development (stage 2) and another in intermediate prenatal development (stage 6). The lowest longitudinal FA values were found in stages 4 and 8-10. Cross-sectional FA showed a single peak during mid-prenatal development (stage 5), followed by a significant decrease from stages 6 to 10 ( Figure 4B). Despite the differences in the trajectories of longitudinal and cross-sectional FA throughout development, our results indicate a trend of stable FA throughout development. The ANOVAs for differences in longitudinal and cross-sectional FA across  developmental stages found statistically significant differences in cross-sectional FA but not in longitudinal FA (Table 3). Our results do not reveal clear similarities between peaks of FA and peaks of humeral growth, indicating that asymmetry does not increase when growth rates are higher (Supplementary Figure 2). Nevertheless, future studies should further explore the possible link between growth rates and magnitudes of asymmetry during development. The scatterplots of longitudinal and cross-sectional asymmetry did not indicate a clear association between cross-sectional and longitudinal FA throughout development, with a slightly negative tendency (Figure 5). A linear regression model of cross-sectional and longitudinal FA confirmed a non-significant negative correlation between both dimensions of humeral FA (R 2 = 0.016, P = 0.979). The cross-sectional and longitudinal FA also did not correlate with the number of specimens and species per stage (cross-sectional FA: specimens R 2 = 0.136, P = 0.159 and species R 2 = 0.183, P = 0.121; longitudinal FA: specimens R 2 = 0.096, P = 0.66 and species R 2 = 0.072, P = 0.546), indicating that sampling heterogeneity did not influence our results.

DISCUSSION
This study is the first to demonstrate the presence of FA and DA during bat prenatal development. Our results showed stable magnitudes of longitudinal and cross-sectional FA across prenatal development, rejecting our hypothesis of decreasing asymmetry throughout ontogeny. However, our results also show that both measures of FA (longitudinal and cross-sectional) did not correlate across development, showing decoupled ontogenetic trajectories, similar to our prediction. Decoupled longitudinal and cross-sectional asymmetry magnitudes also indicate that the developmental pathways regulating symmetrical bone elongation and thickening likely act independently within a single bone. Signaling pathways regulating symmetric growth have been detected in the apical ectodermal ridge and the zone of polarizing activity of the limb bud (Allard and Tabin, 2009;Wolpert, 2010). Sonic hedgehog and fibroblast growth signaling factors play a major role during limb growth, controlling limb polarity (i.e., anterior-posterior limb differentiation), cell proliferation, and symmetric growth (Allard and Tabin, 2009;Wolpert, 2010). Previous studies have found decreasing trajectories of asymmetry across prenatal development in the postcranium of humans and mice, suggesting an interplay between the timing of morphogenesis and growth rates over the course of ontogeny (Hallgrímsson, 1998;Hallgrímsson et al., 2003). The compensatory growth between left and right sides of the body in response to increased asymmetry has been discussed as a mechanism to reduce DI (Aparicio, 1998;Ueti et al., 2015).
FA in bats has been an area of increasing interest in recent decades (Gummer and Brigham, 1995;Juste et al., 2001a,b;Lüpold et al., 2004;Voigt et al., 2005;de Figueiredo et al., 2015;López-Aguirre and Pérez-Torres, 2015;Ueti et al., 2015;Robaina et al., 2017;Castillo-Figueroa, 2018;Monteiro et al., 2019). A set of two studies explored the patterns of FA in insular populations of fruit bat species Eidolon helvum and Rousettus egyptiacus of the Gulf of Guinea (Juste et al., 2001a,b). Juste et al. (2001a) found similar patterns of FA across populations of both species, discussing the interpretation of population asymmetry parameters (i.e., consistent patterns of FA for a set of characters across populations of the same species), and suggested its scalability above the species level. Juste et al. (2001b) found consistent magnitudes of multivariate FA and a significant integration of asymmetry across traits and species, hypothesizing FIGURE 5 | Scatterplot of association between longitudinal (unsigned fluctuating asymmetry, FA) and cross-sectional (unsigned FA) fluctuating asymmetry across prenatal development. Colored dots represent developmental stages, with stages 1-10 representing early to late prenatal development.
high canalization in the developmental pathways controlling phenotypic asymmetry that are shared between the two species.
The presence and magnitude of FA has also been associated with reproductive success and sexual selection in bats (Voigt et al., 2005;Monteiro et al., 2019). Analyzing forearm length asymmetry, it has been suggested that sexual selection favors more symmetrical males in the polygynous greater sac-winged bat Saccopteryx bilineata, canalizing FA (Voigt et al., 2005). The number of offspring produced by males of S. bilineata was found to significantly decrease with increasing forearm asymmetry (Voigt et al., 2005). Increases in forearm asymmetry in the Neotropical frugivorous bat Carollia perspicillata have also been linked with a significant decrease in survival probability and the probability of more than one pregnancy per reproductive season (Monteiro et al., 2019). However, Lüpold et al. (2004) did not find a significant association between FA and other measures of individual fitness and allometry in the bat genitalia. All previous studies analyzing humeral asymmetry in bats were based on adult specimens with fully ossified epiphyses, whereas our study could not analyze cartilaginous osseous tissue.
Studies have also linked the presence and magnitude of FA to anthropogenic perturbations and habitat degradation, analyzing four Neotropical phyllostomid species (i.e., Artibeus lituratus, Artibeus planirostris, C. perspicillata, and Sturnira lilium) (de Figueiredo et al., 2015;Castillo-Figueroa, 2018). Neither study explicitly tested the presence of the three kinds of asymmetry, casting doubt on their interpretation that FA is an accurate index of resilience and adaptation of species to perturbations. Robaina et al. (2017) may have also insufficiently assessed the presence and magnitude of all types of asymmetry before drawing conclusions on the validity of FA to reflect functional importance. Our results suggest that special attention to the statistical framework used to describe the biological and theoretical interpretation of asymmetry in bats is warranted.
Asymmetry across multiple traits has shown decoupled patterns, reflecting functional differences and indicating independent developmental mechanisms controlling phenotypic symmetry across different structures of the body of bats (Gummer and Brigham, 1995;Robaina et al., 2017), turtles (Rivera and Neely, 2020), and birds (Aparicio and Bonal, 2002). Our results of decoupled patterns of longitudinal and crosssectional humeral FA point toward independent trajectories of bone elongation (endochondral ossification) and thickening (intramembranous ossification) while also suggesting that it could be applicable within single structures. Longitudinal and cross-sectional growth are hypothesized to be decoupled ontogenetic processes (Enlow, 1963;Montoya-Sanhueza et al., 2019) and to respond to different selective pressures, with cross-sectional bone deposition potentially associated with biomechanical resistance against torsional and bending stresses (Blackburn, 2011;Perchalski et al., 2018) and bone elongation correlating with maintenance of body proportions within a functional unit (Ueti et al., 2015).
Peaks of longitudinal (stages 2 and 6) and cross-sectional FA (stage 5) at different developmental stages could indicate variations in the timing of ossification onset and/or growth rates across species, following the morphogenetic drift model (Hallgrímsson, 1998). Despite bats showing a general mammalian developmental pattern (López-Aguirre et al., 2019a), heterochronies and diverging allometric trajectories for the ossification of the humeri have been found across bat species (López-Aguirre et al., 2019b). The reported slower ossification of the humerus in yinpterochiropteran bats could indicate that different growth rates at a given developmental stage could result in increased magnitudes of FA (López-Aguirre et al., 2019b). However, our study does not allow for the accurate comparison of growth rates across species, reflecting the challenge of amassing embryonic material for non-model species, as not all species have complete developmental series. Future studies should focus on further exploring the relationship between interspecific variation in growth rates and magnitudes of FA. Despite both longitudinal and cross-sectional FA showing homogenous magnitudes across development, only cross-sectional FA showed statistically significant differences across developmental stages. Increasing postcranial morphological disparity and integration across prenatal development in bats has been reported in ossification sequences and metric growth (López-Aguirre et al., 2019a).
We hypothesize that our results (significant differences in cross-sectional FA across stages) indicate a greater selective pressure to canalize cross-sectional asymmetry as a response to functional demands associated with a newborn bat's ecology. In particular, there may be a greater requirement for symmetrical cross-sectional growth to facilitate early roosting behaviors. In many bat species, newborns attach to the mothers using their feet and thumb (Koyabu and Son, 2014) rather than to immediately fly, with bone elongation continuing during this period. Multiple studies have documented postnatal development of flight in bats, describing the altricial development of the forelimb in bats followed by accelerated bone elongation (Hughes and Rayner, 1993;Kunz and Robson, 1995;Kunz et al., 2009;Lin et al., 2010;Eghbali et al., 2017;Eghbali and Sharifi, 2018). First flights in most species that have been studied occur consistently in synchrony with weaning, usually a couple of weeks after birth once adult body dimensions are reached (Hughes and Rayner, 1993;Kunz and Robson, 1995;Kunz et al., 2009;Lin et al., 2010;Eghbali et al., 2017;Eghbali and Sharifi, 2018). We hypothesize that because self-powered flight is not achieved immediately after birth, bone elongation asymmetry and compensatory growth to optimize wing proportions would be less constrained prenatally (Ueti et al., 2015). Further studies quantifying the developmental trajectories of asymmetry should focus on describing the complementary developmental process (pre-and postnatal). Additionally, testing this hypothesis in mammal species with newborns that also cling to the mothers (e.g., macaques and colugos) would further elucidate whether the patterns found in this study are common across Mammalia.
We found significant support for the presence of FA and DA during the prenatal development of the humerus in bats. We also show that magnitudes of FA remain stable across prenatal development, and we hypothesize this to be evidence of developmental control of FA. Moreover, we find evidence for decoupled patterns of longitudinal and crosssectional asymmetry throughout prenatal humeral development. We suggest that decoupled morphogenetic processes (i.e., bone thickening via intramembranous ossification and bone elongation via endochondral ossification) and functional differences between bone elongation and cross-sectional bone deposition may be associated with the newborn's ecology (i.e., pup roosting behavior and the later acquisition of flight). To our knowledge, this study is the first to analyze asymmetry patterns in the development of bats, providing new information about phenotypic asymmetry and DI in non-model taxa. We highlight the importance of assessing the correlation between FA and DI beyond patterns of total asymmetry FA and DI.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The animal study was reviewed and approved by the Institute of Ecology and Biological Resources, Vietnam Academy of Sciences ethics permit No. 322/STTNSV.

AUTHOR CONTRIBUTIONS
CL-A, SH, and LW designed the study. CL-A gathered and analyzed the data. CL-A, SH, LW, and DK wrote the manuscript. DK and VTT sampled the specimens. DK prepared and scanned the specimens. CL-A processed the scans. All authors read and approved the final manuscript.