Interspecific variation in bristle number on forewings of tiny insects does not influence clap-and-fling aerodynamics

The wings of flying insects show tremendous diversity in shape, size and function. Curiously, the wings of several families of flight-capable insects smaller than fruit flies have independently evolved ptiloptery (Polilov, 2015; Sane, 2016), meaning wings with long setae at the fringes. Though their extremely small sizes (body length <2 mm) make visual observation difficult, tiny flying insects are not limited to just a few outlying examples. Rather, more than 5500 species of thrips (Thysanoptera; Morse and Hoddle, 2006), as well as several hundred species of bristle-winged wasps (Trichogrammatidae, Mymaridae, Mymarommatidae; Heraty et al., 2013), have been identified to date. Despite their importance as biological vectors of plant viruses and as invasive pests of commercially important plants (Ullman et al., 2002; Jones, 2005), we still understand little of the flight mechanics of tiny insects. Owing to the difficulty in acquiring free-flight recordings of tiny insects, several studies have used physical and computational modeling to examine the functional significance of wing bristles (Santhanakrishnan et al., 2014; Jones et al., 2016; Lee and Kim, 2017; Kasoju et al., 2018). While these studies have shown that having bristles aids flight at such small sizes, little is known about the extent of variation in bristled wing morphology among different species of tiny insects. Moreover, it remains unclear whether tiny insects experience selective pressure to optimize the mechanical design of their bristled wings, particularly given the extreme challenges of flight at miniature body sizes.

Pronounced viscous dissipation of kinetic energy occurs at wing length scales on the order of 1 mm, making it difficult for tiny insects to stay aloft. The relative importance of inertial to viscous forces in a fluid flow is characterized using the dimensionless Reynolds number (Re=ρVL/μ), where ρ and μ are the density and dynamic viscosity of the fluid medium, respectively, and V and L are characteristic velocity and length scales, respectively. The length scale has been examined based on wing chord (i.e. L=c; Re_c) and bristle diameter (L=D; Re_b), with Re_c on the order of 1 to 10 and Re_b ranging between 0.01 and 0.07 (Ellington, 1975; Kuethe, 1975; Santhanakrishnan et al., 2014; Jones et al., 2016). Despite the difficulty of sustaining flight at such low Re, entomological studies have reported active flight and dispersal of thrips (Morse and Hoddle, 2006; Rodriguez-Saona et al., 2010). Tiny insects use biomechanical adaptations to overcome the fluid dynamic challenges associated with flight at small scales. These insects operate their wings at near-maximum stroke amplitude using the ‘clap-and-fling’ mechanism, first observed by Weis-Fogh (1973) in Encarsia formosa (Hymenoptera). The use of clap-and-fling has been documented in other freely flying tiny insects, including Thrips physapus (Thysanoptera; Ellington, 1975) and Muscidifurax raptor (Hymenoptera; Miller and Peskin, 2009). Wing rotation during fling has been noted to augment lift via the generation of a leading edge vortex on the wings (Weis-Fogh, 1973, 1975; Lighthill, 1973; Spedding and Maxworthy, 1986; Dickinson et al., 1999; Birch et al., 2004; Miller and Peskin, 2005, 2009; Lehmann et al., 2005; Lehmann and Pick, 2007; Arora et al., 2014; Santhanakrishnan et al., 2018). However, the concomitant generation of large drag force at the start of fling undermines the advantage of clap-and-fling at Re_c relevant to tiny insect flight (Miller and Peskin, 2005; Arora et al., 2014). Previous studies have thus examined the flow structures and aerodynamic forces generated by bristled wings in comparison with solid wings (Sunada et al., 2002; Santhanakrishnan et al., 2014; Jones et al., 2016; Lee and Kim, 2017; Lee et al., 2018; Kasoju et al., 2018; Ford et al., 2019), showing that bristled areas on the wings can reduce the force required to fling the wings apart.

Despite this focus on modeling, morphological variation of bristled wing design in tiny flying insects is far less documented. Jones et al. (2016) examined the inter-bristle gap (G), bristle diameter (D) and wing area covered by bristles in the forewings of 23 species of fairyflies (Hymenoptera: Mymaridae and Mymarommatidae). With decreasing body length (BL), they found that G and D decreased and area occupied by bristles increased. Moreover, Ford et al. (2019) found that the percentage of solid membrane area (A_M) to total wing area (A_T) in the forewings of 25 species of thrips ranged from 14% to 27%, as compared with the percentage of A_M to A_T ranging from 11% to 88% in smaller-sized fairyflies examined by Jones et al. (2016). Yet interspecific variation of G, D, wingspan (S) and number of bristles (n), as well as their concomitant effects on clap-and-fling aerodynamics, are currently unknown.

Such variation in wing morphology across species may arise from many factors. Adaptation drives much interspecific variation (Futuyma and Kirkpatrick, 2017), and many studies have thus focused on the consequences of variation for optimal functional performance. For example, Ford et al. (2019) used physical models to test the aerodynamic consequences of variation in proportion of solid (i.e. compared with bristled) area on wings. They showed that lift-to-drag ratios were largest for bristled wing models with proportions similar to thrips forewings, suggesting that selection may maintain the small range of variation in thrips. Alternatively, variation among species may have little adaptive explanation (Gould and Lewontin, 1979). Contingent factors in evolution may cause distantly related groups to differ, even under the same selective pressures (Gould, 2002; Blount et al., 2018). Thus, high phylogenetic inertia may explain why species from differing clades differ in phenotype (Hansen and Orzack, 2005). Paradoxically, shared evolutionary history can also explain variation among more closely related species. Such species often share factors (e.g. developmental, genetic) that have similar effects on different traits; when one such trait varies among species, the other will likewise vary. For example, shared growth factors underlying different body parts can cause them to covary with body size. If closely related species differ in selection for body size, then they will similarly differ in traits that grow with body size during development. Strong scaling relationships (i.e. allometry) may indicate evolutionary history as a source of interspecific variation (Pélabon et al., 2014). Thus, accounting for phylogenetic relationships and estimating evolutionary inertia can also help explain variation among species.

