Introduction

Evolutionary changes in size are common over the course of animal diversification (Clauset and Erwin, 2008; Hanken and Wake, 1993), strongly impacting many aspects of morphology, physiology, and behaviour. Investigating how changes in body size affect other traits (i.e., scaling) is thus essential to understand the causes of traits variation in animals. Traditionally rooted in the field of morphology and physiology (Gould, 1966; Hill, 1950; Pélabon et al., 2014), allometric studies have broadened to explore how more complex traits scale with animal size. These include locomotor kinematics (e.g. flight, stride or swimming kinematics, Bejan and Marden, 2006; Riskin et al., 2010), and behaviour (e.g. mating or escape tactics, Dial et al., 2008). Although behavioural or biomechanical traits are often more noisy and challenging to measure, their inclusion in allometric studies is crucial to better understand trait evolution (Cloyed et al., 2021). Indeed, most phenotypes involve both morphological and behavioural traits, each responding in potentially different ways to size variation (Clemente and Dick, 2023). Variation in scaling between body size and other traits is generally most apparent between species, and can be partly explained by neutral divergence during phylogenetic history (Blomberg et al., 2003). Understanding the cause of trait variation among species thus requires controlling for the influence of phylogeny on allometric relationships.

In this study, we address the consequences of size variation between species (i.e. evolutionary allometry), and specifically question the relative importance of wing morphology versus wingbeat kinematic changes in mitigating the effect of size on flight ability.

Flight typically combines a suite of morphological and biomechanical traits, providing a relevant context to jointly examine how morphology and locomotor kinematics change with size. The physical laws governing flight allow drawing predictions on the scaling of morphology and kinematics in flying animals.

These physical laws set both upper and lower bounds on the size of flying animals. For large flying animals, the aerodynamic lift for weight support scales linearly with the surface area of the wing and quadratically with the forward speed of the flying animal(Norberg, 2007). Physical scaling laws dictate that, for isometric scaling, wing surface area scales with (body mass2/3) (Schmidt-Nielsen, 1975). Thus, due to constraints for maintaining weight support, larger flying animals exhibit disproportionately larger wings and tend to fly faster (Norberg, 2012). This positive allometry in wing surface area and flight speed with body mass is needed to generate enough lift to stay in the air.

Smaller flying animals such as insects generally fly slower, and thus air movements over the wing are governed by wing flapping, especially in hovering flight. The aerodynamic lift forces produced by flapping wings therefore scale linearly with the second-moment-of-area (S2) of the wing, and quadratically with the angular velocity of the beating wing (C. P. Ellington, 1984). For isometry, the second moment of wing area scales with (body mass)4/3. As this growth factor is larger than one, wing area decreases at higher rate compared to body weight, and therefore the aerodynamic force reduces more rapidly than body weight. Disproportionately larger wings and/or a higher wingbeat frequency is thus expected to evolve in smaller flying insects as a compensatory mechanism. The negative allometric scaling in wingbeat frequency with body mass is supported by the consistent negative correlation between wingbeat frequency and size observed in most flying animals (Norberg, 2007; Rayner, 1988; Riskin et al., 2010; Tercel et al., 2018). A general trend toward disproportionately larger wings in smaller species is however less striking among species (but see Danforth, 1989; C. Ellington, 1984a; Outomuro et al., 2013)), suggesting that morphological changes in response to size variation are not always straightforward.

Although empirical data often align with expectations drawn from physical laws, the scaling of morphology and kinematic with size can sometimes be challenging to predict. Selective pressures associated with ecological specialisation may influence the evolution of morphology and kinematic, and conflict with the constraints imposed by physical scaling laws. Each component may then adjust for the constraints imposed by weight support in varying degree. For example, changes in morphology have been shown to predominantly mitigate the consequences of body size: among hummingbird species of various sizes, changes in wing area – but less in wingbeat kinematics – compensate for weight support (Skandalis et al., 2017). Similarly, wing area but not wingbeat frequency was found to covary with body size among species of bees (Duell et al., 2022; Grula et al., 2021). Such discrepancies with the general scaling expectation may point at a particular selective regime maintaining specialised flight kinematics and ability.

To what extent the evolution of wingbeat kinematics and wing morphology is influenced by physical scaling demands or by ecological specialisation is poorly understood. Here we focused on an ecologically specialised insect group, the hoverflies, to test if the consequences of size variation are mostly mitigated by changes in wingbeat kinematics or in wing morphology.

Hoverflies (Diptera: Syrphidae) are among the most species-rich groups of nectar-foraging insects and exhibit striking variation in size, ranging from 0.5 mg to over 200 mg (Katzourakis et al., 2001). Despite these large size differences, all adult hoverflies show exceptional hovering abilities. This highly specialised flight behaviour may have been promoted because it enhances foraging efficiency, as they hover in front of flowers for nectar feeding. It could also have been favoured because hovering plays an important role in hoverfly mating behaviour: male hoverflies aggregate at specific locations such as under a tree and hover steadily, waiting for passing females to pursue (Gilbert, 1984; Heinrich and Pantle, 1975). Hovering abilities may thus be correlated with mating success (Downes, 1969; Gilbert, 1984). The strong selection acting on the flight of hoverflies may thus interfere with general scaling expectations between body mass, wing morphology and wingbeat kinematics, especially during evolutionary miniaturisation. Tiny hoverflies require relatively large wings or high wingbeat frequencies to maintain weight support during hovering.

Here, we address the outcome of this evolutionary process by examining the scaling relationships between wing morphology and body mass in eight hoverfly species ranging from 5 to 100 mg, and conjointly assessing how wing morphology and wingbeat kinematics scale with body mass. We then combined computational fluid dynamic simulations with quasi-steady aerodynamic modelling to estimate aerodynamic forces produced by the different species in hovering flight. Based on this, we determined how the relative changes in wingbeat kinematics and wing morphology among the eight species contribute to producing weight support during hovering flight. Because flight behaviour is more flexible than morphology (Blomberg et al., 2003), we predict that changes in wingbeat kinematics will prevail over changes in wing morphology. Alternatively, if selection for specialised hovering flight is limiting kinematic flexibility, we expect weight support to be achieved through changes in wing morphology rather than in wingbeat kinematics.

Material and methods

Collecting and identifying hoverflies

We collected hoverflies from the wild, in September 2022 on the campus of Wageningen University & Research in the Netherlands (51°59’01.0”N 5°39’32.7”E; ca. 10 m). A total of eight hoverfly species were captured, including 44 individuals (n=5.5±2.9 individuals per species, mean±standard deviation) (Figure 1). Flies were captured with hand-held nets and brought in the lab to record their flight within an hour after capture.

The studied hoverfly species shown with their phylogenetic relationships. Phylogenetic positions of the studied species in a larger hoverfly phylogeny (91 species) are shown in the bottom left corner. Sample size as number of individuals for each species is in parentheses. Mean values (± standard deviation) for body mass are presented, along with a photography of a representative individual and a wing from each species. Phylogeny was obtained from (Wong et al., 2023). Photographs were taken by the second author.

