Twist–torsion coupling in beating axonemes

Cilia are slender, membrane-enclosed protrusions of eukaryotic cells. Motile cilia drive essential physiological functions such as the breaking of the left–right symmetry in mammalian embryos¹, fluid transport in the human respiratory tract, fallopian tubes and cerebral cavities^2,3,4 as well as the locomotion of single-cell micro-swimmers such as sperm and algae. The internal mechanical core of cilia is evolutionarily conserved among eukaryotes and is called the axoneme. The axoneme consists of nine doublet microtubules (DMTs), which are cylindrically arranged around a pair of singlet microtubules (which scaffold the central apparatus) (Fig. 1a). Dynein motors are distributed along the entire length of the axoneme, arranged in a chiral fashion, and, when active, slide adjacent DMTs. This sliding is constrained at the ciliary base and along the axoneme, which converts sliding into travelling waves of bending that shape the ciliary waveform.

a, A schematic of the axoneme cross-section (viewed from the base) with DMTs numbered 1–9 (ref. ¹⁴). The bridge (dark red) that connects DMT1 and DMT2 may set the bending plane (the rigid-bridge hypothesis). The green and red arrows indicate the material frame of the axoneme (vectors ({mathbf{e}}_{1}) and ({mathbf{e}}_{2}) span the cross-sectional plane, with ({mathbf{e}}_{1}) pointing towards the bridge and ({mathbf{e}}_{2}) being orthogonal to ({mathbf{e}}_{1}), while ({mathbf{e}}_{3}) is normal to the cross-section, pointing along the centre line). Dynein motors (green) are permanently attached to one DMT and transiently interact with their clockwise neighbour DMT. b, A cartoon of a twisted axoneme in which torsion and twist are coupled. Left: torsion quantifies the rotation of the bending plane (rotating green plane) along the 3D centre line (blue). Right: twist quantifies the rotation of the material cross-section of the axoneme (the cross-section orientation is marked by a dashed dark-red line, indicating the bridge). c, The expected movement of an axoneme-attached GNP assuming twist–torsion coupling. In this scenario, the rotation of the bending plane (green) by an angle Δω_3D is exclusively due to axonemal twist, which causes an equal rotation Δω_GNP of the attached GNP. (Rotation angles are measured in a plane (light blue) normal to the x–y plane of the laboratory coordinate system.) d, The alternative scenario of twist-free torsion of the axoneme, where the rotation of the bending plane (green) is caused by bending into different directions relative to the bridge (dashed dark-red line). In this scenario, the rotation of the cross-section (with attached GNP) by an angle Δω_GNP is different from the rotation angle of the bending plane ({Delta omega }_{3{rm{D}}}). The planar axoneme (transparent grey) in c and d depicts a reference shape (with no torsion) of the axoneme and the reference position of the GNP. The insets in c and d show axonemal cross-sections at the indicated light-blue planes. The red and yellow spheres show the GNP in its reference position (planar shape) and its final position (non-planar shape), respectively. Note that the magnitude of deformations is exaggerated for display purposes. For an animated version of Fig. 1c,d, see Supplementary Movie 1.

Full size image

Ciliary waveforms are often non-planar. This non-planarity supports their physiological functions, for example, fluid pumping^5,6,7. In single cells, non-planar beat patterns cause cell rotation and helical swimming, which are key to chemo- and phototactic navigation strategies along chiral paths^8,9,10,11 as well as rheotaxis^10,12,13. How the non-planarity in ciliary waveforms is generated and how it relates to internal deformations of the axonemal structure is not known.

To investigate the relation between non-planarity and structural deformations in axonemal shapes, the concepts of twist and torsion need to be carefully distinguished: (1) twist describes the internal rotation of the axoneme along the arc length and characterizes the structural deformation of the axoneme’s material whereas (2) torsion describes the rotation of the bending plane along the arc length and is a mathematical property of three-dimensional (3D) centre line curves. When twist and torsion are strictly coupled, twist causes an equal amount of torsion for a bent object (Fig. 1b) and both terms are equivalent. This, for example, is the case for an object such as an ordinary plastic ruler that has a highly anisotropic bending stiffness, which sets a preferred bending plane (for example, it only bends in one but not the other direction). The axoneme, however, is a filament bundle for which the degree of anisotropy in bending stiffnesses is unknown, and thus this coupling is non-trivial.

The popular ‘rigid-bridge hypothesis’ asserts that, in Chlamydomonas, the bridge, an internal axonemal component that links DMT1 and DMT2, prevents sliding between these DMTs¹⁴. While the bridge is only present in the proximal part in Chlamydomonas cilia¹⁵, other structural heterogeneities, such as inter-doublet linkers (IDLs), probably also contribute to an anisotropic bending stiffness with equivalent functional implications¹⁵.

If the rigid-bridge hypothesis were true, any rotation of the bending plane of the axoneme (that is, torsion) would necessarily result from an internal rotation of the axoneme, that is, twist (Fig. 1c). On the contrary, if the sliding between DMT1 and DMT2 were fully unconstrained, the bending plane would be free to rotate. This would allow for a scenario in which the axoneme assumes a non-planar shape but does not twist (Fig. 1d). Thus, the question arises of whether twist and torsion are coupled and whether the rigid-bridge hypothesis is true.

To test for twist–torsion coupling, both quantities must be measured in the same system and compared to each other. Whereas several reports indicate torsion in beating cilia and axonemes^{8,9,16,17,18,19,20,21,22,23}, twist cannot be assessed in a straightforward way. Note that torsion measurements alone cannot prove the existence of twist without access to the internal structure of the axoneme. To date, little is known about the twist in beating axonemes. Previous accounts of axonemal twist mostly came from static samples. For example, electron microscopic measurements on Paramecium cilia revealed hetero-chiral twist of 15–20° µm⁻¹ (ref. ²⁴). The only dynamic measurement to date is an observation by Woolley et al.²⁵, where a mitochondrion attached to a beating quail sperm tail served as a tracer particle to follow rotations of the axoneme. Earlier experiments attaching beads to beating sea urchin sperm were limited in resolution and could not resolve rotations²⁶. Neither of these experiments quantified torsion, and thus they cannot test twist–torsion coupling.

