### Study animals

A first breeding population of *P. philippinicum* was obtained from the Audubon Insectarium in New Orleans, Louisiana, USA and shipped to the University of Montana, Missoula, Montana, USA. The insects were housed in a transparent plastic container (50 × 40 × 60 cm) at 22 °C, on 12 h:12 h light:dark cycles, sprayed with water daily (RH = 50–80%), and fed fresh *Rubus idaeus* leaves ad libitum. This population was used to investigate morphological scaling relationships and conduct the flight experiments described below.

A second culture stock of *P. philippinicum* was obtained from Kirsten Weibert (Jena, Germany) and captive bred in the department of Functional Morphology and Biomechanics at Kiel University, Germany. The specimens were kept in a large glass cage with proper ventilation at 20–22 °C (RH = 50–80%; 16 h:8 h light:dark cycle) and fed with fresh blackberry (*Rubus* sp.) and common oak (*Quercus robur* L.) leaves ad libitum. This population was used to investigate attachment pad morphology and attachment forces.

### Sexual dimorphism and scaling relationships

Photographs of the animals in dorsal view were taken using a DSLR camera (EOS 600D, Canon Inc., Tokyo, Japan). Using ImageJ software (v.1.52k) [67], we measured body length (BL, mm), body area (mm^{2}), body circularity (dimensionless, \(\frac{4\pi \times Area}{{Perimeter}^{2}}\)), body aspect ratio (dimensionless, \(\frac{Body\, length}{Average\, body\, width}\)), mean antenna length (mm), mean front femur length (mm), and total wing area (mm^{2}, including both forewings for females or both hindwings for males) in 25 adult males and 19 adult females (Additional file 1: Fig. S1). Antenna length could not be measured in seven males and three females because they were missing flagellomeres on both antennae. Wet body mass (g, measured immediately after flight trials) was obtained using an analytical balance (ME54TE/00, Mettler Toledo, Columbus, OH, USA). We calculated wing loading (N m^{−2}) as wet body mass multiplied by gravitational acceleration (g = 9.81 m s^{−2}) and divided by total wing area. Male flight muscle mass (mg) was obtained by dissecting the muscles out of the metathorax of freshly dead males (n = 23) drying them at 70 °C for 24 h and weighing them with a more accurate analytical balance (UMT2, Mettler Toledo, Columbus, OH, USA).

### Male flight performance

To evaluate male flight performance, here defined as a righting maneuver that required sustained control over long-axis body angle and climb ability, we dropped adult males in the air and recorded their flight trajectories in 2D (Additional file 2: Video S1, Additional file 3: Video S2). Working in Missoula, MT (elevation 978 m above sea level, average air density = 1.07 kg m^{−3}), adult males (N = 16, 0.43 ± 0.006 g) from our American culture population were held by the thorax and dropped by the experimenter at a horizontal body angle from a constant height above a floor (1.5 m). A 2 m-tall slab of ponderosa pine (*Pinus ponderosa*) with natural bark was placed vertically in front of the animal, 2 m away, to serve as a target and landing site (Additional file 1: Fig. S3). The experimental room was largely featureless with white walls and at 26 °C. The floor was covered with thick blankets to avoid injuries if crashing. We recorded the flight of the insects using a high-speed video camera (Photron FASTCAM SA-3, Photron USA, San Diego, CA, USA) sampling at 500 fps with a shutter speed of 1/5000 s and 1024 × 1024 pixel resolution (Photron PFV v.3.20). Because the insects were induced to fly in a trajectory parallel to the plane of the imaging sensor (deviations < 10°), we analyzed the trajectories in two dimensions (Additional file 1: Figs. S3, S4A). Pixels were scaled to metric coordinates using a 50 cm bar held horizontally at the same level of the flight trajectories. We recorded a clay ball in free fall to calibrate the vertical direction. Average vertical ball acceleration was 9.816 ± 0.14 m s^{−2}, less than 1% different from gravitational acceleration (9.805 m s^{−2}) and the camera was oriented to gravity so that the vertical ball drop direction of acceleration was always < 5% of 90°. Each insect was dropped several times (maximum five times) until we obtained two straight trajectories per animal. A resting period of approximately 20 min was left in between each flight to allow recovery. Males were frozen-killed at − 80 °C just after the experiments. After estimating the relative position of their center of mass (see below), males were pinned with their wings fully extended for morphometrics.