Species were first distinguished visually, but their identity was subsequently ascertained using DNA barcoding. This was done by sending a piece of leg of each individual to the Canadian Centre for DNA Barcoding (CCDB), where barcoding was performed following standard automated protocols of the BOLD Identification System (http://www.boldsystems.org/)(Ratnasingham and Hebert, 2007). Sequenced DNA was then compared with a reference library to establish a species match. Sexes were determined visually using criteria from (Gilbert, 2015).

Flight experiments

The flight experiments were carried out in a custom-built octagonal flight arena made of transparent Plexiglas (50×50×48 cm, height×width×length; Figure 2A)(Cribellier et al., 2022). Three synchronised high-speed cameras (Photron FASTCAM SA-X2, set at 5000 frames s−1, spatial resolution of 1084×1084 pixels, shutter speed of 1/10000 s) equipped with a Nikon Sigma 135 mm lens were positioned around the arena to provide a top view, and inclined left and right side view of the flying insects. The camera positioning resulted in a zone of focus of approximately 12 cm3 in the centre of the arena in which body and wing movement were in view and in focus of all three cameras. To increase contrast and therefore facilitate subsequent tracking, the cameras were back-lit using three infrared light panels, placed behind the bottom and lower side walls of the octagon arena, opposite each camera (Figure 2A). Prior to filming, cameras were calibrated with the direct linear transformation (DLT) technique(Hartley and Sturm, 1995) by digitizing a wand of 2.8 cm long moved throughout the zone of focus of the cameras. To track the wand movement, we trained an artificial neural network using the pose estimation Python library DeepLabCut(Nath et al., 2019). Computation of the DLT coefficients was then done using the MATLAB program easyWand (Theriault et al., 2014).

Quantifying the in-flight wingbeat kinematics and wing morphology of hoverflies. (A) Hoverflies were released in an octagon-shaped flight arena. We recorded stereoscopic high-speed videos of the flying hoverflies using three synchronised high-speed video cameras; from the videos we reconstructed the three-dimensional body and wingbeat kinematics. Infrared light panels positioned as the bottom of the arena enabled high contrast between the flying insect and the background. (B) Conventional wing angles were measured at each time step (t=0.0004 s) in the body reference frame. ϕ, wing stroke angle within the stroke plane; η, wing deviation angle out of the stroke plane; θ, wing rotation angle along the spanwise axis; U, air velocity vector relative to the wing; α, angle-of-attack of the wing. (C) Wing morphological parameters including their definitions: wingspan R, wing surface area S, radial position along the span r, local wing chord c at distance r, and the second-moment-of-area S2.

During the flight experiment, all the individuals from a single species were released together in the arena to increase the probability that an individual flew through the intersecting fields of view of the three cameras. For each species, the filming experiment was carried out until a minimum of five flight sequences were obtained, whenever possible. Although this procedure facilitated the recording, this prevented identifying which specific individual or sex was recorded flying, leaving open the possibility of recording the same individual multiple times. Note that this approach aimed at estimating the average wingbeat kinematics per species as we could not directly assign flight measurements at the individual level.

Wingbeats selection

First, we determined the flight speed throughout each recorded flight sequence. We did this by tracking the hoverfly body position throughout the flight sequence using a trained DeepLabCut neural network (Nath et al., 2019). From this we calculated flight speed using a central differentiation scheme. Based on these speed data, we selected the flight sequences with the lowest average flight speed, and then we selected within each of these sequences the wingbeat performed during the slowest part of the flight trajectory. For this selection process, we only included the flight trajectory sections in which the hoverfly was flying straight, to prevent possible variation in wingbeat kinematics stemming from flight manoeuvres. We digitised a minimum of three wingbeats per species, each taken from a different flight sequence. For six of the species we tracked three wingbeats, and for two species (Melanostoma mellinum and Sphaerophoria scripta) we tracked six wingbeats.

Quantifying wingbeat kinematics

We quantified the wing movement using the manual stereoscopic video tracker Kine in MATLAB (Fontaine et al., 2009; Fry et al., 2003) originally designed to track the wingbeat kinematics of fruit flies. The tracker was tailored to hoverflies by building wing models matching the wings of the studied hoverfly species. Wing models were obtained by digitizing the wing outline on high-resolution photographs and consisted of 41 2D coordinates representing the wing shape, with the known hinge and tip positions. To reduce complexity, wings were considered as rigid flat plates, although real hoverflies wings endure significant deformation during the wing stroke (Walker et al., 2010).

We tracked the wingbeat kinematics on down-sampled videos (from 5000 to 2500 frames s−1), as this preserved high enough temporal resolution to capture the wing movement (n=27.9±5.3 frames per wingbeat). The tracking provided us with the position and orientation of the body and wings in each consecutive video frame. The position and orientation of the body were defined in the world reference frame, as a 3D position vector (X=[x,y,z]) and the consecutive Euler angles body yaw, pitch and roll. The angular orientation of the wings was expressed in the hoverfly body coordinate system(Muijres et al., 2014) (Figure 2B), as Euler angles relative to the stroke plane: the stroke angle (ϕ), deviation angle (η) and rotation angle (θ) of each wing (Figure 2B). The stroke plane of each wing was defined as the plane at a 45° angle to the long axis of the body passing through the hinge position of that wing. The three angles (stroke, deviation and rotation angle) were measured separately on the left and right wing and then averaged between the two wings. The time history of the averaged angles was then fitted with a 4th order Fourier series for all three angles (Muijres et al., 2014), and their temporal derivatives were computed analytically (i.e. stroke, deviation and rotation rate). We estimated the absolute angular velocity of the beating wing where and are the stroke rate and deviation rate, respectively. The angle of attack (α) was computed as the angle between the wing plane and the velocity vector of the wing (Figure 2B).

Based on the temporal dynamic of the wing angles, we derived the following parameters: wingbeat frequency (f=1/Δt, where Δt is the duration of the wingbeat), stroke, deviation, and rotation amplitude (defined as Aϕ=|ϕmax −ϕmin|, Aη=max −ηmin|, and Aθ =|θmax −θmin|, respectively), and both mean and peak wing angular velocity ( and ωpeak, respectively). We quantified body kinematics during each wingbeat using the mean flight speed, climb angle, and the pitch angle of both body and stroke-plane. Flight speed U was estimated as the temporal derivative of body positions, and climb angle as γclimb=atan(Uver/Uhor), where Uver and Uhor are the vertical and horizontal component of the flight velocity. Body pitch angle βbody was calculated as the angle between the body axis and the horizontal. The equivalent stroke-plane pitch angle was defined as the pitch angle of the stroke plane relative to the horizontal, calculated as βstroke-planebody–45°.

To verify whether the studied wingbeats were representative of the overall recorded flight behaviour, we calculated the mean wingbeat frequency over the full flight sequence as , where nwingbeats and Tsequence are the number of wingbeats in a sequence and the sequence duration, respectively. This mean wingbeat frequency was then compared with that of the studied wingbeats. To estimate whether the analysed wingbeats were representative for hovering flight we calculated the advance ratio of all wingbeats as , where is the wingbeat-average angular velocity of the beating wing and R is the span of a single wing. Hovering is generally defined as flight with advance ratios J<0.1, i.e. when the forward flight speed is less than 10% of the wingbeat-induced speed of the wingtip (C. Ellington, 1984a).

Finally, although hoverfly flight speed was minimised in the wingbeat selection procedure, remaining variation in body kinematics may still affect the wingbeat kinematics and therefore confound with changes in wingbeat kinematics resulting from scaling effect. We therefore assessed if body kinematics affected wingbeat kinematics by performing multiple regression analyses where each wingbeat kinematic metric was treated as the dependent variable and flight speed and climb angle as model effects.

Quantifying morphology

After being filmed, individuals were sacrificed by freezing them at -20°C, and weighed with a resolution of 0.1 mg using an analytical balance (XSR204 Mettler Toledo). Wings were then detached from the body and photographed using a digital microscope (Zeiss Stemi SV11). From these photographs, we measured the wing area S, wingspan R and weight-normalised wingspan R*, mean wing chord , second-moment-of-area S2, and the normalised second-moment-of-area using WingImageProcessor in MATLAB (available at https://biomech.web.unc.edu/wing-image-analysis/). Weight-normalised wingspan is defined as R*=R/m1/3, and is a parameter describing isometric variations in wingspan. The normalised second-moment-of-area is defined as , and is a non-dimensional parameter defining wing shape changes independent of changes in wing size. Both parameters allowed us to compare variations in wingspan and wing shape between species of different sizes, respectively.

Second, we used a landmark-based geometric morphometric method to quantify more precisely the variation in wing shape between species (Bookstein, 1997). Wing shape outline was determined using 300 semi-landmarks equidistantly spaced along the wing outline. Semi-landmarks are commonly used to describe curvy shapes such as wing outlines, which lack typical identifiable landmarks (Gunz and Mitteroecker, 2013). Semi-landmarks are allowed to slide along the curve to establish geometric correspondence (Gunz and Mitteroecker, 2013). After this alignment, the landmarks are treated as regular ones in subsequent analyses. We placed one fixed landmark at the base of the wing, fixing the overall landmark configuration with respect to this homologous position available for all specimens. All landmarks were digitised using TpsDig2 (Rohlf, 2015). Wing outlines were superimposed by performing a generalised procrustes analysis (Rohlf and Slice, 1990) using the geomorph R package (Adams and Otárola-Castillo, 2013). This procedure isolates the shape information from the extraneous variations, namely: the size, position, and orientation. To visualize wing shape variation, we then performed a principal component analysis (PCA) on the superimposed coordinates. Shape changes carried on the PCA axes were visualised using the plotRefToTarget function in geomorph.

Testing the influence of phylogeny on flight kinematics and morphological traits

To assess whether shared evolutionary history among the studied species influenced variation in morphology and wingbeat kinematics, we tested if closely related species showed more similar values in the measured parameters than distantly related ones. This was done by computing phylogenetic signals with the Blomberg’s K-statistics (Blomberg et al., 2003) on each flight and morphological parameters using the phytools R package (Revell, 2012). Tests were performed on the mean values per species using the most recent phylogeny of hoverflies obtained from (Wong et al., 2023), pruned to match our sample of eight species. Because none of the parameters showed significant phylogenetic signal (i.e. variation in traits was unaffected by shared ancestry in our species sample), all consecutive analyses were conducted without phylogenetic correction. Whenever possible, we thus performed analyses at the individual level instead of using the mean per species.

To assess the influence of phylogeny on the evolution of body size in hoverflies more generally, we also tested the phylogenetic signal on a larger number of species using data gathered from literature. This was done using mean values of thorax width (a proxy of body size in insects) of 29 hoverfly species. Here, the morphological data was obtained from (Gilbert, 1985), and the phylogeny was obtained from (Wong et al., 2023).

To gain insight into whether hoverfly species most likely underwent size reduction or enlargement over the course of their diversification, we visualised evolutionary changes in body size through the phylogeny using the contmap function from the R package phytools (Revell, 2012). The contmap function maps a continuous trait (here body size) onto the phylogeny by estimating the ancestral states at the internal nodes using maximum likelihood and interpolating the states along each edge using equation 2 of (Felsenstein, 1985) (see (Revell, 2013)).

Modelling the aerodynamics of hoverfly flight

To estimate how changes in wing morphology and wingbeat kinematics affect the aerodynamic forces required for flight, we modelled the aerodynamic force production by flapping wings using a simplified quasi-steady approach as (C. Ellington, 1984b)

where F is the aerodynamic thrust force produced during a wingbeat. This aerodynamic force thus varies quadratically with the angular velocity of the beating wing (ω2), and linearly with air density (ρ), the second-moment-of-area of the wing (S2), the wing angle-of-attack (α) and the angle-of-attack-specific thrust coefficient (CFα).

The second-moment-of-area S2 is the only morphological parameter in this equation, thus capturing how wing morphology affects aerodynamic force production, from a first principle. But this second-moment-of-area parameter combines both wing size and shape characteristics (i.e., wing area and its distribution over the spanwise axis (C. Ellington, 1984c)). To disentangle the respective effect of wing shape and size on aerodynamic force production, the second moment-of-area can be decomposed as , where R is the wingspan, is the mean wing chord, and is the normalised second-moment-of-area of the wing. Here, wingspan and average wing chord are size dependent parameters; the non-dimensional second-moment-of-area describes wing shape as the distribution of area along the spanwise wing axis, independently of wing size.

Both angular velocity (ω) and the angle-of-attack (α) of the beating wing are wingbeat kinematics parameters that directly affect aerodynamic force production. A flying insect can vary the angular velocity of its beating wings by changing both its wingbeat frequency and amplitude, as angular velocity scales with the product of wingstroke amplitude and wingbeat frequency (ω∼Aϕ f).

The decomposed model for estimating the effect of wing morphology and kinematics of aerodynamic force production in flapping insect flight is thus

Physical scaling of kinematics and morphology parameters with size

In hovering flight, this aerodynamic thrust force should be directed vertically upwards and equal to the weight of the animal. Thus, to maintain weight support among different sizes of hoverflies, the aerodynamic force should scale linearly with body mass. This cannot be achieved through isometric scaling, i.e. the expected scaling if geometrical similarity is preserved.

For isometric scaling, the morphological parameters of interest scale with body mass as R3m1/3, and , and thus both S2 and F scale with body mass as S2Fm4/3. This value is larger than one, and so with an isometric reduction in size the aerodynamic force for weight support reduces more quickly than body mass. Thus, compared to large species of hoverflies, small species need relatively larger wings or adjust their wingbeat kinematics to maintain weight support.

Using our aerodynamic force model (Eqns 1-2) we can estimate both (1) how wing morphology should scale with size if wingbeat kinematics are kept constant, and (2) how the wingbeat kinematics should scale if wing morphology scales isometrically. (1) To maintain weight support with non-varying wingbeat kinematics, the second-moment of-area should scale linearly with body mass, resulting in negative allometry. This could be achieved by negative allometric scaling in the wingspan, chord and normalised second-moment-of-area (R<m1/3, and ), or a combination of these. (2) To maintain weight support when wing morphology scales isometrically, wing angular velocity should scale with body mass as ω∼m1/6. Again, a hovering hoverfly could achieve this using a comparative adjustment in wingstroke amplitude or wingbeat frequency, or a combination of both. Furthermore, the animal could adjust the aerodynamic force vector using an equivalent adjustment in the wing angle-of-attack α.

To examine how wing morphology and wingbeat kinematics scales with body mass, and whether this scaling is isometrically or allometrically, we regressed body mass (m) against the morphology and kinematics metrics in Eqns 1-2 (, R3, ω2, f, Aϕ and α) using reduced-major-axis (RMA) regressions. RMA is commonly used in allometric studies as it assumes a similar level of error in both variables(Aiello, 1992; Smith, 2009). The symmetric error structure in both variables assumed in RMA makes this method more suited to summarize the relationship between two variables without assuming a causal direction(Warton et al., 2012). In the regression analysis, we used the averaged value at mid-stroke of ω and α, where aerodynamic force production is close to maximum (Tian et al., 2013). Measurements were all log10-transformed prior to regression analysis.

For each morphological metric, we tested if the observed slope significantly deviated from the isometric expectation (expected scaling if geometrical similarity is preserved: Rm1/3; ) using the sma function in the smatr R package (Warton et al., 2012). To better visualize variations in wing shape and size relative to body mass among species, we scaled all wings with the length-scale for isometry (liso=m1/3), and calculated the corresponding weight-normalised wingspan as R* = R/m1/3. We performed additional RMA regressions to test if other kinematic adjustments were associated with changes in body mass, in addition to the ones considered in the aerodynamics model (Eqn 2). This includes regressions of body mass versus the amplitudes of the measured wing angles (Aϕ, Aη, Aθ), and the corresponding peak angular rates .

Computational fluid dynamics simulations

We performed Computational Fluid Dynamics (CFD) simulations to quantify aerodynamic force produced by the different flying hoverfly species. These simulations allowed us to explicitly quantify the effect of size, wingbeat frequency and morphology on aerodynamic force production. Furthermore, we used the CFD results to test the validity of our simplified quasi-steady aerodynamic model (Eqns 1-2).

Each simulation was performed with a single rigid, flat wing with species-specific contours, and the averaged wingbeat kinematics across all species. The computational domain was 8R×8R×8R, and we simulated three consecutive wingbeat cycles. The wing has a thickness of 0.0417R. We assumed hovering flight mode, so no mean inflow was imposed (U = 0 m s-1).

We performed the simulations using our in-house open-source code WABBIT (Wavelet Adaptive Block-Based solver for Insects in Turbulence)(Engels et al., 2022, 2021). The computational approach is based on explicit finite differences that solve the incompressible Navier–Stokes equation, in an artificial compressibility formulation. Wavelets are used to dynamically adapt the discretisation grid to the actual flow at every given time t. This allows us to focus the computational effort where it is required to ensure a given precision, while the grid is significantly coarsened wherever possible. Starting from a coarse base grid, we allow up to seven refinement levels, each one refining the grid by a factor of two. In this work we use the Cohen-Daubechies-Feauveau wavelet CDF42. The grid is composed of blocks of identical size with Bs=22 points in all three spatial directions, yielding an effective maximum resolution of 352 points / R. An additional simulation with eight refinement levels confirmed validity of our results. As demonstrated in numerous validation tests, obtained numerical data can be regarded as quasi-exact in the sense that they are equivalent to an experiment with identical geometry and kinematics. More details on the code can be found in(Engels et al., 2021); the required parameter files to reproduces the simulations presented here can be found in the supplementary. To accelerate the computation, each simulation used up to 600 CPU cores with 2.7 TB of memory in total.

From a fluid dynamics point of view, changing wing length R or wingbeat frequency f solely changes the Reynolds number Re (Batchelor, 1967). For each wing contour, we simulated two different values for Re, one corresponding to a fixed averaged frequency and wing length for all species, and a second one where the wingspan was species-specific. This was initially done to separately assess the combined effect of Reynolds number and wing shape from that of wing shape only. The non-dimensional aerodynamic force was however very similar in both cases, hence the influence of Reynolds number is negligible. In general, the range of Reynolds numbers for the studied hoverfly species is rather narrow (Re ∈ [324,1310]). We finally assessed the contribution of interspecific variation in wingbeat frequency to force production by rescaling all aerodynamic forces with the square of the species-specific and average wingbeat frequency ratio .

Assessing the relative contribution of wing size, shape and kinematics changes to aerodynamic force production

By combining the CFD simulation results with our physical scaling law analysis, we assessed the relative contribution of changes in wing size, shape and kinematics between the studied species on the aerodynamic force production required for weight support.

Because wingbeat kinematics was found to vary little across species and to be uncorrelated with body mass, we based this analysis on the average wingbeat kinematics of all species combined. We did assess the contribution of interspecific variation in wingbeat frequency on aerodynamic force production for weight support, as this varied most between species. For this analysis, we used Eqn 3 to reconstruct the aerodynamic forces produced by each wing type (species), based on both the species-specific wingbeat frequency as well as for the average wingbeat frequency across all species .

We assessed the relative contribution of changes in wing morphology to weight support among species using the following model. According to our quasi-steady aerodynamic model, vertical aerodynamic force for weight support scales with the second-moment-of-area of the wing, which is decomposed into the wing size and shape parameters wing span, mean chord and normalised second-moment-of-area as .

Because isometric scaling results in a reduction in vertical force to weight ratio with decreasing size, smaller insects would require wings with relatively high second-moment-of-area. We estimated how hoverflies achieve this by estimating the scaling between the wing morphology parameters and the vertical aerodynamic force produced by the beating wings, estimated using the CFD simulations. The resulting scaling factors show how variations in wingspan, chord length and normalised second-moment-of-area contribute to producing weight support across sizes. By comparing the scaling factors with those for isometric scaling, we determined how allometric adaptations of wing size and shape alter these aspects among hoverflies of different body sizes.

Results

In this study, we investigated the morphology, and hovering flight biomechanics and aerodynamics of eight hoverfly species, varying in body mass from 5 to 100 mg. The sampled species were relatively well spread throughout the phylogeny of hoverflies (Figure 1). Variation in body mass across species was however unevenly distributed: only one species, Eristalis tenax, weighed more than 50 mg, while the two closely-related species Syrphus nigrilinearus and Eupeodes nielseni had highly similar body masses (Figure 1). Number of males and females from the sampled species was very unbalanced (Table S1), limiting our capacity to investigate for difference between sexes within species. Sexual dimorphism should nonetheless be negligible compared to differences between species.

Influence of phylogeny on the evolution of wingbeat kinematics, wing morphology and body size

There was no significant effect of phylogeny on the flight and morphological parameters among the studied species (n=8) (Table S1). When testing the phylogenetic signal on body size on a larger number of species (n=29), we however detected a marked effect of phylogeny (P=0.005; Blomberg’s K=1.5). The mapping of ancestral state estimates onto the phylogeny suggested that body size has mostly decreased over the course of the diversification of hoverfly species (Figure S1). Although not visible among the eight studied species, variation in body size in hoverflies appears structured by phylogeny overall. Larger sizes are found in the more basal species (e.g. genus Volucella) whereas smaller sizes is found in more recently diverged species (e.g. genus Platycheirus).

No correlation between wingbeats kinematics and body mass

We quantified a total of 33 wingbeats across eight hoverfly species, including at least three wingbeats per species

(Figure 3). The wingbeat frequencies of the studied wingbeats were highly correlated with the average wingbeat frequency calculated on the entire flight sequence (; n=33 flight sequences; r=0.98, P<0.001). Here, the number of wingbeats per flight sequence was 26.2±28.4. This suggests that the digitised wingbeats reliably reflect those of the full flight sequences. Furthermore, the mean advance ratio for the digitised wingbeats was below the generally accepted threshold of 0.1, confirming that the hoverflies were operating in a hovering flight mode.

Relationship between wingbeat kinematics and body mass among hoverfly species. (A-E) Temporal dynamics of the wingbeat kinematics throughout the cycle of all digitised wingbeats. Separate wingbeats are colour coded by body mass (see legend on top), and the black lines show the average wingbeat kinematics for all wingbeats combined. (A-C) Temporal dynamics of the three conventional wingbeat kinematics angle (stroke, deviation, and rotation angle, respectively). (D-E) Angular speed and angle-of-attack of the wings throughout the wingbeat cycle, derived from the stroke and rotation angle, respectively. (F-I) Derived Wingbeat kinematics parameters for each studied species versus the body mass of that species. The kinematics parameters are mean angular wing speed, stroke amplitude, angle-of-attack at mid stroke, and wingbeat frequency. Each data point shows mean ± standard error. None of the wingbeat kinematic parameters were significantly associated with body mass.

The eight hoverfly species exhibited similar wingbeat kinematics (Figure 3A-E) whereby they all moved their wings forward and backwards sinusoidally within the stroke plane (Figure 3A), with a wingstroke amplitude of Aϕ=100°±18° (Figure 3G). Movements out of the stroke-plane, as expressed by the deviation angle, are minimal (Figure 3B). Therefore, the oscillatory angular speed of the wingbeats (Figure 3D) are primarily caused by the wing stroke movement, and peak at approximately 5°×104 s-1 (Figure 3F). During the majority of each wingbeat (forward and backward wingstroke), the wing operates at approximately constant wing rotation angle and angle-of-attack (Figs 3C and 3E, respectively). At the end of the forward and backward wingstroke, the wings rapidly supinate and pronate, respectively. This causes the wing to flip upside down at stroke reversal, allowing it to operate again at the constant wing angle-of-attack of (Figure 3H). The hoverflies flew with a mean body pitch angle of βbody=30°±16°. The corresponding pitch angle of the wing stroke-plane was βstroke-plane=−15°±16°, showing these hoverflies hover with an inclined stroke-plane.

The temporal dynamics of wingbeat kinematics was not visually associated with variation in body mass (Figure 3A-E). Accordingly, none of the derived wingbeat kinematics parameters varied significantly with body mass (Table 1, Figure 3F-I). These include the parameters that directly affect aerodynamic force production (Eqn 1-2, angular velocity ω, angle-of-attack α, stroke amplitude Aϕ and wingbeat frequency f), and the additionally estimated wingbeat kinematic parameters (rotation amplitude Aθ, deviation amplitude Aη, peak rotation velocity and peak deviation velocity ) (Table 1).

Results from reduced major axis regression of hovering-flight wingbeat kinematic parameters against body mass. Test were performed using mean values per species (n=8). See figures 2 and 3F-I for definitions of primary kinematics parameter and visualisations of the data, respectively.

Furthermore, variations in wingbeat kinematics were little affected by the hoverfly flight kinematics. Only flight speed had a positive significant effect on the angular velocity of the beating wings (Table S2). Thus, faster flying hoverflies beat their wings at higher angular velocity. When controlling for this effect of flight on wing angular velocity, the regression between angular wing speed and body mass remained non-significant (P=0.41).

These combined results suggest that variation in wingbeat kinematics between species cannot be explained by the demands of maintaining weight support throughout (isometric) size changes. Adjustments for weight support among hoverfly species thus more likely occur through changes in wing morphology.

Wing morphology covaries with body mass

Wing morphology strongly covaried with body mass: the second moment of wing area (S2), wingspan (R), wing chord , and the non-dimensional second moment of wing area all significantly correlated with body mass (Figure 4). Covariation between the wing morphological parameters and body mass all deviated from the isometric expectation (Table 2). Second moment of wing area, wingspan and mean chord length showed all a slight but significant negative allometry (P=0.002, P=0.015 and P=0.014, respectively), and were thus disproportionately larger in smaller species (Figure 4A-C, Table 2). The deviated more strikingly from the isometric trend (P<0.001, Figure 4D, Table 2), indicating a marked change in wing shape associated with body mass variation. Note that among the smallest species the negative correlation between the and body mass reduces (Figs 4D). In fact, for the species weighing less than 20 mg normalised second-moment-of-area remains mostly constant (correlation between and body mass among species weighing less than 20 mg: r=0.34; P=0.157). This suggests an upper limit to the increase in normalised second-moment-of-area with reducing size.

Results from reduced major axis regression of wing morphology parameters against body mass (n=44 individuals among eight species). See figures 2 and 4 for definitions of primary morphology parameter and visualizations of the data and regressions, respectively.

Wing morphology parameters versus body mass for all studied species, including isometric scaling and the observed reduced major axis regression. (A-D) body weight (abscissa) versus on the ordinate the second-moment-of-area, wingspan, mean chord and normalised second-moment-of-area, respectively. Each data point shows results for a different individual, and is color-coded according to species (see top of A). The expected slopes for isometric scaling are indicated by dashed grey lines, and the best fitting reduced major axis regressions are indicated by continuous black lines (see top of B). All fitted regression slopes are significantly lower than expected under isometry, suggesting negative allometric scaling of wing size and shape with respect to body mass.

Our geometric morphometrics analysis revealed two main axes of wing shape variation among the studied species associated with the principal components PC1 and PC2 that explain 65.24% and 19.07% of the variation, respectively) (Figure S2). Shape variation carried on PC1 was strongly correlated with the as well as body mass (r=−0.78; P<0.001) (Figure 5), depicting the changes in wing shape associated with body mass variation among hoverflies species. This change in reflects a wing surface area mostly located distally in smaller species, whereas most wing surface area is located near the wing base in larger hoverflies species (Figure 5).

Changes in wing shape and size associated with body mass variation and according the geometric morphometrics analysis. (A) normalised second-moment-of-area (S2) vs. the primary Principal Component (PC1) of the geometric morphometrics analysis on the wing outlines (left), and weight-normalised wing outlines (with coordinates X*=X/m1/3) color-coded with body mass (see legend) (right). (B) body mass (m) versus the primary Principal Component (PC1) (left), and weight-normalised wing outlines and wingspan (R*) versus body mass (right). Small transparent data points and large opaque data points show values per individual and per species, respectively. All data points and wing outlines in B are color-coded according to species (see right of B). (A) Variation in the normalised second-moment-of-area (S2) were captured well by PC1 (see also Fig. S2). (B) In larger species, wing surface area was located more proximally (lower values on PC 1) than in smaller species, in which wing area tend to be located more distally (higher value on PC 1). The left and right gritted wing shapes on the PC1 axis show the hypothetical wings for PC1=-0.08 and PC1=0.06, respectively. PC1 explains 65.24% of the variations in wing shape (see also figure. S2). (A-B) Combined with the changes in wing shape (left), weight-normalised wing size and wingspan (R*) tend to be larger in smaller species (right).

The combined correlative and morphometrics analyses thus shows that as their size decreases, hoverflies exhibit wings with both larger normalised second-moment-of-area, and relatively larger wings as depicted by the larger weight-normalised wingspan (Figure 5).

Aerodynamic modelling confirms the main role of wing morphology for weight support

As among the eight hoverfly species, wingbeat kinematics did not vary significantly with body mass (P=0.46), we performed all CFD simulations with the average wingbeat kinematics of all species combined (Figure 6A). We did vary wingbeat frequency (Figure 6B), and the wing shape and size between species.

Aerodynamic forces produced by hovering hoverflies, as estimated using Computational Fluid Dynamic (CFD) simulations. (A-B) For our simulations, we used the species-specific wing shapes and sizes (figure 5), and average wingbeat kinematics (A) and the mean wingbeat frequency (B) across all eight studied hoverfly species. (C) The resulting temporal dynamic of vertical forces throughout the wingbeat cycle, coloured by species (see legend above B). (D) The wingbeat-average vertical force versus body mass, for all simulated hoverfly species operating at both the average wingbeat frequency for all species and their species-specific frequency (square and round data points, respectively). Trendlines for weight support (F=mg), and results of the average wingbeat frequency and species-specific frequency simulations are shown as grey dashed, black dashed and black solid lines, respectively. With the species-specific wingbeat frequency (shown in panel B) hoverflies are closer to producing weight support than for the average wingbeat frequency simulations.

The resulting aerodynamic forces differed between species mostly in magnitude, but the temporal dynamics of force production was strikingly similar (Figure 6C and S3). For all species, the temporal dynamics of the vertical aerodynamic force showed a large force peak preceded by a small peak, during each wing stroke (forward and backward stroke). The back stroke produces a larger vertical force and a more distinct higher harmonic. These dynamics are similar between species, despite their different wing shapes. This suggests that wing shape does not strongly affect the characteristics of aerodynamic force production.

The wingbeat-average vertical aerodynamic forces showed a strong scaling with body mass (Figure 6D). Hence, similar wingbeat kinematics but different wing morphology enable achieving weight support in the studied hoverflies species. For each studied species, we estimated aerodynamic forces with wings flapping at the mean wingbeat frequency of all species, and using their species-specific wingbeat frequency (Figure 6B,D). All species were closer to producing weight support with their species-specific wingbeat frequency than for the equivalent average wingbeat frequency case (Figure 6D). Variation in wingbeat frequency thus appeared necessary for individual hoverflies to achieve weight support, although between species wingbeat frequencies do not correlate with body mass.

The correlation between vertical forces estimated using CFD and the different wing morphology parameters (Figure 7A-D; S2, R, and) were strikingly similar to the comparative correlation between body weight and these morphology parameters (Figure 4). This shows that our simplified aerodynamic model (Eqns 1-2) captures the major effects of wing morphology on aerodynamic force production.

Relationships between wing morphology and vertical force production during simulated hovering flight. (A-D) Wingbeat-average vertical force (abscissa) vs. on the ordinate the dimensional and normalised second-moment-of-area, wingspan, and mean wing chord, respectively. Data points are results for different species, as color-coded according to the legend above panel E). The expected slopes for isometric scaling are indicated by dashed grey lines, and the best fitting reduced major axis regressions are indicated by continuous black lines (see top of A). All fitted regression slopes are significantly lower than expected under isometry, suggesting negative allometric scaling. (E) The relative contributions of different aspects of wing morphology (R, , and ) to total aerodynamic force production (100% for ) in the case of isometric scaling and for the observed allometric scaling (in grey and black, respectively).