To investigate whether twist and torsion are coupled, we use reactivated axonemes, purified from the green alga Chlamydomonas reinhardtii, as a model. We measure the 3D waveforms with defocused dark-field microscopy and use a novel beat-cycle averaging method to achieve high spatio-temporal precision. This allows us to measure torsion reliably. Furthermore, we introduce a rigorous error estimate, providing a space–time map of axonemal torsion inside a region of trust. To measure local rotations of the axonemal cross-section, we attach gold nanoparticles (GNPs) as tracers to the outside of axonemes. By comparing these local cross-section rotations to the local rotations of the bending plane, we show that twist and torsion are coupled. This result relates the geometry of the 3D centre line of the axoneme to its internal structural deformation. Our measurements are consistent with a wave of twist that travels base to tip along the axoneme. In future studies, these twist dynamics can be used to test models that predict the motor activity in the beating axoneme.

High-precision average 3D waveform of isolated axonemes

We measured the 3D shapes of reactivated axonemes, isolated from Chlamydomonas reinhardtii cells, using defocused high-speed dark-field microscopy^17,27 (Fig. 2a–c, Supplementary Fig. 1 and Supplementary Movie 2). First, we used a filament tracking software²⁸ to determine the axonemal centre line in two dimensions (2D), characterized as x, y coordinates along the arc length. To compute the z coordinate of the axonemal centre line at each arc-length position, we exploited the fact that, in dark-field microscopy, the apparent width of an object increases as it becomes defocused. Specifically, we determined the motion blur-corrected full width at half maximum (FWHM) of the dark-field signal of axonemes (measured perpendicular to the centre line; Fig. 2b–d) at each arc-length position and used the relation between the FWHM and the z coordinate as a calibration curve (determined from axonemes immobilized to a glass coverslip; see Fig. 2c,e and Supplementary Fig. 2 for more details). Using this calibration curve, we reconstructed the 3D shapes of reactivated axonemes (Fig. 2f) and obtained 3D waveforms, which comprise the periodic sequence of axonemal shapes. The tracking errors in the x, y and z coordinates were σ_x,y ≈ 7.3 nm and σ_z ≈ 35.6 nm, respectively (Supplementary Fig. 3).

a, A Chlamydomonas reinhardtii cell imaged with phase-contrast microscopy (scale bar 3 μm). b, A defocused dark-field microscopy snap-shot of a reactivated axoneme with an exposure time of 1 ms (scale bar 1 µm). Coloured lines indicate two positions where the axoneme is defocused (red) and where the axoneme is close to the focal plane (blue). c, A schematic of the defocused imaging. Images are recorded by focusing below the axoneme sample (the focal plane is shown as a dashed grey line). d, The measurement of the FWHM in intensity profiles along the coloured lines in b, with lower FWHM values for axoneme parts closer to the focal plane. e, The calibration of the relationship between the FWHM and the distance to the focal plane z. Values were obtained through z scans of axonemes immobilized to the chamber surface closest to the objective (where positive z values correspond to axonemal positions above the focal plane). The calibration curve (green line) was obtained by a smoothing spline fit. The red and blue crosses correspond to the curves and lines in panel d and b. f, An example axoneme shape in 3D, with the z coordinate obtained using the calibration curve in e.

Source data

Full size image

Three-dimensional waveforms from individual axonemes were highly reproducible and thus allow to further increase the precision by averaging. We therefore defined a beat-cycle phase ϕ for each axonemal shape by fitting a periodic model function (Methods and Supplementary Fig. 4a–d). We then averaged the 3D shapes of 17 axonemes (with a total of 3,755 beat cycles and 14 frames per beat cycle on average) with similar phase ϕ, so that we obtained a highly precise average waveform with increased temporal resolution (32 shapes per beat cycle) (Fig. 3a,b). The positional uncertainty of shapes in this average 3D waveform were ({bar{sigma }}_{{xy}}) = 0.19 nm and ({bar{sigma }}_{z}) = 2.20 nm (standard error of mean (s.e.m.); Supplementary Fig. 4e–k), which represent excellent spatio-temporal precision for a beating axoneme or cilium.

a, The x–y projection of the average 3D waveform (an average over 17 reactivated axonemes with a total of 3,755 beat cycles), aligned at the basal position. The colour wheel represents the beat-cycle phase ϕ of each shape. The scale bars in x and y are 1 μm. b, A side view of the average 3D waveform (x–z projection, a rotated by 90° around the x axis with scale bars in x and z of 500 nm (and the z positions enlarged disproportionally), while the z-axis is normal to the boundary plane and points into the observation chamber). Inset: the non-planarity of the waveform, computed as the sum of the squared residuals between each shape of the average waveform and a fitted plane. c, The computation of the torsion τ_n and the 3D curvature к_n from the Frenet–Serret frame with a binormal vector (mathbf{b}) (red), normal vector (mathbf{n}) (green) and tangent vector ({mathbf{t}}) (blue) at subsequent arc-length positions (black filled circles) enumerated by n along the 3D centre line (black line), using the rotation angle of the local bending plane ω_τ and the in-plane rotation angle ψ of the tangent. Additionally, we measure ω_3D as the orientation of the normal vector with respect to the x–y plane of the laboratory frame. d, A map of the 3D curvature к as a function of the beat-cycle phase ϕ and the arc length s. e, A map of the torsion τ as a function of the beat-cycle phase ϕ and the arc length s (where red shows dextral and blue shows sinistral torsion). f, A map of the estimated error σ_τ of torsion as a function of the beat-cycle phase ϕ and the arc length s (calculated using bootstrapping; see Supplementary Fig. 5a–f for details). In d–f, hatched regions indicate where the absolute curvature is below 0.4 rad μm^–1. The complementary, non-hatched region defines a region of trust for estimated torsion.

Source data

Full size image

Non-planarity and torsion of the average 3D waveform

The average 3D waveform shows the strongest non-planarity during the recovery stroke (Fig. 3b; 0.5π < ϕ < 1.5π), which we quantified using a phase-dependent non-planarity measure (Fig. 3b, inset). We further characterize the non-planar centre line in terms of curvature and torsion. To do this, we use the Frenet–Serret frame, which is defined by three orthonormal vectors at each point along the arc length: the tangent vector (mathbf{t}), the normal vector (mathbf{n}) and the binormal vector (mathbf{b}) (Fig. 3c). The 3D curvature κ equals ψ/Δs, where ψ is the angle between subsequent tangent vectors (it mathbf{t}_{n}) separated by an arc-length increment Δs. The torsion τ equals ω_τ/Δs, where ω_τ is the angle between subsequent binormal vectors (mathbf{b}_{n}) separated by an arc-length increment Δs. Non-zero torsion indicates a rotation of the bending plane along the centre line (see also Fig. 1). With this definition, κ and τ have well-defined signs (up to a global choice of reference, see Methods).

