Maximum whole-body force production can influence behavioral outcomes for volant taxa, and may also be relevant to aerodynamic optimization in microair vehicles. Here, we describe a new method for measuring maximum force production in free-flying animals, and present associated data for the wandering glider dragonfly. Flight trajectories were repeatedly acquired from pull-up responses by insects dropped in mid-air with submaximal loads attached beneath the center of body mass. Forces were estimated from calculations of the maximum time-averaged acceleration through time, and multiple estimates were obtained per individual so as to statistically facilitate approximation of maximum capacity through use of the Weibull distribution. On a group level, wandering glider dragonflies were here estimated to be capable of producing total aerodynamic force equal to ∼4.3 times their own body weight, a value which significantly exceeds earlier estimates made for load-lifting dragonflies, and also for other volant taxa in sustained vertical load-lifting experiments. Maximum force production varied isometrically with body mass. Falling and recovery flight with submaximal load represents a new context for evaluating limits to force production by flying animals.
Limits to maximum flight performance may derive independently from anatomical, physiological and aerodynamic constraints on animals, or from their interactions. Experimental manipulations of wing shape (e.g. Vance and Roberts, 2014; Ray et al., 2016) and wing flexibility (e.g. Mountcastle and Combes, 2013), along with hypobaric and hypodense gas manipulations (e.g. Chai and Dudley, 1995; Chai and Dudley, 1996), have demonstrated various limiting factors on force and power production for a number of volant taxa. Such studies may also inform engineering applications, e.g. optimization of aerodynamic forces generated by the wings of microair vehicles. Maximum whole-body force was recognized early on as potentially indicative of limits to animal flight performance (see Plateau, 1865), and has often been studied via vertical load-lifting assays given that the majority of resultant forces during flight act to offset gravity (see Marden, 1987; Dudley, 2000).
Dragonflies (Odonata: Anisoptera) are fast and agile insect fliers, and are known to exhibit free-flight accelerations as high as ∼40 m s−2 (Bomphrey et al., 2016; May, 1991; Rüppell, 1989; Lohmann et al., 2019). Marden (1987) cumulatively applied abdominal loads to 29 individual dragonflies from six different species, and derived maximum vertical forces averaging about 2.8 times the body weight. Methodological constraints associated with cumulative load lifting may, however, have resulted in underestimates of the maximum capacity for force production (see Dudley, 2000; Buchwald and Dudley, 2010). An alternative asymptotic load-lifting method was developed by Chai et al. (1997) for hummingbirds and later modified for orchid bees (Dillon and Dudley, 2004), bumblebees (Buchwald and Dudley, 2010) and tree sparrows (Sun et al., 2016), and has enabled detailed studies of allometric and ecological variability in the maximum vertical forces produced by animal fliers. However, our preliminary application of the asymptotic load-lifting method to dragonflies yielded erratic behaviors and non-vertical flight trajectories.
Instead, here we applied a newly developed experimental protocol to the wandering glider dragonfly (Pantala flavescens), a species renowned for extraordinary feats of flight, including transoceanic migrations (Anderson, 2009) and occurrence at the particularly high elevation of 6300 m in the Himalayas (Corbet, 2004). The experimental protocol involved the dropping in mid-air of dragonflies to which were attached a range of submaximal loads, followed by tracking of the ensuing trajectory. We then estimated time-averaged aerodynamic forces produced by the insects during this loaded fall and their subsequent ascent using acceleration estimates as extracted from the flight trajectory. We hypothesized that sustained forces elicited from these lifting trajectories would exceed those obtained in earlier cumulative load-lifting evaluations of dragonflies (Marden, 1987), and also that free-flight responses to load (and potential motivating cues of optomotor slip) might elicit the generation of weight-specific forces exceeding those recorded for other insects of comparable body mass in asymptotic load lifting. We also obtained repeated measures of performance per individual, and used Weibull distribution modeling (Hagey et al., 2016) to better estimate individual lifting capacity.
MATERIALS AND METHODS
Dragonflies, Pantala flavescens (Fabricius 1798), were captured on the main campus of Beihang University (Beijing, Peoples Republic of China) from 3 to 10 August 2018, and were used in experiments on the day of capture. The male/female ratio of the study population was approximately 3:2, but we pooled data from the two sexes to obtain a species-level assay of maximum performance. Following hypothermic anesthetization of individual dragonflies (Li et al., 2018), body mass and wing morphological parameters (Table S1) were measured. For each individual dragonfly, a small globule of melted soldering tin was then glued onto the ventral segment of the metepimeron nominally beneath the center of body mass of the insect (Fig. 1; Li et al., 2018). The mass of each globule (ml) was specified as a ratio relative to the initial mass of the dragonfly (mb); values for this loading ratio R=ml/mb of 1.0 (N=1), 1.2 (N=3), 1.5 (N=19) and 1.8 (N=7) were used in experiments, and were assigned haphazardly to individual dragonflies. Loads greater than 1.8 times the body mass were not used because preliminary trials (N=5) found no individual that could offset a load twice their body mass during loaded fall after being dropped from rest. In experiments, each dragonfly (within several minutes of morphological measurements) was held by its four wings as folded dorsally, was lifted near the top of the flight chamber, and was then dropped at a nominally horizontal body orientation so as to initiate the loaded fall (Fig. 1; Movies 1 and 2). For filming, only individuals with intact wings and no obvious body damage were used. Two orthogonally arranged cameras (operated at 500 frames s−1) were used to film flight trajectories within an open-top glass cube (volume of 0.6 m3; 0.8 m×0.8 m×1.0 m), using a configuration identical to that of Li et al. (2018). Vertical coverage of the two cameras was restricted to regions at least 0.3 m above the cube's floor, thereby excluding filming of any flights that potentially would invoke the ground effect. Following initial filming, additional trials were performed (with 3 min resting periods between flight events) up to the time when insects failed to sustain body weight during flight. During filming, air temperature inside the experimental chamber was ∼36°C. Before filming and in between trials, dragonflies were housed in mesh cages at ambient air temperature, which was typically 30°C.
time-averaged X-accelerations during forward phase
time-averaged Y-accelerations during forward phase
time-averaged Z-accelerations during pull-up phase
degrees of freedom
horizontal aerodynamic force produced accompanying the production of FV
body weight-specific horizontal aerodynamic force accompanying the production of FV
individual maximum aerodynamic force production
individual maximum body weight-specific aerodynamic force production
maximum vertical aerodynamic force produced during a single flight
maximum body weight-specific vertical aerodynamic force produced during a single flight
- F*X, F*Y, and F*Z
X-, Y- and Z-component of body weight-specific aerodynamic force, respectively
gravitational acceleration (−9.8 m s−2)
initial body mass of dragonfly
mass of globule (load)
number of individuals/sample size
time after release
time period of dragonfly's acceleration from wmin to wmax
- u, v and w
X-, Y- and Z-velocity, respectively
dragonfly's maximum vertical velocity after upward acceleration
dragonfly's minimum vertical velocity during the flight
- x, y and z
X-, Y- and Z-position, respectively
- (X, Y, Z)
global coordinate system
scale parameter of Weibull distribution model
Pixel coordinates approximating the position of the dragonfly's center of mass during flight were tracked using DLTdv5 digitizing tool (Hedrick, 2008) operating in Matlab (126.96.36.1992, R2014a). Using the same three-dimensional reconstruction algorithm as in Li et al. (2018), time series of the pixel coordinates were converted to a global (X, Y, Z) coordinate system in which the positive direction of the Z-axis was defined as vertically upward (see Li et al., 2018), the Y-axis was generally aligned with the dragonfly's sagittal plane anteriorly when initially dropped, and the X-axis pointed in the right lateral direction of the dragonfly. Positional data were then smoothed using a quintic spline filter with the low-pass cut-off frequency specified as 20 Hz, which was approximately the lower limit of dragonfly's wing flapping frequency during the loaded fall. Root mean square level of the removed noise was around 0.7 mm, and smoothed positional data were differentiated to produce time series of velocity vectors (see Fig. S1). The aforementioned data processing procedure was completed using the GCVSPL package (https://isbweb.org/software/sigproc.html; Woltring, 1986). We then identified sufficiently long time periods during which acceleration was relatively constant, here termed quasi-constant acceleration (QCA). Two criteria were used to identify periods of QCA; one was the absolute value of the Pearson correlation coefficient between a velocity component and time over sample periods, for which we set the high-pass threshold to be 0.99. As revealed by the flight trajectory kinematics (see Results and Table 1), dragonflies produced maximum aerodynamic force over a time period of approximately 150 ms. Therefore, we chose a minimum period length of 150 ms (corresponding to 3–5 wingbeats; see Li et al., 2018; Wakeling and Ellington, 1997); both aforementioned criteria were necessary for classification as a QCA period (see Fig. 2). All such periods used for analysis corresponded to at least 76 video frames in duration. Linear equations were fitted using ordinary least squares regressions to the velocity–time curves within QCA periods. The slope of each regression was used as an estimate of the time-averaged acceleration for the corresponding QCA period; time-averaged accelerations for all QCA periods were then compared to obtain the dragonfly's maximum time-averaged acceleration along a coordinate axis during the entire flight (Fig. 2).
We then estimated group-level maximum performance from multiple observations on different individuals (see Hagey et al., 2016, 2017). First, individual maximum capacity based on repeated measures of maximum performance was determined, using Weibull distribution modeling, which substantially reduces the influence of rare observations (see Hagey et al., 2016). For individuals for which maximum force production was evaluated multiple times, maximum capacity was equated to the scale parameter (λ) of the Weibull distribution model, and the estimation error was equated to the standard error of λ. For individuals with only one observation, individual maximum capacity was assumed to equal that of this single performance. Species-level maximum capacity was then estimated as the average of multiple individual capacities (as weighted by their corresponding estimation errors), thus minimizing effects of either small sample sizes for particular individuals or large variation among trials (Hagey et al., 2016). Data from individuals with only one observation were excluded when estimating species-level capacity. Weibull distribution modeling was performed using Matlab (188.8.131.522, R2014a); other statistical tests, e.g. Student's t-tests, were performed using R (3.5.2).
Thirty individuals were filmed in experiments, from which complete trajectories were obtained for 143 loaded falls (referred to as full-trajectory trials; see Fig. 3; Fig. S2 and Table S2). An additional 35 trials were filmed which omitted the initial descent, referred to here as partial-trajectory trials (Table S2). Data from partial-trajectory trials were not included in the graphical illustration of the results (i.e. Figs 3 and 4), but were used to evaluate maximum acceleration. Typically, falling individuals locomoted downward and forward as in diving flight, followed by a pull-up (Fig. 3). Dragonflies veered haphazardly to the left and right when falling (Figs 3 and 4A), but showed generally consistent trends in downward and forward motions (Fig. 4C,E) and in associated speeds (Fig. 4D,F).
Four different temporal phases of the trajectory were identified based on overall trends in the velocity curves. A free-fall phase (0–0.10 s after the drop; phase I in Fig. 4) was marked by rapid increases in the dragonfly's falling speed, but with relatively small variation in horizontal speed. Right after the drop, the falling speed of the dragonfly increased rapidly to 1–2 m s−1 (Fig. 4F). From 0 to 0.10 s after the drop, the vertical velocity curves exhibited a mean slope of −7.9 m s−2, suggesting that acceleration during this phase was mainly due to gravity with no significant aerodynamic force production.
Then, a forward phase of the flight (0.15–0.25 s; phase II in Fig. 4) was marked by a rapid and consistent increase in the dragonfly's forward velocity (up to ∼2 m s−1; Fig. 4D), accompanied by an increase in the vertical acceleration, as indicated by convexity in the vertical velocity curves (Fig. 4F). X-velocity curves exhibited a generally upward trend, but with substantial variation among separate trials. During this forward phase, dragonflies in most trials reached their minimum vertical speed. From 0.15 to 0.25 s after the drop, the averaged X-velocity and Y-velocity curves exhibited a mean slope of 2.1 m s−2 and 8.6 m s−2, respectively. Dragonflies maintained a high horizontal force in this phase to effect forward acceleration, along with increased vertical force production.
During a pull-up phase (0.3–0.45 s; phase III in Fig. 4), vertical velocity rapidly increased concomitantly with reductions in horizontal acceleration. Vertical velocity curves also showed reduced convexity (Fig. 4F), suggesting positive and more constant accelerations compared with those in the forward phase of the flight. The mean slope of the averaged Z-velocity curve from 0.3 to 0.45 s after the drop was 6.8 m s−2, while vertical velocity approached zero (Fig. 4F). Dragonflies thus produced high and sustained vertical forces in parallel with decreased horizontal force.
In the final recovery phase (0.50–0.55 s; phase IV in Fig. 4), when the falling speed of the dragonfly approached zero, vertical and horizontal accelerations simultaneously decreased relative to those in the pull-up phase. From 0.50 to 0.55 s after the drop, the mean slope of the averaged Y-velocity and Z-velocity curves decreased to −5.6 m s−2 and 2.7 m s−2, respectively, indicating reduced vertical force output and slowed forward motion.
In some trials, dragonflies exhibited time-averaged X-accelerations as high as 9 m s−2 during the forward phase (see Fig. 4B, green arrow). For all full-trajectory trials, we calculated time-averaged X- and Y-accelerations during the forward phase, denoted as ax2 and ay2, respectively. The corresponding magnitudes of horizontal acceleration, , were all less than 13.5 m s−2 (Fig. 5). Also, there was a complementary relationship between the square of ax2 and that of ay2, as indicated by a significant linear regression between these data points whenever the magnitude of horizontal acceleration exceeded 9 m s−2 (see Fig. 5). This result suggests that, rather than being associated with an independent lateral motion, X-accelerations simply resulted from a reorientation of the total force vector.
To summarize the general features of these loaded fall trajectories, dragonflies underwent large anterior accelerations following an initial period of free fall. After rapid anterior acceleration, falling speed then tended to decline as a result of increased and sustained vertical force production.
Maximum force production in single flights
The magnitude and direction of time-averaged aerodynamic force production during the four phases were first estimated from slopes of the averaged velocity curves during corresponding time intervals. For a loading ratio (R) of 1.5 (the value most frequently used in this study), body weight-specific force components and the resultant total force during the four phases were calculated (Table 1). Time-averaged Z-acceleration during the forward phase, and X- and Y-accelerations during the pull-up phase, were assumed to be zero, given convexity in the associated sections of the velocity curves (see Fig. 4D,F).
Moreover, to account for the aforementioned discrepancy between individual X-velocity curves and the averaged X-velocity curve during the forward phase of some trials, we compared time-averaged aerodynamic force production during the forward phase with force production during the pull-up phase for each separate flight. We calculated time-averaged Z-acceleration (az3) during the pull-up phase for 102 full-trajectory trials (duration of the other 41 trials being less than 0.45 s), and obtained each trial's body weight-specific resultant force production during the two phases, (1+R)/|g|·‖(ax2,ay2, − g)‖ and (1+R)/|g|·‖(0,0,az3, − g)‖, respectively. A paired bi-directional Student's t-test (P<0.001) demonstrated that dragonflies significantly generated 15% more body weight-specific force (difference of mean body weight-specific force of 0.49, 95% confidence interval: 0.41–0.58) during the pull-up phase relative to the forward phase.
Consequently, the greatest resultant force was produced in the pull-up phase (among the four trajectory phases), and its direction was generally aligned vertically (Table 1). Using the QCA method, we calculated maximum vertical acceleration for 148 trials derived from 30 individuals (consisting of 121 full-trajectory trials from 25 individuals, and 27 partial-trajectory trials from 11 individuals; partial-trajectory trials contained the pull-up phases for maximum vertical acceleration calculations). The other 30 trials of the 178 total trials lacked appropriate QCA periods to calculate maximum force and were not included in the analysis (see Discussion). Corresponding maximum vertical forces FV averaged 1.41×10−2 N (s.d. of 3.0×10−3 N), and their body weight-specific values F*V averaged 4.18 (s.d. of 0.58), with the highest value being 5.81. Horizontal forces produced during the pull-up phase could not be accurately calculated using the QCA method because of considerable variation in X- and Y-accelerations (see Fig. 4B,D). Nonetheless, horizontal forces were approximated for each trial by time-averaging components of horizontal acceleration over the same QCA period. Horizontal forces FH during these periods averaged 3.3×10−3 N (s.d. of 1.9×10−3 N), and body weight-specific values F*H averaged 0.99 (s.d. of 0.52), corresponding to ∼24% of the vertical force and with the highest value being 2.54. There was a non-significant linear relationship between horizontal and vertical forces (F*H = 0.13 F*V+0.43, r2=0.02, P=0.07).
Combining these maximum vertical force estimates with the corresponding horizontal forces, the total aerodynamic force was calculated as their resultant. For the 148 trials, total force exceeded the maximum vertical force only by small amounts (average of 4%, with a maximum of 18%). Given their small contribution to total force as well as their much lower accuracy of estimation, all horizontal forces were subsequently neglected and the vertical force was considered to represent the total force output.
Maximum force production by individuals, and group-level capacity
Neglecting innate individual differences, two external factors (i.e. the loading ratio and the categorical variable as to whether or not the start of the loaded fall trajectory was recorded) were considered to evaluate associated effects on estimates of individual maximum performance. Individual maximum body weight-specific force production (F*max), as determined using Weibull distribution modeling for multiple-observation individuals, was used to compare individuals of different body mass.
For the 143 full-trajectory trials, maximum vertical force production was obtained using the QCA method from 121 trials by 25 individuals (multiple-observation individuals, N=18), for which F*max was then determined. These values of F*max were then segregated into three groups according to their loading ratio [i.e. R=1.0 or 1.2 (N=4), R=1.5 (N=15) and R=1.8 (N=6)], and unpaired bi-directional Student's t-tests were employed to detect differences among groups (Table 2). F*max of individuals with an R of 1.5 did not significantly differ from those with an R of 1.8, whereas individuals with an R of 1.0 or 1.2 exhibited significantly lower values of F*max when compared with the other two groups with higher R (see Table 2). This result implied that individuals with relatively low values of R (i.e. 1.0 or 1.2) may not have exhibited maximum performance and, consequently, 25 full-trajectory trials from these individuals were excluded. These individuals did not produce any partial-trajectory trial.
For the six individuals that produced at least one full-trajectory trial and at least one partial-trajectory trial, F*max was estimated twice for each individual using either only full-trajectory trial(s) or only partial-trajectory trial(s). A paired bi-directional Student's t-test demonstrated that F*max from these two categories did not significantly differ (see Table 2), indicating that whether or not the start of the loaded fall trajectory was recorded did not affect the estimate of individual maximum force production.
Based on the 123 observations (96 full-trajectory trials and 27 partial-trajectory trials) for maximum force production from 26 individuals [R=1.5 (N=19) and R=1.8 (N=7)], the maximum value of individual maximum force production (F*max) for this set of individuals was 1.89×10−2 N (s.e.m. of 3.5×10−4 N), whereas the maximum F*max within the sample was 5.47 (s.e.m. of 0.157). Multiple observations of maximum vertical force were obtained for 23 dragonflies [R=1.5 (N=17) and R=1.8 (N=6)]. Weighting the F*max of these 23 individuals by their estimation errors, we then estimated maximum body weight-specific force production for the entire sample to be 4.27 (s.e.m. of 0.005).
After determining values of Fmax for each of the aforementioned 26 individuals, log-transformed data were regressed against log-transformed body mass to evaluate the allometry of force production (Fig. 6). For the pooled dataset individuals, the ordinary least squares linear regression was significant (logFmax=0.86logmb−1.44; r2=0.71, P<0.001), with 95% CI for the slope between 0.64 and 1.07, indicating non-significant deviation from isometry.
The method used here for eliciting maximum free-flight forces in dragonflies clearly indicates a capacity well exceeding that of steady trimmed flight, but identification of upper limits requires that (a) maximum performance during loaded falls and subsequent recovery were adequately assessed, and (b) the estimated maxima approximate absolute maximum capacities.
The 30 trials with no QCA period present, and thus with no observed maximum force production, are referred to as non-QCA trials; the other 148 trials with a QCA period present are referred to as QCA trials. Values of Cr for QCA trials ranged from 0.19 to 1.78, whereas values for non-QCA trials were all less than 0.67, indicating that the QCA method was sufficient for all trials with a Cr higher than 0.67. For non-QCA trials, the relatively low values of Cr suggest an inability to maintain large vertical force production over a sufficient period of time so as to gain enough positive vertical momentum. We also calculated, for all non-QCA trials, the time-averaged body weight-specific vertical forces produced during acceleration from wmin to wmax. Equating these values to maximum force production during flight, we obtained an estimate of F*max for 13 individuals. A paired bi-directional Student's t-test was then conducted to compare these F*max estimates with those estimates of the same individuals as determined from QCA trials; the former was significantly smaller than the latter, by a value of 0.57 (95% confidence interval: 0.31–0.84; see Table 2). Thus, non-QCA trials can be ignored as dragonflies did not maintain high aerodynamic forces during these flights.
Addressing the second criterion for identification of upper force limits, dragonflies generally extended the duration of vertical force production so as to achieve higher values of Cr, rather than increasing maximum forces per se. This result suggests that dragonflies were near their limits to force production in the recorded flights. For the 123 trials (96 full-trajectory trials and 27 partial-trajectory trials) used to determine group-level force production capacity, we obtained the time period T during which the dragonfly increased its vertical velocity from wmin to wmax, and investigated the relationship between T and Cr, and between F*V and Cr, respectively (Fig. 7). When Cr increased from 0.2 to 1.8, F*V increased by 26% (from 3.85 to 4.85; see Fig. 7A), whereas T increased by 300% (from 0.12 s to 0.48 s; see Fig. 7B). With increasing Cr, the duration of force production thus increased 10 times more than did the magnitude of these forces, suggesting that the latter was close to an upper limit.
Estimated group-level capacity of the wandering glider suggests a time-averaged maximum vertical force of 4.27 times the body weight. This value is substantially higher than that reported for dragonflies in cumulative lifting of abdominally applied loads (i.e. 2.8; see Marden, 1987), and is comparable to transient estimates for free-flying dragonflies in various behavioral contexts (see Bomphrey et al., 2016; May, 1991; Rüppell, 1989; Lohmann et al., 2019). The value estimated here also well exceeds peak values for flapping cycle-averaged aerodynamic forces during voluntary maneuvers of dragonflies [e.g. 1.54 in turning flight (Li and Dong, 2017) and 1.33 in take-off (Li et al., 2018)]. Surprisingly, the value obtained here is much higher than the body weight-specific maximum force capacity reported for vertical load lifting by a range of other flying animals (some of which are of comparable body mass), including a value of 1.52 for bumblebees (Buchwald and Dudley, 2010), 1.8–2.1 for orchid bees (Dillon and Dudley, 2004), 2.2 for fruit flies (Lehmann, 1999), 2.26 for tree sparrows (Sun et al., 2016) and 1.8–3.9 for hummingbirds (Chai et al., 1997; Altshuler et al., 2010). These comparisons suggest that free forward flight with submaximal load may elicit behavioral responses different from those in vertical ascent, and that kinematic constraints identified in vertical and sustained load lifting (e.g. maximum stroke amplitude; see Dudley, 2000) may not pertain to the flight behaviors described here. Detailed high-speed assessment of wing motions would be necessary to assess this hypothesis. The results obtained for the wandering glider also indicated an isometric increase in total force production with increasing body size (Fig. 6), as opposed to negative allometries obtained in interspecific studies of asymptotic load lifting in orchid bees (Dillon and Dudley, 2004) and hummingbirds (Altshuler et al., 2010). Further assessment of intraspecific variation in morphological (e.g. relative muscle mass) and kinematic features (e.g. wingbeat frequency) during loaded flights of dragonflies would be informative in understanding these differing allometries.
In conclusion, we present a new method for assessing maximum force production by free-flying dragonflies, namely submaximally loaded falling, which elicits much higher whole-body forces than have previously been reported. Using a subset of trajectories for which accelerations were both high and sustained (i.e. using the QCA method), and by applying the Weibull distribution to better estimate population-level performance, we estimate maximum aerodynamic force production to be 4.3 times the body weight for this particular dragonfly species. With corresponding modification of the present method, consistent measurements of maximum force production during free flight could be obtained for various volant taxa, some of which have been studied relative to performance limits only in hovering flight using load-lifting methods. The method described here is also amenable to studying flight with ablated wings, and flight in variable-density gas mixtures.
We thank Yingxin Du, Wenqian Wu and Chengwang Zhang for assistance with experiments and data analysis, Yu Zeng for helpful comments, Victor M. Ortega-Jimenez for instructions on GCVSPL, and Travis J. Hagey for advice on Weibull distribution modeling.
Conceptualization: G.S., R.D., T.P., Q.L.; Methodology: G.S., R.D., T.P., Q.L.; Software: G.S., M.Z., L.P.; Validation: G.S., M.Z., L.P.; Formal analysis: G.S., R.D., T.P., Q.L.; Investigation: G.S., M.Z., L.P.; Resources: G.S., R.D., T.P., M.Z., L.P., Q.L.; Data curation: G.S., T.P.; Writing - original draft: G.S., R.D., T.P.; Writing - review & editing: G.S., R.D., T.P., M.Z., L.P., Q.L.; Visualization: G.S., R.D., T.P.; Supervision: R.D., T.P., Q.L.; Project administration: G.S., R.D., T.P., M.Z., Q.L.; Funding acquisition: G.S., T.P., Q.L.
This work was supported by National Natural Science Foundation of China (grant nos 51636001 and 51706008), National Science and Technology Major Project (2017-II-0005-0018), and Aeronautics Power Foundation of China (grant no. 6141B090315). G.S. was funded by China Scholarship Council.
The authors declare no competing or financial interests.