We continued to assess the relative contribution of the different wing morphology parameters (R, and ) to vertical force production (Figure 7E). This shows that for both the isometric scaling and the observed allometric scaling, variations in wingspan affect the changes in vertical force the most (75% and 74%, respectively). In contrast, the allometric changes in , and especially in, provided most of the allometric change in vertical force production compared to the isometric case. In fact, for our observed allometric case, the changes in mean chord contributed relatively less to vertical force production than for the isometric case (a reduction of 9%). This corresponds to the 11% increase in relative contribution of allometric changes in to the associated aerodynamic forces (Figure 7E). This shows that allometric changes in wing shape, rather than wing size, between species of different sizes have the largest relative effect on achieving weight support.

Discussion

Allometric changes in wing morphology enables miniaturised hoverflies to maintain weight support during hovering

In this study, we investigated whether variation in body mass among hoverfly species is associated with changes in wingbeat kinematics or in wing morphology. Both components directly influence the aerodynamic forces generated during flapping flight (C. Ellington, 1984b), and therefore could accommodate for weight support during evolutionary changes in size occurring throughout the diversification of hoverfly species. Comparing how wingbeat kinematics and wing morphology scale with body mass across hoverfly species, we showed that wing size and shape rather than wingbeat kinematics scale with body mass. This suggest that wing morphology might respond stronger than wingbeat kinematics to the constraints imposed by weight support in hoverflies.