Mathematically, the torsion τ is only defined at arc-length positions where the curvature κ is non-zero. In practical terms, this implies that torsion estimates become unreliable when the curvature is close to zero. To identify regions where we can reliably determine torsion, we locally quantify the variability of torsion estimates as a function of the absolute curvature |κ| (Supplementary Fig. 5a–f). We define a ‘region of trust’ for reliable torsion estimates for curvature values above |κ| > 0.4 rad μm⁻¹ where the estimated torsion error σ_τ is 4.8° μm⁻¹ or less (Fig. 3d–f, non-hatched region in (ϕ, s) space). From here on, we only consider torsion values within this region of trust. Torsion can be either negative (blue, sinistral, that is, left-handed) or positive (red, dextral, that is right-handed) and ranges from –24.8° μm⁻¹ to +22.7° μm⁻¹. At a given arc-length position, the torsion changes dynamically and occasionally even switches sign during one beat cycle (Fig. 3e and Supplementary Fig. 5g–i). This becomes especially apparent at arc-length positions s = 0.25L or 0.5L. To characterize dynamic torsion, we estimate the peak-to-peak amplitude of the phase-dependent torsion for different arc-length positions. We find that this amplitude is approximately constant along the arc length with an average of 20.0 ± 7.6° μm⁻¹ (mean ± s.d., N = 24; Supplementary Fig. 6b). Although this value only provides a lower bound as we only consider torsion within the region of trust, it shows that torsion changes with the beat-cycle phase, which means that torsion is dynamic. For a given beat-cycle phase, torsion also changes as a function of the arc length (Fig. 3e and Supplementary Fig. 5j–l). For example, at ϕ = 0.5π and ϕ = π, the sign and magnitude of torsion changes along the arc length. Generally, we observe that the region of trust in (ϕ, s) space is traversed by a minimum of torsion (Fig. 3e, blue, sinistral), reminiscent of a travelling torsion wave propagating from base to tip.

Cross-section rotation measured using GNPs

To test whether a rotation of the local bending plane (torsion) is accompanied by a rotation of the axonemal cross-section (twist), we first determined both rotations in a common frame of reference, the laboratory frame. There, a rotation of the local bending plane is characterized by an angle ({Delta omega }_{3{rm{D}}}) and a local rotation of the axonemal cross-section is characterized by an angle Δω_GNP (Fig. 1c). To measure Δω_GNP, we attached GNPs to the surface of beating axonemes, which were then imaged at 5,000 frames per second (f.p.s.) (Fig. 4 and Supplementary Movie 3). Using comparative waveform analysis and hydrodynamic modelling, we verified that the GNP attachment does not influence the axonemal beat (Supplementary Figs. 7 and 8). We determined the projected distance d_C between the GNP and the axoneme centre line in 2D images with nanometre precision (Fig. 4a, Methods and Supplementary Fig. 9). We find that d_C oscillates at the frequency of the axonemal beat with peak-to-peak amplitudes ranging from 13 nm to 124 nm depending on the azimuthal and longitudinal positions of the GNPs on the axoneme (see Fig. 4b for a typical example). Using simulated data, we confirmed that the observed peak-to-peak amplitudes are larger than the potential curvature-dependent tracking errors with expected magnitude of 8 nm (or less) for curvature values of 0.8 rad µm⁻¹ (or less) (Supplementary Fig. 10). Furthermore, we confirmed that d_C errors due to hydrodynamic drag forces acting on the GNP or due to possible defocusing as result of z movement of beating axonemes are both negligible (Supplementary Figs. 11–14). To reduce noise, we computed the beat-cycle average d_C(ϕ) as a function of the beat-cycle phase ϕ for each GNP individually, similar to the computation of the average 3D waveform (Fig. 4c, red line). The change of d_C during the beat cycle is indicative of a rotation of the local axonemal cross-section. From the d_C(ϕ) profile, we calculate the cross-section rotation angle ω_GNP(ϕ) with respect to the laboratory frame, using the known radii of the axoneme (r_axoneme = 100 nm)^29,30 and the GNP (r_GNP = 25 nm) (Fig. 4c, right axis). This defines the peak-to-peak rotation amplitude Δω_GNP of the oscillating rotation angle (Fig. 4d). Finally, we determine Δω_GNP(s) as a function of the arc length s by combining the results from all the GNPs (Fig. 4e). In the following, Δω_GNP(s) and Δω_3D(s) are used to compare the rotation of the axonemal cross-section with the rotation of the bending plane.

a, A high-speed dark-field microscopy image (exposure time 193 µs) of a reactivated axoneme with an attached GNP (scale bar 1 μm; right panel: zoom-in). We fit a Gaussian model describing the intensity profile of the GNP and the nearby axoneme (Methods and Supplementary Fig. 9) to precisely measure the centre line (blue) and the GNP position (red cross) to calculate their projected distance d_C in the 2D image (pixel size 73 nm). b, An exemplary plot of the time-dependent distance to the centre line d_C (Methods and Supplementary Figs. 15 and 16). c, The distance to the centre line d_C as a function of the beat-cycle phase ϕ (grey dots show pooled data from 67 beat cycles, while the red line shows the phase average), from which the rotation angle ω_GNP of the axonemal cross-section (in the laboratory frame) and its peak-to-peak amplitude Δω_GNP were computed (for more examples, see Supplementary Fig. 17). d, A visualization of Δω_GNP. e, The peak-to-peak amplitude Δω_GNP for axonemes, each with one or two GNPs attached at different arc-length positions (grey-filled circles) and averages after binning the arc-length positions (black-filled circles with whiskers, mean ± s.e.m., dashed lines indicate bin boundaries). The number of biological replicates in each bin was n = 9, 6, 3 or 2. Error statistics for bins with fewer then two data points are not displayed.

Source data

Full size image

Indication for twist–torsion coupling