In this study, we quantified variation in morphology across species of bristle-winged insects and addressed the factors potentially driving this variation. We first measured wing morphology from 59 species of thrips and fairyflies. We then conducted phylogenetic regressions of key variables on body length and we quantified evolutionary inertia. Using the morphological data as a guide for biologically relevant variation, we then fabricated physical bristled wing models varying in G, D, maximum wingspan (S_max) and n. These physical models were comparatively tested using a dynamically scaled robotic platform mimicking the portion of clap-and-fling kinematics where wing–wing interaction occurs. Aerodynamic force measurements and flow field visualization were conducted to identify the functional significance of the above bristled wing design variables. Because of the high variation in n and G/D despite the extreme aerodynamic demands of flight at small size, we hypothesized that at Re relevant to tiny insect flight, dimensionless aerodynamic forces generated by clap-and-fling would be minimally impacted by variation in n and G/D within their biological ranges. If true, tiny flying insects may not experience selective pressure to further functionally optimize the mechanical design of their bristled wings.

We measured average BL, A_T, S_max, n, G and D from published forewing images of thrips and fairyflies, whose size ranged from 0.1 to 2 mm in BL. In the Supplementary Materials and Methods, we detail our criteria for choosing published forewing images for measurement. Based on these criteria, we selected forewing images of 16 thrips species for measuring S_max, A_T and n, and of 22 different thrips species for measuring G and D (Mound and Reynaud, 2005; Mound, 2009; Zhang et al., 2010; Riley et al., 2011; Cavalleri and Mound, 2012, 2014; Ng and Mound, 2015; Masumoto et al., 2013; Minaei and Aleosfoor, 2013; Zamar et al., 2013; Dang et al., 2014; Cavalleri et al., 2016; Lima and Mound, 2016a,b; Mound and Tree, 2016; Wang and Tong, 2016; Goldaracena and Hance, 2017; PaDIL: http://www.padil.gov.au). The thrips species considered here encompassed three different taxonomic families. In addition, we selected 21 fairyfly species for measuring S_max, A_T and n (Huber et al., 2006, 2008; Huber and Baquero, 2007; Lin et al., 2007; Huber and Noyes, 2013), largely overlapping those of Jones et al. (2016), who presented data on G and D for 23 species.

We measured bristled wing morphological variables from these images using ImageJ software (Schneider et al., 2012). S_max was defined to be the distance from the center of the wing root to the tip of the bristles, following Fig. 1A. Average wing chord (c_ave) was calculated by measuring A_T using the same procedure as in Jones et al. (2016) and Ford et al. (2019), then dividing A_T by S_max. G/D ratio was calculated from the measurements of G and D in the forewing images. BL measurements were made either from images (where available) or from the text of the article containing the image. A full list of species, corresponding measurements and publication sources of the original images is provided as Appendix S1 in Figshare (https://doi.org/10.6084/m9.figshare.14478108.v1).

We accounted for shared evolutionary history among species in our regressions by using phylogenetic generalized least squares (PGLS; Martins and Hansen, 1997). Regressions were fit with the maximum-likelihood value of λ (Pagel, 1999), the phylogenetic signal of regression residuals. This procedure best balances species similarity owing to shared history and shared adaptation (Hansen and Orzack, 2005), which improves statistical inference (Revell, 2010). Moreover, λ can be used as a metric of the role of evolutionary history in a fitted relationship (Hansen and Orzack, 2005).

Phylogenetic data for our study species were scarce. Only nine of our 59 species of thrips and fairyflies were included in published phylogenies, and these nine were scattered across published trees (Munro et al., 2011; Buckman et al., 2013; Lima and Mound, 2016a,b; Pereyra et al., 2019). Thus, we simulated many possible phylogenies for our study species and conducted comparative analyses across these trees. This procedure allowed for both integration over phylogenetic uncertainty (Martins, 1996) and assessment of the sensitivity of our results to any specific potential phylogeny (Losos, 1994). Herein, we briefly summarize our procedure for simulating phylogenies. We refer readers to the Supplementary Materials and Methods for detailed simulation methods, justification and discussion of why phylogenetic regressions should be robust to variation or error in phylogeny.

We constrained our simulated trees to fit current taxonomic knowledge, as adding some phylogenetic structure increases accuracy over completely random approaches (Housworth and Martins, 2001; Martins, 1996; Martins and Housworth, 2002; Symonds, 2002). This meant, for example, that all species of a given genus were each other's closest relatives in every simulated tree. For thrips, taxonomic information was extracted from the comprehensive Thrips Wiki (https://thrips.info/wiki/; accessed 15 March 2021). Fairyflies are likely a polyphyletic group of two families in two superfamilies of wasps (Mymarommatoidea: Mymarommatidae and Chalcidoidea: Myrmaridae; Huber, 1986; Davis et al., 2010; Munro et al., 2011); in simulations we assumed these two families to be each other's sister taxon. Genera for these two families were extracted from taxonomic accounts (Gibson et al., 2007; Huber, 2005, 2017; Lin et al., 2007; Poinar and Huber, 2011). Phylogenies were simulated in the package phytools v.0.7-70 (Revell, 2012) in R v.4.0.2 (https://www.r-project.org/). We simulated 10,000 trees, then pruned each tree to only include the species for which we had phenotypic data, which varied based on the response variable. All tree simulation R code, taxonomic information and resulting trees are included in Figshare as Appendices S2–S4 (https://doi.org/10.6084/m9.figshare.14478108.v1).

Regression analyses were conducted on logged variables, as is standard in body-size scaling analyses (Voje and Hansen, 2013; Pélabon et al., 2014; Glazier, 2021). For each simulated tree, we compared four nested models: (1) a null model with only an intercept; (2) a simple model of regression in which both thrips and fairyflies shared all parameters; (3) a model in which both groups shared a scaling slope but had different intercepts; and (4) a full model in which both groups differed in slope and intercept. These models thus allowed us to estimate scaling relationships between variables and ask whether such relationships differed in thrips and fairyflies (Gartner et al., 2010; Moen et al., 2016). All regressions were estimated in the package phylolm v.2.6.2 (Ho and Ané, 2014). We compared models for each tree with the corrected Akaike's information criterion (AICc) and its associated weights (Burnham and Anderson, 2002). We used the model weights to calculate model-averaged regression parameters, adjusted R² and λ values (Burnham and Anderson, 2002; Posada and Buckley, 2004). We then averaged these values across trees, as well as the AICc values and model weights. Assuming that each randomly resolved tree is equally likely, such means represent values integrated over phylogenetic uncertainty (Martins, 1996). We also calculated the 95% confidence intervals of slopes, accounting for both estimation and phylogenetic uncertainty (Martins, 1996). Finally, we calculated the proportion of trees in which a scaling model (i.e. models 2–4) had the highest weight. This proportion reflected the effect of phylogenetic structure on finding a non-zero scaling relationship (Losos, 1994).

Variation in wingspan, bristle number and G/D at different body lengths motivated our subsequent physical model experiments. However, we designed these models at a chord-based Re, rather than body length. Moreover, our experiments held two variables constant (e.g. wingspan and bristle number) while varying a third (e.g. G/D). Thus, we also examined PGLS correlations between these variables, likewise calculating means across the simulated phylogenies, as above. We estimated these correlations using custom R code from Moen et al. (2013), following Rohlf (2006). All R code for regression and correlation analyses, as well as for producing the resulting figures, is provided in Appendices S5–S6 on Figshare (https://doi.org/10.6084/m9.figshare.14478108.v1).

The bristled wings tested in this study were simplified to rectangular shape with constant wing chord (c in Fig. 2A) to minimize variability in confinement effects along the wingspan from the tank walls. The percentage of A_M/A_T in all the models was maintained at 15%, which is in the range of A_M/A_T of thrips and fairyflies (Ford et al., 2019). Bristle length (L_b, see Fig. 2A) and w were maintained as constants on either side of the membrane for all 23 wing models tested. The values of constants c, L_b and w are provided in Table S1.

Scaled-up physical models were used in this study to examine the roles of bristled wing geometric variables on clap-and-fling aerodynamics at Re_c=10. We used this approach to overcome the difficulty of resolving the flow around and through a bristled wing on the scale of 1 mm length. As we did not match the values of dimensional geometric variables to those of real insects, we used geometric similarity to match non-dimensional variables (n, G/D) in all the physical models to be in the range of tiny insects. As n depends on G, D and S per Eqn 1, the choices of non-dimensional variables include n, G/D, G/S and D/S. We chose G/D to match Jones et al. (2016). In addition, to understand the isolated role of each dimensional variable, we tested scaled-up models varying in G, D and S. For each condition, we maintained the two other dimensional variables as constants and also matched the non-dimensional variables (n, G/D) to be within their biologically relevant ranges identified from morphological analysis. For details on the fabrication details of bristled wing models, refer to the Supplementary Materials and Methods.

The dynamically scaled robotic platform used in this study (Fig. 3A,B) has been described in previous studies (Kasoju et al., 2018; Ford et al., 2019) and experimentally validated against results in Sunada et al. (2002) corresponding to a single wing in translation at varying angles of attack (in Kasoju et al., 2018). For more details on the robotic platform and justification of our forewing approach, refer to the Supplementary Materials and Methods.

Free-flight recordings adequate for characterizing instantaneous wing kinematics are unavailable for most species of tiny insects. Thus, we used a modified version of 2D clap-and-fling kinematics developed by Miller and Peskin (2005). The simplified kinematics used here do not capture: (1) 3D flapping translation during the downstroke and upstroke, and (2) wing rotation at the end of the downstroke (‘supination’). In real insects, the flapping cycle includes the combination of wing revolution (which we referred to as ‘3D flapping translation’ following terminology in Sane, 2003), wing rotation and elevation with respect to the root of the wing. In our study, the wings rotated and translated along a horizontal line with no change in elevation or stroke angle (Fig. 3C,D). ‘Wing rotation at the end of downstroke’ refers to the ventral stroke reversal (supination) at the end of downstroke that is observed in 3D flapping flight. In this study, a ‘stroke cycle’ is defined as a clap stroke and fling stroke (the latter corresponding to pronation or dorsal stroke reversal) and does not include the ventral stroke reversal occurring towards the end of downstroke. Similar or modified forms of these kinematics have been used in several other studies (Miller and Peskin, 2004, 2009; Santhanakrishnan et al., 2014; Arora et al., 2014; Jones et al., 2016; Kasoju et al., 2018; Ford et al., 2019; Kasoju and Santhanakrishnan, 2021). Fig. 2B shows the motion profiles prescribed for a single wing, where dimensionless velocity (instantaneous wing tip velocity U divided by steady translational velocity U_ST) is provided as a function of dimensionless time (τ) during rotational and translational motion. τ was defined as τ=t/T, where t represents instantaneous time and T represents time taken to complete one cycle of clap-and-fling. The motion profile for the other wing was identical in magnitude but opposite in sign, so that the wings would travel in opposite directions. Both wings moved along a straight line (no change in elevation and stroke angles). Schematic diagrams of the clap phase (Fig. 2C) and fling phase (Fig. 2D) are provided to show the direction of motion and wing position at the start and end of each portion of each half-stroke. The wings were programmed to start from an initial position corresponding to the start of the clap phase, and this was followed by the wings moving toward each other until the start of the fling phase, after which the wings moved apart from each other. The distance between the wings at the end of the clap phase was set to 10% of chord length, which we justify in the Supplementary Materials and Methods. In addition, the wingbeat kinematics are undescribed for most species of tiny insects and are likely variable across species (Lyu et al., 2019). For the present study, we prescribed 100% overlap between rotation and translation during both clap and fling, meaning that the wings translated during the entire rotational time. This was because previous studies (Arora et al., 2014; Kasoju and Santhanakrishnan, 2021) have shown that high overlap between rotational and translational motions significantly increases the aerodynamic forces (both lift and drag).

Two-dimensional time-resolved particle image velocimetry (2D TR-PIV) measurements were conducted to characterize the flow generated during clap-and-fling motion by bristled wing pairs along the chordwise plane (data acquired along a horizontal plane shown in Fig. 3A). Two-dimensional TR-PIV based two-component velocity vector fields were also used to determine the strength (i.e. circulation) of the leading edge vortex (LEV) and the trailing edge vortex (TEV). Two-dimensional phase-locked PIV (2D PL-PIV) measurements were conducted to characterize flow leaked along the span of bristled wings (data acquired along two vertical planes shown in Fig. 3C). For more details on validation of 2D flow simplification, the experimental arrangements and processing steps used for 2D TR-PIV and 2D PL-PIV measurements, refer to the Supplementary Materials and Methods.

In some cases, it may be possible to directly estimate the reverse (i.e. leaky) viscous volumetric flow rate in the direction opposite to bristled wing motion from the 2D PL-PIV data. However, we were not able to calculate this flow rate directly because high-magnification images would be needed to resolve flow through inter-bristle gaps (i.e. on the order of a few millimeters). This conflicted with our desire to use lower magnification in order to resolve flow across the entire wingspan (i.e. 10× greater than G) for calculating Q_viscous across a bristled wing.

Most variables showed considerable diversity across species. In thrips, S_max ranged from 305 to 1301 μm and bristle number (n) ranged from 44 to 161 (Fig. 1B,C). In fairyflies, S_max ranged from 180 to 1140 μm and n ranged from 32 to 104 (Fig. 1B,C). S_max increased with body length with negative allometry, meaning that larger individuals had relatively shorter wings than smaller individuals (Fig. 1B, Table 1). Most model weight across phylogenies indicated support for a model with the same slope and intercept for thrips and fairyflies (Table S2). n increased with body length similarly in both groups (Fig. 1C, Table 1), though there was nearly equivalent support for similar versus differing intercepts in the groups (Table S2). The latter meant more bristles at the same body length in fairyflies (Fig. 1C). In both S_max and n, however, we found that AICc model weight was concentrated on the two models with the same slopes for the two groups, which suggests similar scaling relationships. In contrast, while the inter-bristle gap to bristle diameter ratio (G/D) decreased with body length across both groups (Fig. 1D), the model with the most weight had a different slope and intercept for the two groups (Table S2). G/D more strongly decreased with increasing body length for the larger-sized thrips species (Fig. 1D, Table 1). The model in which both groups shared a slope and intercept also showed high statistical support across trees (Table S2). Regardless of the optimal model, these results mean that larger animals have more tightly packed bristles, with less leakage. Phylogenetic signal (λ) was close to 1 in S_max (i.e. residual species similarity reflects phylogeny), nearly 0 in n (i.e. similarity is independent of phylogeny) and intermediate in G/D.

Overall, our results suggest that both groups follow shared trends in bristle variables with BL across bristle-winged insects. Yet only BL strongly predicted S_max, with R²_adj almost two times lower for both n and G/D (Table 1). These latter results made us predict that variation in these latter two variables would have less aerodynamic consequences than S_max, motivating our robotic model experiments. Given weak correlations among S_max, n and G/D (Table S3, Fig. S1), we probed the effect of varying each of these variables while holding the other two constant.

For all the wing models tested, C_D and C_L were observed to follow the same trends over time during both clap and fling (Fig. 4A,B). Peak C_D occurred during fling (τ∼0.6) in all wing models (Fig. 4A). This time point corresponds to the end of rotational acceleration and translational acceleration (Fig. 2B), such that the wing pair would experience larger viscous resistance. C_D was found to drop after τ∼0.6 until the wing rotation ended (τ∼0.73) for all the wing models (Fig. 4A). Just before the C_D reached the negative value at the end of fling where the wings decelerate, we observed C_D to plateau from τ∼0.73 to 0.84 (Fig. 4A). This time corresponds to steady translational motion of the wings (Fig. 2B), where the wings translate with constant velocity at 45 deg angle of attack. Most of the drag during a cycle was generated in fling. Temporal variation of C_D was lower during clap half-stroke (τ=0–0.5) as compared with fling (Fig. 4A).

Three positive C_L spikes were observed in all wing models (Fig. 4B): (1) τ∼0.6 in fling, similar to that of peak C_D; (2) start of clap (τ∼0.16); and (3) end of clap (τ∼0.38). τ∼0.16 corresponds to the end of translational acceleration at a 45 deg angle of attack and τ∼0.38 corresponds to the end of rotational acceleration during clap (Fig. 2B). Peak C_L occurred during fling for all wing models. Unlike for drag force, both clap and fling half-strokes contributed almost equally to lift generation.

Both C_D and C_L decreased with increasing G and decreasing D (Fig. 4i,ii). Increasing S increased both C_D and C_L (Fig. 4iii). When increasing n for constant G/D, both C_D and C_L were found to increase (Fig. 4iv), particularly at the beginning of the fling phase. In contrast, increasing G/D for constant n decreased both C_D and C_L (Fig. 4v). Across all the wing models tested, we observed noticeable negative lift towards the end of fling. This is due to the wings not coming to complete rest and performing stroke reversal to position the wings for clap for the next cycle.

Cycle-averaged force coefficients (⁠⁠) were used to examine how each geometric variable impacted aerodynamic forces in a complete cycle (Figs 5 and 6). Individually increasing G and D showed little to no variation in when considering the standard deviations (Fig. 5A,B). decreased with increasing G and showed little to no variation with increasing D (Fig. 5A,B). Both and increased with increasing S from intermediate to large values of S (Fig. 5C). increased with increasing n (Fig. 6A). increased with n, most notably at n>88, though it plateaued between some consecutive values (Fig. 6A). Increasing G/D showed little to no variation in and when considering the standard deviations (Fig. 6B), though extreme values of G/D slightly differed.

Spanwise distribution of horizontal velocity (u) was examined near the instant of peak C_D (τ∼0.63) from 2D PL-PIV velocity fields (Fig. 7A). Looking at the extremes of each test condition, u increased with: (1) decreasing G; (2) increasing D; (3) increasing S; (4) increasing n; and (5) decreasing G/D. This reveals how each variable (i.e. G, D, S, n, G/D) differentially affects flow through a bristled wing. Similar to C_D, Le was observed to peak during fling. During the fling half-stroke, Le peaked either at τ∼0.56 or τ∼0.63 for all wing models (Fig. 7B) where the wings were near the end of rotational acceleration (Fig. 2B). Similarly, wing deceleration during fling from τ∼0.69 to τ∼0.88 resulted in a drop in Le (Fig. 7B). During steady wing translation from τ∼0.75 to τ∼0.82, Le was found to almost plateau in all wing models.

Le was larger in early clap (τ∼12.5) right after the wing pair started from rest, with minimal time for boundary layers around each bristle to be well developed. Thereafter, Le decreased with increasing clap duration until τ∼0.38, corresponding to the end of rotational acceleration (Fig. 2B). This latter observation in clap is in direct contrast to the peak in Le during fling, which was observed at the end of rotational acceleration. This disparity can be explained by examining the prescribed wing motion. In clap, wings were prescribed to translate first at a 45 deg angle of attack and then rotate. This provides ample time for the generation of shear layers around the bristles that block inter-bristle flow (see Kasoju et al., 2018 for a detailed discussion). Both rotation and translation started simultaneously in fling, necessitating more time for shear layers to develop around the bristles.

Maximum Le (Le_max) increased with increasing G and decreasing D (Fig. 7Bi,ii). However, changes in Le were comparatively small for the range of variation in G and D tested in this study. Similar to force coefficients (Fig. 4iii), increasing S did not show any particular trend for Le (Fig. 7Biii). However, if we look at the extreme wingspans (67.5 and 94.5 mm), Le was found to increase with increasing S. Increasing n for constant G/D was found to decrease Le. Changing G/D for constant n showed little to no Le variation.

Velocity vector fields overlaid on out-of-plane vorticity contours (ω_z) showed the formation of LEV and TEV over the wing pair during clap and fling half-strokes (Movies 1–3). Vorticity in the LEV and TEV increased near the end of clap and in early fling, when the wings were in close proximity of each other (Fig. 8B–D). This suggests that wing–wing interaction plays an important role in LEV and TEV formation, which in turn impacts force generation. Circulation (Γ) of both the LEV and TEV showed little to no variation with changing G, D and S. Peak Γ for both the LEV and TEV occurred in fling (τ=0.65), near the end of both translational and rotational deceleration (Fig. 2B). This was followed by a decrease in Γ of both LEV and TEV with increasing fling time (Fig. 8B–D). Γ of the LEV and TEV increased slowly in time during clap and reached a maximum near the end of the clap (τ=0.35), corresponding to the start of translational deceleration and end of rotational acceleration. The latter was identical to the instant at which peak Γ occurred in fling.

From the prescribed kinematics (Fig. 2B), peak rotational acceleration started early in fling, whereas it started later in clap. This could be why Γ peaked early in fling and later in clap. This suggests that wing rotation plays a dominant role in LEV and TEV development. Also, both wings are in close proximity during the later stages of clap and early stages of fling, suggesting the importance of wing–wing interaction in LEV and TEV development. Thus, wing rotation in concert with wing–wing interaction augments LEV and TEV circulation during both clap and fling half-strokes.

Recent studies have shown that bristled wings provide drag reduction in clap-and-fling at Re_c relevant to tiny insect flight (Santhanakrishnan et al., 2014; Jones et al., 2016; Kasoju et al., 2018; Ford et al., 2019). However, n, S_max and G/D had not been measured in different families of tiny insects, and their individual effects on aerodynamic forces were unclear. From our analysis of variation across thrips and fairyflies, we found that S_max and n increased with BL in both groups. We also found that G/D decreased with BL in both groups, but more strongly in thrips. Within the biologically relevant range of n and G/D, we found that: (1) increasing G provides more drag reduction as compared with decreasing D, (2) changing n for constant G/D has little variation on lift generation for n<100, and (3) changing G/D for constant n minimally impacts aerodynamic forces. The minimal influence of n and G/D on clap-and-fling aerodynamics, despite broad biological variation, suggests that tiny insects may experience lower biological pressure to functionally optimize n and G/D for a given wingspan.

Variation among related species can stem from many factors: evolutionary history, correlated response in selection to other traits, physical constraints associated with body design and function, and adaptation to variation in body size, ecology or environment (Gould and Lewontin, 1979; Alexander, 1985; Taylor and Thomas, 2014). In the case of bristled wing morphology of tiny insects, most studies have examined physical constraints and adaptation, that is, whether interspecific variation has consequences for flight aerodynamics, possibly driven by variation in body size. For example, Ford et al. (2019) reported a narrow range of A_M/A_T (14–27%) across 25 thrips species, but much higher variation across fairyflies. In both groups, A_M/A_T showed a strong, positive relationship with body length. At Re_c relevant to tiny insect flight, they found the highest aerodynamic efficiency (lift-to-drag ratio) for A_M/A_T in the range of thrips forewings and lower aerodynamic efficiency outside the range, perhaps facilitating flight in the larger-bodied thrips.

In this study, we found that both S_max and n increased with increasing BL in thrips and fairyflies (Fig. 1B,C). Interestingly, the ranges of S_max largely overlapped across fairyflies and thrips, despite differences in BL (most thrips BL>1 mm; all fairyfly BL<1 mm). This suggests that there could be a limit to the aerodynamic benefits of increasing wingspan. Moreover, we found that phylogenetic signal in the regression residuals (λ) was high for S_max on BL (Table 1), which explained the high R² value despite much scatter about the regression line (i.e. phylogeny explained much of the residual variation in Fig. 1B). In other words, closely related species were similar in the way they deviated from the regression line (Revell, 2010), which suggests that underlying growth factors in common with BL may be ultimately driving variation in wingspan across closely related species. If selection favors a change in body size, then wingspan may similarly change.

Values of n were concentrated in the range of 60–90 for the species of thrips and fairyflies that we examined, corresponding to a large BL range of 300–1700 μm. Moreover, the relationship between n and BL was relatively weak (R²_adj=0.350; Table 1). These observations led us to hypothesize that n may not need to be optimized to fall within a narrow range for a given BL toward improving aerodynamic performance. Consistent with this hypothesis, our robotic models showed insensitivity of aerodynamics to this range of n. The weak phylogenetic signal in regression residuals (Table 1) suggests little influence of evolutionary history (Hansen and Orzack, 2005). Therefore, the factors affecting the evolution of bristle number remain unclear.

Jones et al. (2016) previously showed no relationship between G/D and BL in fairyflies. However, our analyses suggest that there is an overall reduction in G/D with size in bristle-winged insects, with a steeper decline in thrips (Fig. 1D, Table 1). This difference in our results and those of Jones et al. (2016) stemmed from both our use of phylogenetic analyses and from including the larger thrips, which revealed an overall trend across taxa. That said, this pattern was still relatively weak (R²_adj=0.376; Table 1), with much variation in G/D at a given BL. Previous studies have reported that both lift and drag forces increase with decreasing G/D (Jones et al., 2016; Kasoju et al., 2018). This result could explain the more steeply negative relationship between G/D and BL in thrips, the larger of the two groups: as body mass increases, more lift is necessary to allow flight. Yet the high variation in G/D at long BL in fairyflies raises a question as to whether their G/D needs to be optimized for improving aerodynamic performance. In particular, we currently lack observations of fairyflies in free flight and thus do not know how or to what extent they use flapping flight. An intriguing possibility is that fairyflies facultatively parachute, and their wing structure better reflects the selective demands of that behavior. Thrips have been observed to facultatively parachute (Santhanakrishnan et al., 2014), increasing the probability that fairyflies do so as well.

Physical model studies of flapping flight match Re_c of the experiments to biological values to achieve dynamic similarity. Specific to the bristled wings of interest to this study, dynamic similarity of inter-bristle flow characteristics also necessitates matching Re_b to be in the range of tiny flying insects. When both Re_c and Re_b are matched between a physical bristled wing model to those of tiny insects, the scale model will produce non-dimensional forces similar to those of real insects. This is the major reason for presenting forces in term of non-dimensional coefficients throughout this study.

It has been reported that thrips (Kuethe, 1975) and the wasp Encarsia formosa (Ellington, 1975) operate at Re_b=10^–2 and 10⁻¹, respectively, and both at Re_c∼10. With the exception of Jones et al. (2016), the majority of modeling studies of bristled wing aerodynamics (Sunada et al., 2002; Santhanakrishnan et al., 2014; Lee and Kim, 2017; Lee et al., 2018; Kasoju et al., 2018; Ford et al., 2019) only matched Re_c∼10 without matching Re_b to be relevant to tiny insects. Matching Re_b ensures that the flow through bristles of a model (and hence Le) would be similar to that of real insects. Considering that lift and drag are known to be impacted by the extent of leaky flow (Kasoju et al., 2018), we matched Re_b to fall within 0.01 to 0.1 in the majority of our physical models.

Previous studies proposed that the substantial drag reduction realized with bristled wings in clap-and-fling is due to fluid leaking through the bristles (Santhanakrishnan et al., 2014; Jones et al., 2016; Kasoju et al., 2018). We found that Le peaked at τ∼0.56 or τ∼0.63 (Fig. 7B) for each condition of varying G and D, corresponding to the beginning of the fling phase. Interestingly, both C_D,max and C_L,max were observed between the same two time points, showing the importance of Le on dimensionless aerodynamic forces.

Previous studies of flow through bristled appendages found that Le is a function of both G and D (Cheer and Koehl, 1987; Hansen and Tiselius, 1992; Leonard, 1992; Loudon et al., 1994; Koehl, 1995). These studies also found that Le can change drastically for Re_b between 0.01 and 0.1, which is in the range of Re_b for tiny insects. We calculated Re_b for each wing model using D as the length scale in Eqn 2. Within the biological Re_b range (0.01–0.1), average force coefficients (⁠⁠, ⁠) showed no variation when varying D (Fig. 9A,B). For experiments of varying G, we maintained D and S as constants. The calculated Re_b was identical for these tests and within the biological Re_b range. Therefore, for a constant Re_b, can be varied significantly by varying G while maintaining minimal changes in (Fig. 9A,B).

Increasing Re_b via varying D showed opposite trends in C_D,max and Le_max (Fig. 9E,G). Within the biological Re_b range, increasing D decreased Le_max and increased C_D,max. Similarly, for a constant Re_b, increasing G increased Le_max and decreased C_D,max. These changes in leakiness for varying G and D are in agreement with previous studies (Cheer and Koehl, 1987; Loudon et al., 1994). Collectively, for Re_b in the range of tiny insects (0.01–0.1), we found that varying G provides drag reduction (C_D,max and ⁠) as compared with varying D, by augmenting Le. Tiny insects could possibly meet their flight demands by modulating the inter-bristle gap. Ellington (1980) observed that dandelion thrips (Thrips physapus) open their forewing setae prior to take-off, suggesting that modulation of G may be possible when preparing for flight.

Little to no variation in for both conditions (varying G and D) is attributed to formation of shear layers around the bristles that lowers the effective gap, resulting in the bristled wing behaving like a solid wing (Lee and Kim, 2017; Kasoju et al., 2018). Miller and Peskin (2005) proposed that LEV–TEV asymmetry plays a critical role in lift generation in clap-and-fling at Re_c∼10. For varying G and varying D, we observed LEV circulation (Γ_LEV) to be larger compared with TEV circulation (Γ_TEV) for most of the clap-and-fling cycle (Fig. 8B,C). The implication of this asymmetry on lift generation can be seen by examining time-variation of C_L (Fig. 4Bi,ii), where positive C_L was observed for most of the cycle. Both Γ_LEV and Γ_TEV peaked at the same time point at which we observed peak C_L.

Several studies examining the aerodynamic effects of varying S have reported contradictory findings. Although some studies found little variation in force coefficients (Usherwood and Ellington, 2002; Luo and Sun, 2005; Garmann et al., 2013), others have postulated that longer wingspans are detrimental for force generation (Harbig et al., 2013; Han, Chang and Cho, 2015; Bhat et al., 2019). All these studies considered solid wings at Re_c>100. Our study is the first to report the effect of varying S on the aerodynamic performance of bristled wings performing clap-and-fling at Re_c=10. Within the biological Re_b range, both and were found to increase with S (Fig. 9A,B). In addition, C_D,max and Le_max increased with increasing S (Fig. 9E,G).

The increase in G when increasing S is expected to increase Le and lower drag. However, we found that increasing S increased both Le and drag. Increasing S increases the wing surface area, which can explain the increase in drag. In addition, increasing G also increases Le. We speculate that the increase in Le with increasing S would minimize the increase in drag that would be expected from increasing wing surface area. Separately, varying S showed little change in Γ_LEV and Γ_TEV (Fig. 8D), which resulted in small changes in C_L (Fig. 4Biii). Within the biological range of n, G/D and Re_b, we postulate that larger S might be particularly beneficial to tiny insects when parachuting (Santhanakrishnan et al., 2014), as larger drag can slow their descent.

substantially increased with increasing n, while showed minimal variation for n≤88 and then increased with further increases in n (Fig. 6A). Wing models with n≤88 showed better aerodynamic performance in terms of force generation as compared with models with n>88. Interestingly, forewing morphological analysis showed that values of n were concentrated in the region of 30–90 for thrips and fairyflies. Moreover, generated for this dominant range of n was larger than generated for n=6 and 16. Thrips have been observed to intermittently parachute (Santhanakrishnan et al., 2014), likely to lower the energetic demands of flapping flight and potentially also during wind-assisted long-distance dispersal (Horridge, 1956). During parachuting, these larger drag forces can assist them in migrating longer distances (Morse and Hoddle, 2006). In addition, our morphological measurements showed that n varied from 32 to 161 across species, so lower n may better assist in generating lift needed for active flight, whereas larger n may better generate drag needed for passive dispersal via parachuting. Currently, it is unknown whether species with larger n tend to parachute more often.

Large variation in C_D,max and Le_max with n (Fig. 9F,H) showed the influence of the number of bristles on aerodynamic performance. Le_max decreased with increasing n, while C_D,max increased with increasing n. This suggests that changing n can aid or hinder aerodynamic performance by altering the leaky flow through the bristles. However, within the biological range of Re_b and n, only marginal changes in in comparison with were observed (Fig. 9C,D). This suggests that for a fixed S and G/D, tiny insects may experience reduced biological pressure to fit a particular number of bristles for adequate lift generation. This inference is also supported by the broad interspecific variation in n (Fig. 1C).

Within the biological Re_b range, C_D,max and Le_max were found to minimally change with increasing G/D (Fig. 9F,H). Also, varying G/D within the biological Re_b range produced little to no variation in and ⁠. Note that for varying G/D within the biological Re_b range, the inter-bristle gap in the corresponding physical models was nearly identical, which likely explains the minimal change in Le_max. From these results, we summarize that within the biological range of Re_b, G/D variation for a fixed S, n and G results in little variation in aerodynamic force generation.

Morphological measurements showed that G/D in thrips decreased with increasing BL, while the relationship was shallower for fairyflies. This dissimilar result in fairyflies and thrips raises a question regarding our use of static wing images for G/D measurements as opposed to free-flight wing images. We were restricted to using static forewing images owing to the lack of free-flight wing images of tiny insects with adequate (i.e. high) magnification. It is unknown at present whether tiny insects can modulate G/D during free flight, as such a capability could permit them to tailor aerodynamic forces in relation to ambient conditions (e.g. temperature, humidity, wind speed) and associated energetic costs.

We see many directions for future work. First, many bristle-winged insects show asymmetry in wing shape (Fig. 1; Jones et al., 2016). We did not consider the effects of the asymmetry in L_b on either side of the forewing (i.e. leading edge and trailing edge) or of bristle angle relative to the horizontal. Asymmetry in L_b within the biological Re_b range may not noticeably affect clap-and-fling aerodynamics, because damage may occur to the wing bristles during an insect's life and biological systems are often robust to such perturbations. Nonetheless, this may be a worthwhile direction for future work. Similarly, our physical models did not account for variation in wing shape and were simplified to a rectangular planform. There is much additional diversity in wing shape, especially when comparing fairyflies (teardrop-shaped) with thrips (smaller chord relative to span; Ford et al., 2019). At Re_c=10, changes in wing shape did not significantly affect the trend of aerodynamic force generation in time during clap-and-fling (comparing lift and drag force generation of rectangular bristled wing pairs used in Kasoju et al., 2018 to approximately elliptical bristled wing pairs used in Ford et al., 2019). However, the possible effects of wing shape on flying in bristle-winged insects – particularly across body sizes – would be valuable to study. Finally, the bristles on the wings of these insects are considerably flexible, yet we suspect them to behave stiffer in motion owing to high viscous forces. This was also evident with the stainless-steel wires that we used as bristles. Although these wires looked very flexible in air, the wires did not flex when tested in glycerin. We chose bristles that did not flex during motion because no quantitative data are available on flexibility of bristles in tiny insects. Based on published high-speed video of thrips (Santhanakrishnan et al., 2014; Cheng and Sun, 2018; Lyu et al., 2019), it is evident they flex their wings along the spanwise direction when flinging their wings apart at the start of downstroke. Because the variability in the wing flexibility along the wingspan has not yet been characterized in any published study, we used rigid wing models. Future studies are needed to document interspecific diversity in wing shape and flexibility to examine how they might affect aerodynamic forces.

Our analysis of forewing morphology in thrips and fairyflies showed similar scaling relationships between the two groups in the variables tested (n, G/D and S_max). Within the biologically relevant range of Re_b (0.01–0.1) for tiny insects, we observed that increasing the inter-bristle spacing (G) for fixed bristle diameter (D) decreased drag forces significantly. This was supported by a significant increase in leakiness observed during early fling. However, changes in average lift forces were minimal, suggesting that having the capability of increasing the inter-bristle spacing during free flight could help these insects to overcome large drag forces with minimal changes in lift force. We also found that varying bristle diameter (D) had no effect on aerodynamic force generation, and varying the non-dimensional inter-bristle gap-to-diameter ratio (G/D) showed no significant influence on aerodynamic force generation. Finally, although we found that drag forces significantly decreased with decreasing number of bristles (n), lift force only minimally changed for n<100. At n>100, we observed a significant jump in lift forces. Considering the broad variation of n (32–161) observed across species, the lack of change in lift forces for n<100 suggests that tiny insects may experience less biological pressure to optimize n for a given wingspan. Alternatively, stabilizing selection may maintain species within a range of values that does not affect flight performance.

The authors would like to thank the monitoring editor Sanjay Sane and three anonymous reviewers for their constructive comments.