Changes in wing morphology associated with body mass variation are commonly observed in flying animals (Le Roy et al., 2019; Norberg, 2007; Rayner, 1988; Tercel et al., 2018), indicating that wing and body morphology evolve together as an integrated phenotype. Increasing size is usually not only associated with larger wings, but also accompanied by changes in wing shape (Danforth, 1989; Outomuro et al., 2013). Teasing apart the changes in wing morphology resulting from constraints imposed by body mass from that of other selective pressures is however not straightforward. It is generally facilitated by comparing species with similar ecology, because this limits the effect of contrasted ecological pressures at play in different species (García and Sarmiento, 2012). Hoverfly species present a fairly similar ecology in that all adults are specialised flower-visitors (Klecka et al., 2018). Comparable selective regime are thus likely acting on the evolution of flight performance and associated flight kinematics and morphology in adult hoverflies.

Shared ancestry can also affect the evolution of wing morphology and its effect should thus be controlled to understand how body mass may constrain the evolution of other traits (Blomberg et al., 2003). Among the eight species studied here, variation in the measured morphological traits was not explained by phylogenetic distance. The observed changes in wing morphology accompanying body mass variation can thus be assumed free from the influence of phylogeny. This situation however does not apply to all hoverfly species. Indeed, the test of phylogenetic signal on an enlarged number of species revealed that phylogeny explain body mass variation in hoverflies more broadly (Figure S2). Changes in wing morphology resulting from constraints imposed by body mass are therefore expected to be confounded with the effect of phylogeny in a sample including more species. The association between body mass and phylogenetic history was further highlighted by the ancestral state estimation of body size performed on 29 species. This showed that most recently diverged hoverfly species tend to be smaller, and most likely derive from an ancestor of larger size. Throughout the diversification of hoverflies, evolutionary changes in size thus most plausibly showed a decreasing trend, resulting in evolutionary miniaturisation. This is in line with the general trend observed in insect more generally (Polilov, 2015). Hence, our findings are better discussed in the perspective of decreasing size.