If twist and torsion are strictly coupled, the measured torsion waves would imply that there are equal twist waves and that the axonemal cross-section and the bending plane would rotate together (Fig. 1c). We quantify the local bending plane rotation angle ω_3D(ϕ, s) from the average 3D waveform and compute the rotation amplitudes Δω_3D(s) (Fig. 5a, red curve; see also Supplementary Fig. 6a for a map of ω_3D(ϕ, s)). Because we measure Δω_3D(s) within the region of trust, it provides a lower bound for the rotation amplitude of the bending plane within the beat cycle. In the scenario of twist–torsion coupling, Δω_3D(s) should be equal to Δω_GNP(s), which we also evaluate in the region of trust (Fig. 5a, dashed blue curve). In the opposite scenario of twist-free torsion, there would be no twist and, hence, no axonemal cross-section rotation due to twist (Fig. 1d). However, an attached GNP could still show a rotation with peak-to-peak amplitude Δω_{no twist}(s) relative to the laboratory frame, due to possible axoneme rolling. To estimate ({Delta omega }_{{rm{no; twist}}}left(sright)), we first compute the hypothetical material frame of the axoneme, assuming that the axoneme has zero twist. We compute this twist-free frame by mathematically ‘untwisting’ the Frenet–Serret frame; we compute a new frame that rotates along the axoneme length at a rate that exactly cancels the torsion of the Frenet–Serret frame (for details, see Supplementary Fig. 18 and Supplementary Information). For every phase of the beat cycle, this twist-free frame is only defined up to a constant rotation around the axonemal centre line, similar to the motion of a universal joint. A hydrodynamic argument is needed to determine this constant rotation. According to the Helmholtz minimization principle³¹, any body immersed in a viscous fluid moves in such a way that the total rate of hydrodynamic dissipation is minimal. Hence, also the contribution due to a constant rotation of the twist-free frame must be minimal. This allows us to compute this rotation for every phase of the beat cycle (for details, see Supplementary Fig. 18b–d). The rotation of the twist-free frame directly provides the rotation of the GNP in the twist-free scenario, which we characterize by the peak-to-peak amplitude ({Delta omega }_{{rm{no}}, {rm{twist}}}left(sright)) (Fig. 5a, grey curve; see also Supplementary Fig. 18).

a, The rotation amplitudes Δω_GNP (light-blue circles), Δω_3D (rosé circles) and Δω_{no twist} (grey circles), together with the corresponding means ± s.e.m. in consecutive arc-length bins of length 0.25L (blue, red and black). The number of biological replicates for Δω_GNP in each bin was n = 9, 6, 3 and 2. Error statistics for bins with fewer then two data points are not displayed. b, A scatterplot of the differences between Δω_GNP(s) and Δω_3D(s) (measured in the region of trust) and between Δω_GNP(s) and (Delta {omega }_{{rm{no; twist}}}(s)) (measured in the region of trust and in the entire beat cycle). These differences are approximately normally distributed (Kolmogorov–Smirnoff test, with (Delta {omega }_{{rm{GNP}}}left(sright)-Delta {omega }_{3{rm{D}}}(s)): P = 0.67 and (Delta {omega }_{{rm{GNP}}}left(sright)-Delta {omega }_{{rm{no; twist}}}(s)): P = 0.44 in the region of trust and P = 0.70 in the full beat cycle). The average difference between (Delta {omega }_{{rm{GNP}}}(s)) and (Delta {omega }_{3{rm{D}}}(s)) is 5.6 ± 3.5° (mean ± s.e.m., n = 20), with a displayed 95% confidence interval of [−1.6 to 12.9], and not significantly different from zero (Student’s t-test, two-tailed, P = 0.12). The average difference between (Delta {omega }_{{rm{GNP}}}(s)) and (Delta {omega }_{{rm{no; twist}}}(s)) in the full beat cycle is 21.9 ± 3.2° (mean ± s.e.m., n = 20), with a displayed 95% confidence interval of [15.3 to 28.6], significantly different from zero (Student’s t-test, two-tailed, P < 0.001). The average difference between (Delta {omega }_{{rm{GNP}}}(s)) and (Delta {omega }_{{rm{no; twist}}}(s)) in the region of trust is 14.9 ± 2.9° (mean ± s.e.m., n = 20), with a displayed 95% confidence interval of [8.9 to 20.8], significantly different from zero (Student’s t-test, two-tailed, P < 0.001). Note that the rotation amplitude differences are robust against variability in the axoneme diameter (Supplementary Fig. 19). c, The rotation angles ({omega }_{{rm{GNP}}}(phi ,s)) (blue line), ({omega }_{3{rm{D}}}left(phi ,sright)) (red line, shown for the region of trust) and ({omega }_{{rm{no; twist}}}left(phi ,sright)) (black line) for three example arc-length positions (s = 0.28L, 0.34L and 0.59L). Since the sign and the offset between the curves are unknown, these were determined by minimizing squared residuals. d, A scatterplot of the unsigned Pearson correlation coefficient |r| calculated between ({omega }_{{rm{GNP}}}(phi ,s)) and ({omega }_{3{rm{D}}}(phi ,s)). The 16 filled circles indicate significant correlations (P < 0.05). The four open circles indicate non-significant correlations (P > 0.05).

Source data

Full size image

When we compare the above-described rotation amplitudes, we find that ({Delta omega }_{{rm{GNP}}}left(sright)) matches ({Delta omega }_{3{rm{D}}}left(sright)) but not ({Delta omega }_{{rm{no; twist}}}left(sright)) (Fig. 5a). Specifically, the mean values of ({Delta omega }_{{rm{GNP}}}left(sright)) and ({Delta omega }_{3{rm{D}}}left(sright)), calculated for each arc-length bin, are not significantly different (Student’s t-test, two-tailed, unpaired, α = 0.05). Likewise, the average difference ({Delta omega }_{{rm{GNP}}}(s)-{Delta omega }_{3{rm{D}}}(s)) within the region of trust is not significantly different from zero (Fig. 5b). In contrast, the rotation amplitude calculated for the twist-free torsion scenario ({Delta omega }_{{rm{no; twist}}}) is consistently lower than the experimental measurement of ({Delta omega }_{{rm{GNP}}}) (Fig. 5a). Specifically, the mean values of ({Delta omega }_{{rm{GNP}}}left(sright)) and ({Delta omega }_{{rm{no; twist}}}left(sright)), calculated for each arc-length bin, are significantly different (Student’s t-test, two-tailed, unpaired, α = 0.05), except for the most distal quartile, where n = 2 for ({Delta omega }_{{rm{GNP}}}). Additionally, we quantify the difference ({Delta omega }_{{rm{GNP}}}({rm{s}})-{Delta omega }_{{rm{no; twist}}}({rm{s}})) within the region of trust, which is significantly different from zero (Fig. 5b). When we consider ({Delta omega }_{{rm{GNP}}}) and ({Delta omega }_{{rm{no; twist}}}) not only in the region of trust but in the entire beat cycle, we find that the mean values of ({Delta omega }_{{rm{GNP}}}left(sright)) and ({Delta omega }_{{rm{no; twist}}}left(sright)) in each quartile are significantly different and that this difference deviates from zero even more strongly (Fig. 5b). Thus, our data are not consistent with the hypothesis that the 3D waveform is generated without twist but consistent with twist–torsion coupling.