Video digitization was done by tracking morphological landmarks using the open source video analysis tool *DLTdv5* by T. Hedrick implemented in MATLAB (R2016b, MathWorks, Natick, MA, USA) [68] and the obtained data was analyzed using R (v 3.6.1) [69]. In *DLTdv5*, we used autotracking mode (predictor tool: extended Kalman) and manual tracking when the autotracking mode was unreliable. We marked the position of the head and of the terminal abdominal segment on each frame. Body pitch (°) was calculated for every frame by calculating the angle between the horizontal and the line linking the position of the head and that of the terminal segment (Additional file 1: Fig. S4B). Typically, in the first phase of the fall (free fall), body pitch decreases (i.e., the insect rotates forward, eventually diving head first), before the insect opens its wings (t_{1}) and actively corrects (t_{2}) and stabilizes its body pitch (phase 3) (Additional file 1: Fig. S4B). Body pitch was smoothed using a Savitzky–Golay filter with a polynomial order of 3 and a window size of 71 (*sgolayfilt*: ‘signal’, *function*: ‘R package’). The beginning of phase 2 (t_{1}) was determined as the time corresponding to the minimal body pitch (Additional file 1: Fig. S4B). The end of phase 2 (t_{2}) corresponded to the time when the insect’s pitch stabilized—i.e., when the rotational velocity (° s^{−1}) of body pitch reached a local minimum after a large peak corresponding to phase 2 (Additional file 1: Fig. S4B). We calculated the average rotational velocity during phase 2 (ω) as:

$$\omega =\frac{Body\, pitch \left({t}_{2}\right)-Body\, pitch ({t}_{1})}{{t}_{2}-{t}_{1}}$$

(1)

We used ω to quantify torsional agility—an important aspect of maneuverability—as it reflects how fast the animal can rotate to correct its body pitch in the air from a free falling, head first, position to a stable flight body pitch. This correction occurred over several wingbeats and was therefore distinct from within and among-wingbeat oscillations which could have indicated a lack of longitudinal stability [41]. The 2D position of the body center of mass was estimated using images taken in lateral view of the males (freshly dead) orthogonally balancing on a horizontal razorblade. Its position relative to the two landmarks (i.e., head and terminal segment) was then calculated which enabled us to define the center of mass of the individual on each flight video. Trajectories of the center of mass were analyzed using the package “trajr” in R [70]. Raw trajectories were smoothed using a Savitzky–Golay filter with a polynomial order of 3 and a window size of 31 (*TrajSmoothSG*: ‘trajr’; Additional file 1: Fig. S4A). Horizontal, vertical and composite velocities and accelerations were then computed on the smoothed trajectories (*TrajDerivatives*: ‘trajr’; Additional file 1: Fig. S4C, E). For each trial, we defined transient and steady states for both vertical and horizontal velocities (Additional file 1: Fig. S4D, E). In both cases, the transient state corresponded to the free fall and maneuver of the insect in the air during which body velocity greatly varied. The steady state started when velocity stabilized and acceleration started oscillating around 0 m s^{−2} (Additional file 1: Fig. S4D, E). We extracted the mean vertical and horizontal velocity (m s^{−1}) during their respective steady states. These measures were used to quantify the capacity of the insect to fly forward and ascend. On the videos, the position of the tip of the wing closest to the camera was also manually marked on frames corresponding to the end of upstroke and downstroke as this was sufficient to measure wing beat frequency (Hz) and stroke amplitude (°). Average wing beat frequency was calculated after the animal reached a stable body pitch (i.e., after t_{2}; Additional file 1: Fig. S4B). For each animal, we measured wing length and the position of the attachment of the wings relative to the head and tip of the abdomen using photographs and ImageJ. We then determined the position of the wing attachment point on each frame on the videos using this relative position between our two landmarks. The amplitude of wing strokes was calculated using trigonometry (Additional file 1: Fig. S7).

### Adhesive pads and substrate attachment performance

We used scanning electron microscopy (SEM) to observe the tarsi of the metathoracic leg of adult males and females, measure attachment pad areas, and describe the microstructures on these pads. Tarsi of the right metathoracic leg were cut off from 20 adult males and 20 adult females and fixed in 2.5% glutaraldehyde in PBS buffer for 24 h on ice on a shaker, dried in an ascending alcohol series, critical-point dried and sputter-coated with a 10 nm layer of gold–palladium. To obtain overview images, we used a rotatable specimen holder [71] and the scanning electron microscope (SEM) Hitachi TM3000 (Hitachi High-technologies Corp., Tokyo, Japan). The micrographs for visualization and measurements were taken at an acceleration voltage of 15 kV. The attachment microstructures on the tarsi of both sexes were further examined using the SEM Hitachi S4800 (Hitachi High-Technologies Corp., Tokyo, Japan) at 7 kV of acceleration voltage. Processing of the raw micrographs and measurements of projected attachment pad area (mm^{2})—i.e., the surface area of the tarsus specialized for adhesion and friction [72, 73]—were done using Photoshop CS6 (Adobe Systems Inc., San José, CA, USA).