As size decreases, wings become relatively less effective in producing the aerodynamic forces required for weight support (Schmidt-Nielsen, 1975). The constraints imposed by weight support during hovering flight should then promote the evolution of disproportionately large wings, or an increased wingbeat frequency or amplitude in smaller hoverflies species. We found that none of wingbeat kinematic parameters scaled with body mass but that wings indeed tended to be disproportionately larger in smaller species (i.e. negative allometry).

Here, all the wing morphology parameters associated with aerodynamic force production exhibited significant negative allometry, including wingspan R, mean chord , and both dimensional and dimensionless second-moment-of-areas (S2 and , respectively). Changes in multiple aspects of the wing morphology may hence contribute to mitigating the negative effect of size reduction on flight ability. A negative allometry between wing and body size has been found in other insect groups: for solitary bees this was found at both the interspecific and intraspecific level (Grula et al., 2021), and at the intraspecific level in rhinoceros beetles (Kawano, 1995). Disproportionately larger wing in larger species (positive allometry) was however found in dragonflies (May, 1981). Smaller insects may not systematically show disproportionately larger wings because of the diversity of mechanisms mitigating weight support involving either morphological or kinematic changes, but also because of contrasted flight modes across taxa (e.g., gliding flight in dragonflies).

The covariation between wing shape and body mass is most apparent among the hoverfly species with body mass above 20 mg, as for species with body masses below this value the normalised second-moment-of-area of the wing remains almost constant (Figure 4D; i.e. Platycheirus clypeatus, Melanostoma mellinum and Tropidia scita in our study). The flattening of the to body mass curve with decreasing size suggests an upper limit in the normalised second-moment-of-area among hoverfly species of approximately 0.58 (Figure 4D). Thus, for small hoverfly species the redistribution of wing area towards the wingtip might be limited, potentially due to structural constraint. The bending moment on the wing hinge and consequently the flight power muscles scale with the third-moment-of-area of the wing , which scales with the product of and R (Muijres et al., 2017). Thus, an increase in comes at a stringent flight muscle force requirement; this might explain the observed plateau in for small hoverfly species. More research is needed to quantify these apparent conflicting evolutionary driving forces on wing shape.