Time correlation in cross-section rotation and bending plane

If twist and torsion are coupled, not only the peak-to-peak amplitudes of the GNP rotation and the bending plane rotation but also their dynamics during the beat cycle should agree. The direct comparison between ({omega }_{{rm{GNP}}}left(phi ,,sright)) and ({omega }_{3{rm{D}}}(phi ,,s)) within the region of trust (Fig. 5c) shows that their dynamics agree. We computed the correlation between ({omega }_{{rm{GNP}}}(phi ,s)) and ({omega }_{3{rm{D}}}(phi ,s)) (unsigned Pearson coefficient r) and find a significant correlation for 16 out of 20 GNP positions (P < 0.05, mean correlation 0.82; Fig. 5d). In the basal region (below s = 0.2L, corresponding to the most proximal 2.5 µm of the axoneme), we find three GNPs (out of nine GNPs measured at similar location) with non-significant correlations. This we attribute to measurement uncertainties in this region where the rotation amplitudes are low. Taken together, the temporal correlation of the axonemal cross-section rotation and the bending plane rotation further supports the hypothesis of twist–torsion coupling.

Torsion waves in beating axonemes shape their 3D waveform

We measured the 3D waveform of isolated Chlamydomonas axonemes with unprecedented spatio-temporal resolution (32 averaged shapes in a beat cycle, ({bar{sigma }}_{{xy}}) = 0.19 nm and ({bar{sigma }}_{z}) = 2.20 nm) using a combination of high-speed 3D shape reconstruction and beat-cycle averaging. This waveform is made available with this publication (Source Data). Our approach enables the measurement of high-resolution torsion profiles and provides rigorous error bounds. Within a beat cycle, the torsion oscillated with an amplitude of 20.0° µm⁻¹ (with error below 4.8° µm⁻¹), consistent with a hetero-chiral torsion wave (combining left- and right-handedness) propagating from the basal end of the axoneme to its distal tip. Our results are in agreement with earlier studies on reactivated Chlamydomonas axonemes²¹. Torsion generates non-planar shapes. In isolated axonemes, we found that the non-planarity is maximal during the recovery stroke (ϕ = π) with strong negative torsion in the bend region. This observation is consistent with measurements of the non-planarity in intact Chlamydomonas cilia³².

It is known that Chlamydomonas cells rotate during swimming, which is important for their phototaxis³². The non-planarity of the 3D waveform measured here may contribute to this cell rotation. However, the observation that Chlamydomonas cells rotate predominantly during the power stroke of their cilia³³ (when the beat is mostly planar) and not during the recovery stroke (when we observe strongly non-planar shapes) suggests another relevant contribution, possibility related to twist near the axonemal base.

Twist and torsion are coupled in beating axonemes

By comparing the cross-section rotation (measured by GNPs attached to the surface of axonemes) with the rotation of the local bending plane (measured from the 3D shapes), we provide evidence for the coupling of torsion and twist. We find that the bending plane and axonemal cross-section rotate together. This implies a highly anisotropic bending stiffness of the axoneme, which strongly supports the rigid-bridge hypothesis, but not the opposite hypothesis of twist-free beating. We believe that structural heterogeneities in the Chlamydomonas axoneme together generate this anisotropic bending stiffness. Among those structures is the 1–2 bridge¹⁴, which crosslinks DMT1 and DMT2 in the proximal quarter of the axoneme and the IDLs^15,34. IDL3, in particular, is present along the entire length of DMT1 in Chlamydomonas and is thought to rigidly cross-link DMT1 and DMT2. Additionally, outer arm dynein is missing between this pair of DMTs¹⁴. These structural inhomogeneities may restrict sliding and contribute to setting the bending plane. This is consistent with measurements of inter-doublet sliding in reactivated Chlamydomonas axonemes using cryo-electron tomography, where the sliding amplitude was low between DMT1 and DMT2 but high normal to the plane of the bridge³⁵.

Similar and additional inter-doublet links are found in other organisms. For example, in sea urchin, Trypanosoma, paddle cilia, gill cilia of a lamellibranch mollusc, mouse and human sperm flagella, a bridge structure connects DMT5 and DMT6^15,36,37,38. Therefore, we expect that structural heterogeneities, which generate an anisotropic bending stiffness, are a conserved feature of the axonemal structure.

A recent theoretical study investigated numerically how an anisotropic bending stiffness and non-zero twist rigidity impacts simulated cilia beat patterns³⁹. Already for moderately anisotropic bending stiffnesses (with a ratio of bending stiffnesses equal to 2.5), beat patterns became more planar, with bending predominantly in the soft direction and non-zero twist. In contrast, if the imposed bending stiffness was isotropic, symmetric helical beat patterns emerged that had almost zero twist. These two numerical cases correspond to the limit cases considered here: the limit case of highly anisotropic bending stiffness and low twist rigidity corresponds to the scenario of perfect twist–torsion coupling, while the limit case of isotropic bending stiffness and large twist rigidity corresponds to the twist-free scenario. We anticipate that the numerical findings from Rallabandini et al.³⁹ are generic and should thus qualitatively hold also in other proposed models of motor control in the axoneme. Our data show that reactivated Chlamydomonas axonemes beat consistently with the first limit case of anisotropic bending stiffness but not the second limit case of isotropic bending stiffness. In the presence of such an anisotropic bending stiffness, the observed hetero-chiral torsion waves are indicative of equal hetero-chiral twist waves. Our analysis provides the missing link between torsion (characterizing centre line shapes) and twist (characterizing mechanical deformations).