To measure attachment forces (mN) in both pull-off (adhesion) and traction (friction) directions, we used 20 adult males (Mean ± S.D. = 0.46 ± 0.02 g) and 20 adult females (5.22 ± 0.31 g). A horsehair was glued to the metanotum of each insect and then attached to a 100 g force transducer (FORT100, World Precision Instruments, Sarasota, USA, linearity error: < 0.1%, resolution: 0.01%), connected to a BIOPAC model MP100 and TCI-102 system (BIOPAC Systems, Inc., Goleta, CA, USA), and mounted on a motorized micromanipulator (DC 3001R, World Precision Instruments Inc.). Maximum adhesion forces were recorded by vertically pulling the insects off a horizontal glass plate until they detached from the glass plate [45, 74]. The micromanipulator was moved upwards with a speed of 200 µm/s at a step size of 10 µm until the specimen was detached from the surface as indicated by an instantaneous drop in force. Maximum friction forces were recorded by horizontally pulling the insects backwards with the same retraction velocities as above, until detaching them from the glass plate [45, 75]. A glass plate was used as the substrate for the attachment force measurements, to eliminate mechanical interlocking of the claws with surface irregularities of rough substrates, or penetration of soft substrates. As glass is smooth on the microscopical level, this substrate enables estimation of the traction and pull-off performance of the attachment pads themselves on a standardized level, without the influence of substrate irregularities. Force–time curves were obtained using Acqknowledge 3.7.0 (BIOPAC Systems Inc., Goleta, CA, USA) and the maximum peaks were extracted as maximum adhesion forces, or maximum friction forces respectively. Each of the 20 males and females were measured three times in both directions on a glass plate and the average of the three measurements was used as the individual maximum adhesion/friction force. The order of the individuals was randomized and the substrate was cleaned between every measurement. The experiments were conducted at 20–23 °C and 50–60% relative humidity.

Static safety factors (i.e., SF_{static} = \(\frac{Attachment \,force}{Body\, weight}\)) were computed for each individual. Following the methods of Higham et al., 2017 [47], we estimated impact forces (F_{i}, in N) during landing using a model to eventually compute dynamic safety factors (i.e., SF_{dynamic} = \(\frac{Friction\, force}{{F}_{i}}\)). In our landing model, we assumed that males would stop immediately after landing on the tip of a leaf without slippage. F_{i} was calculated using the work-energy principle:

$${F}_{i}=mg+ \frac{m{v}^{2}}{2\delta }$$

(2)

where m = mass of the insect (kg), g = acceleration of gravity (9.81 m s^{−2}), v = landing speed of the insect (m s^{−1}) and δ = deflection of the leaf (m) (Fig. 4).

We used a simplified model of a leaf to estimate a range of values for the deflection of a leaf upon landing of a male leaf insect and explore potential values of SF_{dynamic}. The leaf and petiole were considered as a uniform cantilever beam with impact forces applying at the tip at length L. We considered the leaf to be initially horizontal. Deflection was calculated as a function of the weight of the insect. Given the relatively light weight of male leaf insects, and in contrast with [47], who were considering geckos (i.e., roughly 28 times heavier), we only accounted for small leaf deflections (deflection angle θ < 5°) and therefore used the following equation to estimate the deflection of the tip of the leaf δ [76]:

$$\delta =\frac{mg\,{L}^{3}}{3EI}$$

(3)

where L is the total length of the leaf (m) and EI is its flexural stiffness (N m^{2}). We used the range of values estimated by [47], for leaf length and corresponding EI (i.e., L_{1} = 23 cm, EI_{1} = 2.67 \(\times\) 10^{–3} N m^{2} and L_{2} = 86 cm, EI_{2} = 28.48 \(\times\) 10^{–3} N m^{2}) to account for the diversity of leaf mechanical properties and thus explore the possible range of SF_{dynamic}.