Here we further investigated the relative contribution of the allometric changes in the different aspect of wing morphology to force production by decomposing analytically the allometric component of R, , and (Figure 7). This highlighted that allometric changes in resulted in the largest additional force production compared to isometric scaling. This change in wing shape that most strikingly deviated from isometry was emphasised in our geometric morphometrics analysis. Indeed, in our Principal Component Analysis (PCA) performed on superimposed wing shape coordinates, the primary axis PC1 explained a striking 65% of the variations in wing shape, and was strongly correlated with normalised second-moment-of-area . This PC1 shape variation depicts a shift in the spanwise distribution of wing area, from more proximally located in larger species to more distally located in smaller species (Figure 5). Interestingly, similar changes in wing shape associated with decreasing body size have been found among species of Hymenoptera (Danforth, 1989). An illustration of the most pronounced end of this pattern is the highly spatulated wing shape found in tiny wasp species (e.g. in the Mymaridae or Mymarommatidae families). As highlighted by our combined quasi-steady model and CFD results, such wing shape confers increased lift force production during flapping flight (C. Ellington, 1984b), and hence also contribute to maintaining weight support with decreasing body size. The reduction in body size occurring through the diversification of hoverfly species may thus have favoured wing shapes in smaller species that produce aerodynamic forces more effectively.