Relation to past work

The structural similarities among motile axonemes across species (for example, the bridge structure) strongly suggest that twist–torsion coupling may be a general feature^40,41. To measure the 3D waveforms of beating cilia and flagella, previous studies employed stereographic imaging⁴², digital in-line holography^16,43, multi-focal 2D dark-field microscopy¹⁷, defocused 2D bright- and dark-field microscopy^18,27 as well as multi-focal phase contrast microscopy²¹. However, these studies were limited in spatio-temporal accuracy, which made it difficult to quantify torsion reliably. Nonetheless, the consistent observation of non-planar shapes of Chlamydomonas cilia³², Paramecium cilia⁴², sperm flagella^16,17,18,44 and Malaria and Trypanosoma parasites^43,45 indicates the existence of torsion. Assuming that twist–torsion coupling is a general feature of beating axonemes, previous torsion measurements can be interpreted as indirect twist measurements.

Our observation of twist in beating axonemes is in agreement with twist measurements in fixed axonemes from Paramecium⁴⁶, Trypanosoma⁴⁷ and Ciona sperm⁴⁸. The observation of the lateral movement of mitochondria during the beat of surface-attached quail sperm²⁵ is consistent with twist waves that we find in Chlamydomonas axonemes.

Twist generation in the axoneme

Our finding of hetero-chiral twist waves raises the question of how twist is generated. The axoneme has an inherently chiral architecture where dynein motors attached to each DMT exert active forces on the clockwise neighbouring DMT (when viewed from the base, Fig. 1a). Active dyneins slide adjacent DMTs. Sliding is restricted at the base. Bending and twist are then generated according to the following geometric concepts: (1) If there is a difference in DMT sliding on opposite sides of the axoneme cross-section, the axoneme bends⁴⁹; (2) If there is net sliding, defined as the sum of signed sliding displacements around the circumference of the axoneme, the axoneme twists⁵⁰. We illustrate these geometric concepts in three examples: (1) If dyneins induce active sliding on only one side of the axoneme and DMTs on the opposite side are free to slide in the opposite direction, the net sliding is zero and there is bending but no twist; (2) If dyneins are active on both sides of the axoneme and induce the same amount of sliding, the net sliding is non-zero and there is twist but no bending; (3) If dyneins induce an unequal amount of sliding on opposite sides of the axoneme, there will be both bending and twist. Due to the chiral arrangement of the dyneins that walk towards the microtubule minus ends (which is towards the base), dynein-generated sliding always induces dextral twist.

How sinistral twist is generated is unclear. Potentially, the axoneme might be already sinistrally twisted in the absence of dynein forces. Interestingly, structural sinistral twist was reported for non-motile 9 + 0 axonemes of human islet cilia⁵¹. Sinistral structural components were also reported for motile 9 + 2 axonemes, such as the central apparatus^24,52, or the heads of the radial spokes⁵³. On the other hand, axonemal twist could also be generated by the side-stepping of dyneins²⁰. In gliding assays, it was observed that dyneins can rotate microtubules clockwise when viewed from the minus end^54,55. We argue that such a rotation generated by side-stepping would translate into sinistral twist in the axoneme. Taken together, hetero-chiral twist waves could emerge if the axoneme is twisted sinistrally by default and gets periodically unwound and even dextrally over-twisted by dynein-generated sliding.

Mechanisms of motor control

The precise mechanism of dynein coordination driving the axonemal beat remains a matter of debate^{56,57,58,59,60,61,62,63,64}. Almost all theoretical models assume that motor activity is regulated by mechanical deformations such as curvature, sliding or DMT spacing. A recent theory suggested that twist contributes to regulate dynein activity²⁰. In short, in an axoneme that is both bent and twisted, DMT spacing increases on one and decreases on the opposing side of the axoneme, which could regulate dynein activity. This and other theories prompt precise measurements of time-resolved mechanical deformations of the axoneme, such as the twist investigated here. Specifically, by providing evidence for twist in bent axonemes, our study presents experimental support for a twist-assisted curvature control model as proposed in ref. ²⁰. Our quantitative measurements of axonemal twist waves with peak-to-peak amplitudes of 20° µm⁻¹ imply an accumulated twist angle of approximately >15° over distances of 2 µm. This is in agreement with the theoretical prediction by ref. ²⁰, where an accumulated twist angle of 14° was sufficient to generate transversal forces strong enough to change the DMT spacing. In conclusion, our data inform models of motor coordination in beating axonemes.

Summary and outlook

We provide strong evidence for twist waves in beating Chlamydomonas axonemes. This twist generates non-planar 3D waveforms, which can contribute to fluid pumping, cell rotation and helical swimming, required for the physiological function of cilia and flagella. Twist is a dynamic mechanical deformation generated by active motor forces. While motors twist the axoneme, this deformation may regulate their activity through a mechanical feedback loop.

Experiments

Isolation and reactivation of Chlamydomonas axonemes with GNPs