Finally, landing speed was estimated as a function of the body mass of the individual. Using our experimental flight videos (see above), we extracted the instantaneous speed of the experimental males right before making contact with the wood slab (i.e., the landing target) and built a linear regression between landing speed and body mass, including individual ID as a random factor (*lmer*: ‘lme4’, [77]). We found a significant positive effect of body mass on landing speed (see “Results”, Additional file 1: Fig. S6) and used the fixed-effect estimates (intercept and slope) from this model to predict the landing speed of the males for which we empirically measured attachment forces, given their body mass. This estimated landing speed was used in Eq. 2.

### Computational fluid dynamic simulations

To investigate how lift and drag produced by the insect were affected by size and shape, we generated 3D models of males of varying size and shape and predicted these forces during a steady and horizontal flight using computational fluid dynamics (CFD). We created a reference 3D surface model of an adult male body using photogrammetry. We pinned the body of a dead specimen in a flight posture (i.e., legs extended, forewings opened perpendicularly, antennae oriented 50° up, hindwings removed). We did not model flapping aerodynamics as fully integrating the complex three dimensional trajectories and aeroelastic deformations of the wings into a CFD model (e.g., [78]) was beyond the scope of this study. The aerodynamic interactions of the flapping wings with the body of insects in slow flight has been shown to be negligible (~ 5%) in [79]. Once dried, the individual was vertically mounted onto a pin on a custom-made turntable. 2D images using a DSLR camera (EOS 600D, Canon Inc., Tokyo, Japan) equipped with a macro lens (Canon EF 100 mm f/2.8 Macro USM), were then obtained from 100 different orientations (Additional file 1: Fig. S8). The 3D model was then reconstructed from these multiple images using Autodesk ReCap Pro 2019 (v5.0.4.17, Autodesk Inc., San Rafael, CA, USA) and subsequently smoothed and rewrapped using Autodesk Meshmixer 2017 (v11.5.474, Solid accuracy: 402, cell size: 0.202, density: 219, offset: 0.25, min thickness: 0.14 mm). From this reference model, we built four additional and artificial models using the “Move” tool in Meshmixer to either manually extend or shrink the abdominal lobes and therefore manipulate abdominal shape. These artificial shapes purposely spanned a wider range of body aspect ratios than the one found in actual males *P. philippinicum* (male natural range: 4.89–6.42, female natural range: 2.28–2.84, model range: 2.28–9.47). The model with the lowest aspect ratio displayed a female-like abdominal shape while the model with the highest aspect ratio had no abdominal expansions. The models were further scaled to a body length of 53 mm (mean male BL in our Montana population = 52.6 ± 0.25 mm) using Autodesk fusion 360 (v2.0.8335). From each of these five meshes, we created four additional models (25 models in total) respectively scaled to a factor 0.85, 0.95, 1.05 and 1.15. The insect models were tilted at a 44° body pitch in our control volume. This angle was determined using our flight experiments and a LMM with stable body pitch as the response variable, horizontal and vertical body velocity as main effects and individual ID as a random factor. Using the parameters estimated by this model, we predicted a body pitch of 43.9° for a vertical velocity of zero and a mean horizontal velocity of 1.57 m s^{−1} (i.e., the average horizontal velocity calculated from our flight trials, after the animal had stabilized its body pitch: 157 ± 6.8 cm s^{−1}).

We constructed a control volume around these body meshes in Autodesk CFD 2019 (v19.2), that provided numerical solutions to the Reynolds-averaged Navier–Stokes equations [50, 80]. A fluid volume was built around the mesh with walls far enough from the model mesh to avoid any reflection effects (1 × 0.5 × 0.5 m) (Additional file 1: Fig. S9A). The fluid was assigned the default properties of air in CFD 2019 (density at sea level = 1.205 kg m^{−3}, viscosity at 20 °C = 18.2 μPa s^{−1}). The phasmid models were then assigned properties of hardwood which, for mass-less and stationary models, should have no impact on results. The input flow on the anterior end of the control volume was set to 1.57 m s^{−1}. We held the air velocity around in the insect constant as we were only considering horizontal flight and average horizontal velocity after reaching a steady state did not significantly increase with size in our flight trials (Table 1, Additional file 1: Fig. S5A). We applied a zero-pressure condition on the opposing end of the volume. A slip/symmetry condition was applied to all other fluid boundaries. We automatically meshed the domain around the phasmid model, applied a surface refinement, and locally defined a non-uniform mesh refinement region (0.7 × 0.2 × 0.2 m) with a mesh size reduced to 75%, around and behind the model to better capture the resulting wake. We ran steady-state simulations using the turbulence model k-epsilon. The maximum number of iterations was set to 3000 although the simulations were stopped when they reached convergence according to the default convergence detection parameters of the CFD software (mean = 849 ± 50 interactions). The adaptive meshing tool was used to insure mesh optimization for our models and mesh independence of the results. The simulation was first run with the meshing parameters described previously. Then, the solution results of this simulation were automatically used to refine the mesh in high velocity gradient regions and rerun the simulation. We enabled the ‘flow angularity’ option to improve mesh resolution in areas with a lot of flow separation, the ‘free shear layers’ and ‘external flow’ options to refine the mesh in areas of strong velocity gradients. We ran three such cycles for each model. Final mesh sizes averaged 865,446 ± 74,771 nodes and 4,411,596 ± 379,780 elements (Additional file 1: Fig. S9). Finally, to help evaluate the validity of our simulations, we placed a sphere with the same Reynold’s number as the one calculated for the original phasmid model (Re = 5558) in a similar control volume with the exact same settings as our insect simulations. We found a drag coefficient (C_{d}) of 0.643. This is very close to the value predicted from experimental data (C_{d} = 0.652) and which was determined using the Eq. 8.83 in Morrison, 2013 [81].