Factors decoupling wingbeat kinematics from the influence of body mass

Our study revealed that wingbeat kinematics is relatively stable across hoverfly species, with none of the measured kinematics parameters significantly correlated with body mass. Performing computational fluid dynamic simulations using the average wingbeat kinematics for all species, but incorporating differences in wing morphology, demonstrated that weight support was indeed achieved across species. Important however is that variation in wingbeat frequency – albeit uncorrelated with body mass – has proven crucial for correctly estimating the aerodynamic force necessary to achieve weight support. Adjustments in wingbeat frequency thus appear to play a role in achieving weight support during hovering for individual hoverflies. The hovering hoverflies might thus primarily use wingbeat frequency to finely tune the vertical force required for weight support. Unlike birds and bats, flying insects cannot morph actively their wings to control aerodynamic force production (Harvey et al., 2022). As a result, flying insects rely purely on wingbeat kinematics changes to trim aerodynamic forces and torques during steady flight (Muijres et al., 2017), and to adjust these forces and torques during manoeuvring flight (Hedrick et al., 2009; Muijres et al., 2014). Our CFD simulation results suggest that hovering hoverflies use relatively simple wingbeat frequency modulations to compensate for the continuous daily fluctuations in body mass, as a result of food and water intake and excretion.

The lack of covariation between in-flight wing movement and body mass goes against the prediction that wingbeat kinematics is more likely to be adjusted with varying body mass due to the greater flexibility of flight as compared to morphology. This finding is also surprising given that size variation is generally accompanied with changes in wingbeat kinematics in flying animals (Norberg, 2007; Rayner, 1988; Riskin et al., 2010; Tercel et al., 2018). A previous study dissecting the flapping wing kinematics of miniaturised insects found an increased wing deviation in smaller species enabling higher wing velocity and hence aerodynamic force production(Lyu et al., 2019). Such a change in wing kinematics associated with decreasing size can be explained by the increased viscous drag constraint that tiny flying insects experience, requiring higher aerodynamic force generation. Here we did not find an increased wing deviation angle in smaller hoverfly species, most-likely because our study focus on a too narrow range of size.