Axonemes from Chlamydomonas reinhardtii wild-type cells (cc-125 mt+) were isolated and reactivated following the experimental procedures described in ref. ⁶⁵. All reagents were purchased from Sigma Aldrich, unless stated otherwise. In brief, cells were grown in TAP medium under constant illumination by light-emitting diode light pads (light therapy lamp HST-001, 10,000 lux, 10 W, 4,500 K) and air bubbling at 22–24 °C for 3–4 days to final density of (3–7) × 10⁶ cells ml^–1. Cilia were isolated by the dibucaine procedure, separated from the cell bodies by centrifugation (2,400g, 25% sucrose cushion) and demembranated in HMDEK (30 mM HEPES-KOH, 5 mM MgSO₄, 1 mM DTT, 1 mM EGTA, 50 mM potassium acetate, pH 7.4) augmented with 1% (v/v) IGEPAL and 0.2 mM Pefabloc SC. The membrane-free axonemes were resuspended in HMDEKP (HMDEK with 1% (w/v) polyethylene glycol, 20 kDa) with 30% sucrose, 0.2 mM Pefabloc added, and stored at −80°C. Prior to reactivation, axonemes were thawed at room temperature, then kept on ice for at most 3 h. The reactivation was performed in flow chambers (with a depth of 100 µm) built from easy-cleaned glass and parafilm⁶⁵. For the experiments with GNPs, 1 µl of thawed axonemes was mixed with 1 µl of GNP solution (diameter 50 nm, 3.5 × 10¹⁰ ml⁻¹) and incubated for 10 min on ice. Axonemes with or without GNPs were diluted in HMDEKP reactivation buffer containing 1 mM adenosine triphosphate (ATP) and an ATP-regeneration system (1 unit per millilitre creatine kinase, 5 mM creatine phosphate) used to maintain the ATP concentration. The axoneme dilution was infused into a glass chamber, which was blocked with casein solution (from bovine milk, 2 mg ml⁻¹) for 10 min. Before imaging, the flow chamber was sealed with twinsil-speed (picodent) to avoid evaporation. The sample was equilibrated on the microscope for 5 min before data were collected for a maximum time of 60 min on the side of the glass chamber. Note that, when viewed in an inverted microscope, axonemes on the bottom and top surface of the observation chamber show opposite sense of rotation, although their swimming paths and beats have the same chirality⁶⁶ (Supplementary Fig. 1). We further note that we cannot exclude a bias in the reactivation of axonemes towards axonemes originating from cis-cilia, similar to previous studies in cell models⁶⁷.

GNP attachment to the axoneme

The GNPs attach non-specifically and at random positions on the axoneme. On the basis of previous work, we expect that the gold surface will bind by a sulfur–gold reaction to cysteines of the tubulin dimers, with rupture forces in the nanonewton range^68,69. With respect to this tight binding, the forces that act on the GNP during axonemal beating are low (<100 fN; Supplementary Table 1). Further, we found no significant difference between waveforms of axonemes with and without attached GNPs (Supplementary Fig. 7).

Imaging of axonemes

Reactivated axonemes were imaged by dark-field microscopy set up on an inverted Nikon ECLIPSE Ti2-E microscope, using a Nikon 100× iris oil (0.5–1.2 numerical aperture (NA)) lens in combination with a 1.5× or 1.0× tube lens and a Nikon 1.3 NA oil condenser. Images were recorded with a pco.dimax CS4 high-speed camera. The z scan images for 3D calibration were recorded with a pco.edge 4.2 camera. In both cases, the sample was illuminated using a Sola SE2 light engine (Lumencor) combined with a 496 LP filter (Semrock, Brightline). Movies of reactivated axonemes with GNP were recorded with the 1.5× tube lens and an NA setting of 1.1–1.2 on the objective at a frame rate of 5,000 f.p.s. In total, we recorded more than 150,000 image frames from a total 20 GNPs (5,000–10,000 images each) attached to 19 beating axonemes (one with two GNPs), corresponding to a total of more than 1,750 beat cycles. For the reconstruction of the 3D waveform, we recorded movies at 1,000 f.p.s. using a 1.0× tube lens and an objective NA of 1.0. For 3D reconstruction, we first focused below the axoneme, so that the axoneme never entered the focal plane. We then recorded a total of 53,000 defocused images from a total 17 axonemes (3,000–5,000 each), corresponding to a total of more than 3,750 beat cycles.

Data analysis

Calibration of defocused dark-field microscopy and 3D shape reconstruction

We reconstructed the 3D shape of the axonemal centre line in single defocused images of beating axonemes from the x–y position measured by 2D filament tracking and the z position derived from the FWHM of intensity distributions measured normal to the long axis of the axoneme. All analysis was done with MATLAB (version R2020a) unless stated otherwise. Specifically, we subtracted the static background (the median intensity in each pixel calculated for the entire movie) using FIJI. We then determined the 2D centre line using the filament tracking software FIESTA²⁸. We interpolated the 2D centre line using a smoothing spline fit (MATLAB) to gain equidistant points with a spacing of 110 nm along the centre line. Subsequently, we determined the intensity distributions along line scans normal to the local tangent of the 2D centre line at each control point of the 2D axoneme centre line and measured the local FWHM by Gaussian fitting. We corrected the observed FWHM on the basis of the transverse distance travelled to eliminate contributions from motion blur (Supplementary Fig. 2). To relate the FWHM to an absolute z position with respect to the focal plane, we performed the following calibration: We non-specifically immobilized axonemes (without ATP) on a dichlorodimethylsilane-coated glass coverslip. Subsequently, we recorded a z scan with a range of ±10 µm relative to the z position at which the axoneme was in focus. To obtain a calibration curve, we measured the FWHM as a function of the distance to the focal plane (the z position where the FWHM is minimal; the calibration curve is shown in Fig. 2e, smoothed using a spline fit with smoothing parameter of 0.9). To account for the refractive index mismatch between the immersion oil and water (reactivation buffer), we corrected our calibrated z positions by a factor of 1.137, which is the ratio between the two refractive indexes (n_oil = 1.515, n_water = 1.333). We estimated the localization uncertainty of 3D shapes obtained in single images as σ_xy = 7.3 nm and σ_z = 35.6 nm (Supplementary Fig. 3).

Beat-cycle phase and beat-cycle averaging

To determine the beat-cycle phase ϕ of a given shape, we fitted the arc-length profile of the 2D tangent angle (time average subtracted) with a sinusoidal function (Supplementary Fig. 4a–d). By this definition, the recovery stroke included shapes with the beat-cycle phases [0.5π ≤ ϕ < 1.5π]. The remaining shapes belonged to the power stroke (Supplementary Fig. 4d). For details on the determination of the beat-cycle phase, see Supplementary Fig. 4a–d. To calculate a high-resolution 3D waveform of reactivated Chlamydomonas axonemes, we used beat-cycle averaging. To do this, we subdivide the beat cycle into 32 bins (with equal width of ∆ϕ = π/16) and group shapes of similar beat-cycle phase. For each bin, we calculated the average shape with 30 points along the arc length (∆s_2D = L/30). To do this we: (1) averaged the profiles of the 2D curvature along the arc length ({kappa }_{2{rm{D}}}({s}_{2{rm{D}}})) and the profiles of the z position along the arc length (zleft({{rm{s}}}_{2{rm{D}}}right)) of all the recorded axoneme shapes, and (2) calculated the average 2D position (xleft({s}_{2{rm{D}}}right)) and (yleft({s}_{2{rm{D}}}right)) from the average curvature ({kappa }_{2{rm{D}}}({s}_{rm 2D})) by integration. We obtained the averaged x, y, z positions along the 2D centre line for each beat-cycle bin, which together compose the average 3D waveform. For details on the averaging method, see Supplementary Fig. 4e–j. We estimate the localization uncertainty of the average 3D shapes obtained through beat-cycle averaging as σ_xy = 0.19 nm and σ_z = 2.20 nm (Supplementary Fig. 4k).