The weight of the models (mN) with a non-modified abdominal shape (reference models) was determined from their BL using the linear regression built between male body mass and BL in our American population. To estimate the weight of the models with artificial abdominal shapes, we first measured their abdominal area relative to that of the reference model of identical BL and measured the average weight of the leaf-like abdominal expansions by cutting these extensions from five freshly frozen-killed males and measuring their areas and mass (89.7 ± 0.6 g m^{−2}). The CFD simulations estimated the aerodynamic forces (drag and lift, mN) that apply to the rigid insect body flying horizontally and steadily.

For each model, we measured the projected frontal area on a plane perpendicular to the air flow. We then calculated their coefficient of drag (\({C}_{D}= \frac{2 {F}_{drag}}{\rho {v}^{2}A}\)) and lift (\({C}_{L}= \frac{2 {F}_{lift}}{\rho {v}^{2}A}\)) using the model frontal area (A), the mass density of air (= 1.20473 kg m^{−3}), the velocity of the insect (v = 1.57 m s^{−1}) and the drag or lift force estimated from the CFD simulations. Lift to drag ratios (\(\frac{{C}_{L}}{{C}_{D}}\)) were also calculated for each model.

### Scaling of muscle power available and power required for flight

From theory, we estimated the scaling relationships of the power available (P_{a}) and the power required for flight (P_{r}) with body size using our empirical data. As leaf insects are slow flyers (average flight speed = 1.70 ± 0.07 m s^{−1}, mean ± SE), P_{r} mostly corresponds to the induced power P_{ind}—i.e., the cost for producing lift [20, 51]. P_{ind} is the product of the net required force from the wings to maintain the animal in the air and of the induced velocity in the wake. Following [82], we assumed that, under isometry, the induced velocity in the wake and weight-specific power required for slow flight should be proportional to the square root of wing disc loading (DL, N m^{−2})—i.e., body weight (W_{B}) divided by wing disc area (A_{WD}). Wing disc area—i.e., the area swept out by the wing during a wing beat cycle and through which air is accelerated downward to develop lift force—is determined by wing length (L_{w}) and stroke amplitude (\(\theta\), °) (Eq. 4).

$$DL=\frac{{W}_{B}}{{A}_{WD}}= \frac{{W}_{B}}{2\pi {{L}_{w}}^{2} \frac{\theta }{360}}$$

(4)

Thus, among geometrically similar animals, the power required for slow flight should scale as BL^{7/2} (Eq. 5).

$${\mathrm{P}}_{\mathrm{r}}\propto {W}_{B} ({DL)}^{1/2}\propto {W}_{B} \left({\frac{{W}_{B}}{2\pi {{{L}_{w}}^{2}} \frac{\theta }{360}}}\right)^{1/2}\propto {BL}^{3}\left({\frac{{BL}^{3}}{B{L}^{2}}}\right)^{1/2}\propto {BL}^{7/2}$$

(5)