While only a handful studies have examined the effect of body mass on detailed wingbeat kinematics in insects (such as wingbeat amplitudes and angle-of-attack Belyaev et al., 2014; Lyu et al., 2019; Meresman and Ribak, 2017), the scaling between body mass and wingbeat frequency has been frequently assessed (Bartholomew and Casey, 1978; Byrne et al., 1988; Casey et al., 1985; De Nadai et al., 2021; Dudley, 1990; Tercel et al., 2018). In theory, an increase in wingbeat frequency is expected in smaller insects because flapping wing-based aerodynamic force production decreases faster than body mass. As such, it is puzzling that the reduction in size among our studied hoverfly species is not associated with increased wingbeat frequency.

Interestingly, inconsistent results on the scaling between body mass and wingbeat frequency have been found among other groups of flying animals. Although predominantly negative (e.g. in butterflies (Dudley, 1990), in orchid bees (Darveau et al., 2005), in mosquitoes(De Nadai et al., 2021), and in a comparative study across insect orders (Tercel et al., 2018)), several studies found weak or no relationship between body mass and wingbeat frequency (e.g. in moths (Bartholomew and Casey, 1978), whiteflies (Byrne et al., 1988), bees (Duell et al., 2022; Grula et al., 2021)). Consistent with our findings, a study comparing the wingbeat kinematics during tethered flight among hoverfly species spanning a range of body mass similar to that investigated here found no significant effect of body mass on wingbeat kinematics (Belyaev et al., 2014).

A possible explanation to the heterogeneous trends found among insect groups is that the extent to which body mass constrains the evolution of wingbeat kinematics depends on species ecology. Adaptation to a particular ecological niche may drive the evolution of specialised flight abilities, promoting a specific wingbeat kinematics. Selective pressures exerted on these wingbeat kinematics may sometimes overcome the constraints imposed by weight support. In that case, relatively constant wingbeat kinematics would be maintained over a range of body mass, and changes in wing morphology alone are expected to adjust for weight support. Consistent with this hypothesis is the decoupling of wingbeat frequency from body mass variation, which has been highlighted in species where hovering flight ability is crucial for food acquisition. These include bees (Duell et al., 2022; Grula et al., 2021), sphingid moth (Bartholomew and Casey, 1978), hummingbirds (Skandalis et al., 2017), and hoverflies (Belyaev et al., 2014) as also in the present study. The absence of correlation between body mass and wingbeat kinematics observed in hoverflies may thus be due to their highly specialised hovering flight abilities.

The selective pressures promoting the evolution of unique hovering flight abilities in hoverflies are not precisely known but hovering steadily is likely advantageous for navigating between flowers when foraging for nectar. Good hovering flight performance may also have been promoted because it increases mating opportunities during the leck gathering of males (Downes, 1969; Gilbert, 1984; Heinrich and Pantle, 1975). Selection for increased foraging efficacy and/or mating opportunities may have maintained consistent wingbeat kinematics across hoverfly species despite variation in body mass.

Conclusion

Altogether, our results suggest that the hovering flight abilities of hoverflies may stem from highly specialised wingbeat kinematics that have been conserved throughout their diversification. In contrast, changes in wing morphology have evolved with the aerodynamic constraints associated with the apparent evolutionary miniaturisation of hoverflies. The existence of selective pressures maintaining stable wingbeat kinematics despite variation in body mass implies that high hovering flight performance can only be achieved over a limited range of wingbeat kinematic parameters, which remains to be ascertained.

Our study moreover highlights that unlike many other Diptera, hoverflies make use of a highly-stereotyped inclined stroke-plane wingbeat kinematics for hovering, where during the dorsal-ventral wingstroke the wing moves not only from back to front but also downwards (Mou et al., 2011). This wingbeat kinematics pattern may contribute to the specialised hovering abilities of hoverflies, as it is consistently observed among the eight species studied here, with their stroke-plane pitch angle during hovering being approximately orientated at −30°. Furthermore, despite significant changes in wing morphology and wingbeat frequency among the studied species, the temporal dynamics of aerodynamic force production is strikingly similar between them (Figures 6C and S3). This suggests that not only the wingbeat kinematics is conserved between species, but also the underlying aerodynamics including the relative use of different unsteady aerodynamic mechanisms (Sane, 2003). Thus, the apparent highly-specialised flight style of hoverflies is maintained among several species with a large range of body masses, and robust against variations in wing shape and wingbeat frequency.

List of symbols and abbreviations

Acknowledgements

The authors would like to thank Ilam Bharathi for helpful discussions about insect flight aerodynamics, and Remco Pieters for his help in building and operating the experimental setup.

Funding

This work was supported by an NWO Vidi research grant to F.T.M. (I/VI.Vidi.193.054). The study was also provided with computer and storage resources to T.E. by GENCI at IDRIS thanks to the grant 2023-A0142A14152 on the supercomputer Jean-Zay’s CSL partition.

Supplementary material

Ancestral state estimates for body size plotted onto a phylogeny of hoverflies. Ancestral states were estimated and plotted along the tree branches following method described in Revell (2013). Thorax width data (a proxy of body size) were obtained from Gilbert (1985) and the phylogeny from Wong et al. (2023). n=29 species.

Result of geometric morphometrics analysis on the wing outlines of eight hoverfly species. Data points show results for different individuals, coloured by species (see legend on top). (A) The two first principal components (PC1 and PC2) of the geometric morphometrics analysis, together with the associated shape changes and percentage of the variation explained by the principal components. Shown wing shapes represent the extreme values associated with each PC axis. (B) PC1 and PC2 (first and second row, respectively) plotted against body mass, highlighting that the change in wing shape carried on PC1 is associated with variations in body mass.

Temporal dynamic of the vertical aerodynamic forces produced during the wingbeat of the eight studied hoverfly species, estimated using Computation Fluid Dynamics (CFD) simulations. Data for the eight species are color-coded according to the legend on top. All simulations were run with the mean wingbeat kinematics of all hoverflies combined (Fig 6A), but with species-specific wing morphology (see legend on top). The aerodynamic forces in the different panels are scaled and normalized using four different methods: (A) aerodynamic forces in mN as directly estimated using the CFD simulations, and where all wings operate at the mean wingbeat frequency of all hoverflies ; (B) aerodynamic forces as produced by each hoverfly beating its wings at the species-specific wingbeat frequency; (C) aerodynamic forces normalized with its wingbeat-average force, ; (D) aerodynamic forces produced using the species-specific wingbeat frequency, and normalized with the mean weight of the specific hoverfly species F/mg.

Number of female and male individuals in the studied species.

Phylogenetic signal computed on morphological and flight traits on the eight studied species.

Results from multiple regressions testing the effect of body dynamics on the wingbeat kinematics parameters.