Measurement of the distance between GNPs and the axoneme centre line

We precisely measured the distance between the GNP and the axoneme 2D centre line d_C in each recorded image. To do this, we used FIESTA²⁸ to roughly measure the position of the axoneme-bound GNP in the image. Using this position, we selected a region of interest of 25 × 25 pixels, which contained the GNP in the centre. Within this box, we fitted a 2D intensity model, which was the sum of a symmetric 2D Gaussian function and a Gaussian wall function. The centre of the symmetric 2D Gaussian described the centre position of the GNP, while the centre line of the Gaussian wall described the axoneme centre line. We used both to calculate the normal distance between the GNP and the axoneme centre line, which we called the distance to the centre line d_C. We subtracted an orientation-dependent systematic error from the measured d_C (as detailed in Supplementary Figs. 15 and 16). To calculate the average d_C as a function of the beat-cycle phase for a single GNP, we performed beat-cycle averaging as described above. We used the average d_C(ϕ) to calculate the angular position of the GNP in the cross-section as observed from the laboratory coordinate system ({omega }_{{rm{GNP}}}). Note that, since d_C is only a two-dimensional projection, we cannot measure whether the GNP is attached above or below the axoneme, which would change the sign of ({omega }_{{rm{GNP}}}). We quantified the local cross-section rotation by the peak-to-peak amplitude (Delta {omega }_{{rm{GNP}}}). We determined the approximate point of GNP attachment along the arc length of the axoneme via a segmented line measurement tool (FIJI) and normalized this distance by the axoneme length.

Calculation of the 3D curvature, torsion and ω
_3D

A space curve r(s) is parameterized by the arc length s with three unit vectors of the Frenet–Serret frame: the tangent vector t(s), the normal vector n(s) and the binormal vector b(s) at each arc-length position s. These three vectors are mutually orthogonal with b(s) = t(s) × n(s). If the curvature κ(s) is non-zero at r(s), and thus the local bending plane is well defined, b(s) is perpendicular to this plane, which is spanned by t(s) and n(s). Here, the curvature κ is defined as the rate at which the tangent vector t(s) rotates (within the bending plane) as a function of s. The torsion τ(s) is defined as the rate at which the bending plane rotates. Mathematically

The measured 3D shapes represent discretized space curves for which we calculate approximations for the unit vectors of the Frenet–Serret frame t_s, n_s, b_s. We approximate (1) the local tangent vector at a given point with the mean of the two adjacent secant vectors of the discretized space curve, (2) the local binormal vector as the cross-product of the two adjacent tangent vectors and (3) the local normal vector as the cross-product of the binormal vector and the tangent vector. We require that the z component of the binormal vector is always positive. Thereby, we can define a unique sign of the curvature for the entire 3D shape. Using the approximation of the Frenet–Serret frame, we calculate the torsion and the curvature as the arc-length derivative of the signed rotation angle of the local bending plane ω_τ and the signed in-plane rotation angle ψ, respectively. Additionally, we obtain the local bending plane orientation with respect to the laboratory x–y–z coordinate system ω_3D as the angle between the normal vector and the x–y plane. We present these 3D shape parameters as functions of the 3D arc length s and the beat-cycle phase ϕ as (kappa (phi ,s)), (tau (phi ,s)) and ({omega }_{3{rm{D}}}left(phi ,sright)) (Fig. 3d,e and Supplementary Fig. 6a).

The measurement uncertainty and the choice of the region of trust for torsion

To estimate the uncertainty of the torsion measurement, we used a bootstrapping approach. Our experimental dataset contained 52,966 shapes (from 17 axonemes). From those shapes, we drew 1,000 sets of 52,966 shapes with replacement. For each set, we calculated the average 3D shape as well as the torsion (tau (phi ,s)). The standard deviation over all 1,000 computed torsion maps was used as an estimator of the measurement uncertainty of the torsion. To determine the region of trust of the torsion measurement, we used the argument that the normal (and binormal) vector of the Frenet–Serret frame is ill defined in regions of low curvature. Since torsion is calculated using those vectors, the variability of the torsion increases with decreasing absolute curvature (Supplementary Fig. 5d–f). On the basis of this dependence, we defined a conservative region of trust in which the variability is approximately constant, which is the case for values of local curvature of |κ| > 0.4 rad µm⁻¹ and above. The same region of trust was used for the rotation angle ω_3D of the local bending plane relative to the laboratory frame.

Sign conventions

We consider the x–y plane viewed along the negative z vector of the laboratory frame, that is, from inside the chamber (Fig. 3a).

We defined the sign of the curvature κ such that the time-averaged axonemal shape projected on the x–y plane (static waveform component⁷⁰) has positive curvature, corresponding to bending in the clockwise direction (when traversed from base to tip). Positive torsion corresponds to dextral (right) handedness (as shown in Fig. 1b) and negative torsion to sinistral (left) handedness.

The sign of ({omega }_{3{rm{D}}}) and ({omega }_{{rm{no; twist}}}) was chosen to be opposite to the sign of the z component of the normal vector at the corresponding arc-length position (Supplementary Fig. 6a). The sign of d_C is defined with respect to the proximal–distal polarity of the axoneme (positive/negative when the GNP is on the left/right side of the 2D centre line, when it is traversed from base to tip) (Supplementary Fig. 9). The sign of the measured ({omega }_{{rm{GNP}}}) is arbitrary for individual GNPs since we do not know on which side of the axoneme the GNP was bound. To compare the temporal dynamics of ({omega }_{{rm{GNP}}}(phi ,s)) and ({omega }_{3{rm{D}}}left(phi ,sright)) in Fig. 5d, we use the unsigned Pearson correlation coefficient, which is independent of whether the GNP was bound to the upper or lower side of the axoneme.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.