In slow flight, P_{a} is the product of the net wing force and of the tangential velocity of the tip of the wing—i.e., angular velocity (rad s^{−1}) × wing length (m)—or the product of the muscle work—i.e., force \(\times\) distance of contraction—and of flapping frequency (Hz). We accepted the assumption that, for geometrically and dynamically similar organisms, force is proportional to the cross-sectional area of the muscles, which scales as BL^{2}, and distance of contraction scales as BL^{1}. Thus, work scales as BL^{3} and is therefore directly proportional to W_{B} [51, 83]. Flapping frequency is predicted to scale as BL^{−1} when the animal is using maximal or near-maximal effort [84]. Therefore, under isometry, P_{a} is expected to scale as BL^{2}.

To estimate the scaling exponent of the tangential velocity of the tip of the wing (i.e., angular velocity (rad s^{−1}) × wing length (m)) and of wing disc loading (Eq. 4) with L in leaf insects, we built LMMs with log_{10} tangential velocity or log_{10} disc loading as the response variable, log_{10} BL as the fixed effect and male ID as a random factor. Tangential velocity of the tip of the wing did not significantly scale with BL or differ from isometrical expectations (β = − 1) in our flight experiments (β = − 2.43 ± 1.60, \({\chi }^{2}\) = 2.46, df = 1, p = 0.12). Similarly, the observed scaling exponents of body and flight muscle mass did not significantly differ from isometric expectations (Additional file 1: Table S1). Coupling this with the aforementioned assumptions about muscle force, we estimated P_{a} scaled isometrically in leaf insects (β = 2). Wing disc loading (DL) positively scaled with BL (\({\chi }^{2}\) = 8.99, df = 1, p = 0.003). While, under isometry, we expected a scaling exponent of 1 (Eq. 4), we found that DL increased disproportionately with BL (β = 4.96, 95% CI [1.99, 7.94]). This is a consequence of the reduced wing stroke amplitude seen in larger individuals (Table 2). As the induced velocity (V_{i}) in the wake scales proportionately with the square root of wing disc loading, we estimated that P_{r} scaled more steeply with size than expected under isometry (β = 5.48 vs 3.5) (Eq. 6)

$${\mathrm{P}}_{\mathrm{r}}\propto {W}_{B} ({DL)}^{1/2}\propto {BL}^{3} {({BL}^{4.96})}^{1/2}\propto {BL}^{5.48}$$

(6)

Consequently, ∆P (P_{a} − P_{r}) decreased with body size more rapidly than would be expected under isometry.

### Statistical analyses

All statistical analyses were run in R version 3.6.1 [69] and all statistical tests were two-sided. For all linear models, we systematically checked the normal distribution of the residuals and the absence of any specific patterns in their distribution.

We tested for sex differences in mean BL, body mass, body aspect ratio and antenna length using Wilcoxon–Mann–Whitney tests (*wilcox.test*: ‘stats’). To test for sex differences in the scaling relationships (i.e., in slope and intercept) of the various morphological traits, attachment and aerodynamic force measurements with BL, we built ordinary least square regressions [85] including each log_{10}-transformed trait as response variable and log_{10} BL, sex and their interaction as predictor variables (*lm*: ‘stats’). Type-I ANCOVAs were used to determine significance of the fixed effects (*anova*: ‘stats’). Departure from isometry which corresponds, on a log–log scale, to a slope of 1 for linear measurements, 2 for areas and 3 for masses [86], was tested using 95% confidence intervals (CI) around the estimated regression slopes (*confint*: ‘stats’). We similarly built linear models to investigate the effect of body mass, leaf size (i.e., small or large) and their interaction on the estimated landing impact forces in males and dynamic safety factors.

To test for effects of body size, wing size and body shape on male flight performance, we built linear mixed models (LMM) (*lmer*: ‘lme4’). Response variables were either rotational velocity (ω), mean vertical or horizontal velocity, wingbeat frequency, or wing stroke amplitude. Body mass, wing area and body aspect ratio were mean-centered and standardized (\(\mu\) = 0 and σ = 1, *scale*: ‘base’) and were included as main fixed effects. Individual ID was included as a random factor to account for replications of each individual. Likelihood ratio tests were subsequently performed sequentially to assess the significance of the fixed effects (*anova*: ‘lme4’). For response variables significantly affected by body mass and wing area and to illustrate their combined effect, we built and plotted similar LMMs but with wing loading as the only fixed effect.

Following our CFD simulations, we tested for the effects of body size and shape on lift to drag ratio, C_{D}, C_{L}, and relative aerodynamic forces applying on the body using linear models including BL and body aspect ratio as explanatory variables. Sequential ANOVAs (type I) were subsequently performed to assess significance (*anova*: ‘stats’). In models including relative force estimates as response variables, variables were log_{10}-transformed to compute scaling exponents [86].