Investigation of multiple-dynein transport of melanosomes by non-invasive force measurement using fluctuation unit χ

Pigment organelles known as melanosomes disperse or aggregate in a melanophore in response to hormones. These movements are mediated by the microtubule motors kinesin-2 and cytoplasmic dynein. However, the force generation mechanism of dynein, unlike that of kinesin, is not well understood. In this study, to address this issue, we investigated the dynein-mediated aggregation of melanosomes in zebrafish melanophores. We applied the fluctuation theorem of non-equilibrium statistical mechanics to estimate forces acting on melanosomes during transport by dynein, given that the energy of a system is related to its fluctuation. Our results demonstrate that multiple force-producing units cooperatively transport a single melanosome. Since the force is generated by dynein, this suggests that multiple dyneins carry a single melanosome. Cooperative transport has been reported for other organelles; thus, multiple-motor transport may be a universal mechanism for moving organelles within the cell.

triphosphatase that relieves the autoinhibition of the motor, which is critical for avoiding unnecessary consumption of ATP when the motor is not bound to an SVP 23 -had fewer FPUs than their wild-type counterparts 12 .
In this study, we investigated χ for the transport of melanosomes, organelles filled with the melanin pigment, in zebrafish melanophores. Melanosomes are transported by microtubule motors, kinesin-2 and cytoplasmic dynein, and the actin motor myosin-V 24 . They disperse or aggregate in response to hormones in melanophores. The primary physiological purpose for movement of melanosomes in animals is colour change. Particularly, we focused on the aggregation process of melanosomes transported by dynein. There are several advantages of measurement of χ in this transport system. Firstly, because melanin pigment is black and easily observed by bright-field microscopy, the recording rate can be increased up to 800 frames per second (fps). Such a high-speed recording enables accurate measurement of fluctuation in the position of melanosomes. Secondly, we can decrease the number of dynein motors using the inhibitor ciliobrevin, which was recently identified 25 and investigated for many biological phenomena related to dynein [25][26][27][28][29] , thereby allowing observation of the behaviour of fluctuation unit χ in response to the decrease in motors. Thirdly, the force generation mechanism of dynein is still controversial and therefore worth studying. The stall force value for a single dynein molecule was reported to be 1 pN 9,10,30 , though a few studies have reported values of 5-10 pN [31][32][33] . Fourthly, the behaviour of χ remains to be evaluated in non-neuronal cells. When considered together, melanosome transport is a suitable system to explore χ.
In our experiments, the quantal behaviour of χ was observed in melanosome transport by dynein. We concluded that several FPUs carry a single melanosome similarly to cargo transport in neurons [11][12][13] . Further, we found that the number of FPUs measured by χ was dramatically decreased by the addition of ciliobrevin, as expected. We anticipate that the non-invasive force measurement using the fluctuation unit (χ) will be applied to a wide range of cargo transport in eukaryotic cells to elucidate the physical mechanisms in which some human diseases are deeply rooted.

Results
observation of melanosome transport. Melanophores removed from zebrafish scales by enzymatic treatment were cultured in a glass-bottom dish. After 1 day of culture, the motion of melanosomes was observed by bright-field microscopy. Following addition of hormone (epinephrine) to the dish, melanosomes aggregated to the centre of a melanophore. Because the minus ends of microtubules are attached to the centrosome, the motion towards the nucleus is caused by the action of dynein. A schematic of this transport is represented in Fig. 1.
Although they were crowded at the start of migration, individual melanosomes could be tracked after several tens of seconds, when nearly all melanosomes were aggregated. The direction of movement was set as the plus X direction (e.g. red line in Fig. 2a). Except for melanosomes close to the nucleus where cells had greater thickness, we considered that the motion in the Z direction was less than that in the X-Y plane since the images of melanosomes were focused during the runs. The centre position (X) of the melanosome was calculated from the recorded images, and the time course of X was obtained (Fig. 2b). Directional movement was defined as movement at a velocity greater than 100 nm/s, as in a previous study on pigment granule transport 5 . The fluctuation in melanosome position during directional transport, which was observed at a high recording rate of 800 fps (inset in Fig. 2b), was mainly due to thermal noise, stochastic stepping of motors accompanied by ATP hydrolysis, and collision of melanosomes with other organelles and cytoskeletons. Schematic illustration of melanosome transport in a melanophore. A melanosome (black circle) is transported by multiple dynein motors (grey) along a microtubule (pink) in a melanophore. During the run at a velocity v, the force (F) generated by motors is equal to the drag force (Γv), where Γ is the friction coefficient of the melanosome. Addition of hormone (epinephrine) caused melanosomes to aggregate in the centre of the cell (bottom). Individual melanosomes could be tracked under a microscope. The schematic was created by the authors.
www.nature.com/scientificreports www.nature.com/scientificreports/ Calculation of χ for melanin pigment transport. The constant velocity segment (red area in Fig. 2b) of the position (X) of a melanosome lasting 0.5-1 s was identified by visual inspection for calculation of the fluctuation unit χ (equation (6)). The effect of this selection on the calculation of χ was later evaluated by the bootstrapping method. Based on this segment, we calculated ΔX = X(t + Δt) − X(t) (e.g. inset in Fig. 2b) and its probability distribution (P(ΔX)) ( Fig. 2c for the case Δt = 37.5 ms), where the value of Δt ranged from 1.25 to 37.5 ms. P(ΔX) was well fitted by a Gaussian function (solid black curve in Fig. 2c). The Gaussian function (equation (5)) was used to fit P(ΔX) after evaluating the Brownian noise power spectrum density S(f) (equation (11)) of X for small f (Supplementary Fig. S1).
χ was calculated using equation (6) for each Δt (Fig. 2d); χ as a function of Δt converged to a constant value (χ*) after relaxation (∼20 ms), which is attributable to the microscopic environment around the melanosome as well as the enzymatic cycle time 13 . Note that from a physical point of view, the convergence implied that the macroscopic motion of X whose time scale is about 20 ms obeys the simple equation (3). Then, equation (7) was considered for the constant value χ*.
The error of χ* according to range selection of the constant velocity segment was estimated to be 15% based on the principle of bootstrapping. We first selected 10 different partial segments from the original constant velocity segments; χ was then calculated for each partial segment (thin curves in Fig. 2d). It should be noted that the length of partial segments was half that of the original segment.
Melanosome transport by multiple dynein motors. The procedures carried out in Fig. 2b-d for a single melanosome were repeated for 62 different melanosomes from two different melanophore preparations. The calculated values of χ for the 62 melanosomes are shown in Fig. 3a. These data were classified into four groups by applying affinity propagation (cluster) analysis (Fig. 3b). When we examined the q value-the sole parameter of this analysis-a cluster number of 4 was the most stable within the range 0 ≤ q ≤ 1 (Fig. 3c). The histogram of χ* (equation (10)) showed multiple peaks.  (6)). (a) Typical motion of a melanosome after adding hormone (epinephrine). The direction of motion (red) is the positive X direction. (b) Time course of the centre position (X) of the melanosome. When the recording rate was sufficiently high (800 fps), melanosomes were observed to fluctuate while directional motion was driven by dynein motors (inset in b). We focused on the constant velocity segment (red area). c Probability distribution (P(ΔX)) of ΔX in the case of Δt = 37.5 ms, where ΔX = X(t + Δt) − X(t) (inset in b) was calculated for the constant velocity segment. P(ΔX) was fitted with a Gaussian function (black curve). (d) χ (equation (6)) calculated using the fitting parameters (a,b) of the Gaussian function (equation (5)) plotted as a function of Δt (thick black curve). After relaxation time, χ reaches a constant value (χ* (equation (10)). The thick red curve represents χ after applying a smoothing filter (equation (9)). The thin curves (n = 10) represent χ calculated from different partial segments cut from the original constant velocity segments and were used to estimate the error of χ (15%).
The quantal behaviour of χ shown in Fig. 3a is understood as follows. The stall force acting on a single cargo was determined to be quantal using optical tweezers in living cells, reflecting the number of motors hauling the cargo 7,8 . From the schematic in Supplementary Fig. S2, it is understood that the drag force (F) becomes quantal like the stall force when the viscosity effect is large; this is thought to occur during the transport of melanosomes in melanophores. Thus, based on the assumption F ∝ χ* (equation (8)), the multiple peaks in the histogram of χ* (Fig. 3d) were considered to reflect multiple FPUs as well as the peaks of the stall force distributions. It was then verified whether the multiple peaks indeed represent multiple FPUs in the next paragraph. As circumstantial evidence, when the number (N i , where i = 1, 2, 3, and 4) of elements belonging to each FPU was investigated, N i /N 1 (inset in Fig. 3d), whose behaviour was similar to that observed in melanosome transport in Xenopus melanophores 5 , supported the meaning of multiple peaks of χ* and the underlying assumption (equation (8)).

Addition of ciliobrevin D.
To further verify whether χ can detect the number of FPUs carrying a melanosome, the number of active dyneins was intentionally decreased by adding the dynein inhibitor ciliobrevin D 25 . Figure 4a shows the time course of a melanosome in the presence of 40 μM ciliobrevin D. When we examined the time course-whose velocity was similar to that seen when the concentration of ciliobrevin ([ciliobrevin]) was 0 μM (Fig. 2b)-we found that the fluctuation (ΔX) increased (inset in Fig. 4a). This change in ΔX indicated that the fluctuation reflects the effect of the decrease in dynein caused by ciliobrevin. In Fig. 4b, the velocity values of constant velocity segments were compared between cases of [ciliobrevin] = 0 and 40 μM. As was reported for peroxisome transport in Drosophila S2 cells 25 , the high velocity portion was decreased in the presence of ciliobrevin.
The results of χ are shown for [ciliobrevin] = 10 μM (59 melanosomes from two melanophores), 20 μM (61 melanosomes from three melanophores), 30 μM (51 melanosomes from two melanophores), and 40 μM (48 melanosomes from two melanophores) (Fig. 5a-d). The number of FPUs tended to decrease with increasing [ciliobrevin]. Note that the clustering of χ was carried out by affinity propagation (Supplementary Fig. S3). The number and distribution of peaks corresponding to χ appeared to be correlated with the number of dynein motors hauling a melanosome.

Discussion
We observed the aggregation of melanosomes transported by dynein in zebrafish melanophores by bright-field microscopy at a high recording rate (800 fps). This high rate was needed to accurately measure fluctuations in the position of melanosomes. The motion of the melanosomes was quantified as fluctuation unit χ (equation (6)) [11][12][13] , which was determined on the basis of the fluctuation theorem of non-equilibrium statistical mechanics 17,21,22 . The values of χ calculated for 62 different melanosomes revealed their quantal behaviour (Fig. 3a). When equation (8) holds, the quantal behaviour of χ indicates the existence of multiple FPUs. Because force is generated by dynein, the existence of several FPUs indicates that a single melanosome is transported by multiple dynein motors. When the number of active dynein molecules was decreased by treatment with ciliobrevin D 25 , the number of FPUs estimated by χ also decreased. Indeed, χ was correlated with the number of motors, supporting equations (7) and (8). Thus, χ qualifies as an indicator of force generated by motors based on its behaviour in our experiments.
While χ was found to be related to the number of FPUs (Fig. 5), the components of one FPU are still unknown. A recent cryo-electron microscopy study revealed that two dimers of dynein were connected (i.e. four monomers) via adaptor proteins, suggesting collective force generation by dynein 34 . Although the force generated by a single dynein molecule is small 9,10,30 , there is a possibility that dynein molecules collectively exert a large force similar to that generated by the kinesin dimer. A future study will be necessary to elucidate the number of dynein monomers that constitute one FPU in the case of melanosome transport.
It was previously shown that velocity distributions of melanosomes in Xenopus melanophores had several peaks, reflecting multiple motors 5 . The velocity v becomes quantal when the force F of motors is quantal, and the friction coefficient Γ is almost the same for all melanosomes through the relation F = Γv. However, the velocity distributions measured in our experiments did not show multiple distinct peaks ( Supplementary Fig. S4) because Γ appeared to be distributed over a wide range. We speculated that Γ had a spatial dependence in our experiment because the velocity of melanosomes increased near the nucleus where the cell thickness increased. Because we used primary cultures of melanophores from different zebrafishes instead of a melanophore cell line as in a previous study 5 , and the melanophores in each experiment differed genetically, we observed variability in the viscosity of cytosol, which was influenced by melanophore structure. Additionally, in the previous study 5 , actin filaments were depolymerised with latrunculin B, which may have contributed to the constant value of Γ in the previous study.
In a previous study on melanosomes in Xenopus melanophores 5 , it was presumed that melanosome velocity doubled with the number of FPUs. To verify this assumption in our system, we investigated changes in velocity over a long time course (Supplementary Fig. S5). Based on the measurement of χ, we found that the velocity change was accompanied by a change in the number of FPUs, as expected from the previous study. On the other hand, for some melanosomes the velocity changed, while the number of FPUs did not ( Supplementary Fig. S5), possibly because of a change in friction coefficient Γ during the runs rather than a change in the number of FPUs; Γ may exhibit spatial dependence in our case since melanosomes near the cell centre often had a higher velocity.
The noise property of in vivo cargo transport is another important feature of the cellular environment. In a previous study on melanosome transport 35 , it was reported that the microenvironment around melanosomes is www.nature.com/scientificreports www.nature.com/scientificreports/ viscoelastic, based on the power spectrum densities of positions of moving melanosomes, which did not show Brownian noise behaviour. In this study, we focused only on the constant velocity segments (∼1 s) of melanosomes (Figs 2b and 3a) and found that the power spectrum densities for the segments showed Brownian noise behaviour for small frequencies (Supplementary Fig. S1). This Brownian noise property of X in constant velocity segments was also observed for mitochondrial transport 16 and endosome transport 13 but was found to not hold when the trajectory of a cargo included pauses and bidirectional motion 16 . As a result of the Gaussian noise property in equation (3) supported by the power spectrum densities for the segments, the quantitative analysis of force could be possible. Focusing on a constant velocity segment of in vivo cargo transport appears to be important for www.nature.com/scientificreports www.nature.com/scientificreports/ quantitative analysis of its dynamics as well as qualitative analysis of subdiffusive properties of long trajectories which include the characteristic complex motion of cargo in cells.
Determination of the proportionality constant (k B T eff ) between χ* and F is critical to estimate the in vivo force value produced by dynein motors as 1 FPU. In this study, we failed to determine the constant. In the case of axonal anterograde cargo transport by kinesin which was studied previously 13 , the proportionality constant was estimated by comparing the χ*-v relation with the known force-velocity relation of kinesin suggested by the in vitro stall force experiment 36 . The force-velocity relation could be fitted to the χ*-v relation of the neuronal cargo by a single value of k B T eff 13 . Unlike the case of anterograde transport, because the value of stall force produced by dynein motors estimated in the single-molecule experiments is controversial, with reports of 1 pN 9,10,30 to 5-10 pN [31][32][33] , the calibration of F from the χ*-v relation (Supplementary Fig. 6) was difficult in this study. Thus, it is essential to investigate equations (7) and (8) in combination with in vitro and in vivo stall force experiments using optical tweezers in future studies. The force-velocity relation of multiple dynein motors may lead to understanding of why χ* takes the discrete values of 0.5, 2, 4, and 6 for 1, 2, 3, and 4 FPUs (Fig. 4a), which are different from the discrete values of χ* in the case of axonal retrograde transport (0.1, 0.2, and 0.3) 13 . The different structures between the narrow axon and the body of eukaryotic cells, which may have wider spaces for easy cargo movement during transport, may affect the values of χ*. Moreover, T eff may have minor dependence on force, causing the exponential increase in χ* as a function of FPU in the case of dynein, noting that such force dependence of T eff was reported in a driven colloidal particle system 37 . We need to understand both quantitative and qualitative behaviour of T eff to calibrate force from the fluctuation unit χ in future. Melanophore observation by microscopy. After 1 day of culture, melanophores were observed by bright-field microscopy at room temperature. Prior to observation, the melanophores were cultured in L15 medium containing 10 μM α-melanocyte-stimulating hormone for 30 min, then dispersed by washing with L15 medium without α-MSH. After adding epinephrine (final concentration: 10 μM) to the dish, images were acquired using a 100× objective lens and a sCMOS camera at 800 fps. To observe the motion of melanosomes over a wide range, a 0.5× lens was used.

Methods
The centre position (X) of each melanosome was determined from the recorded images using ImageJ software (National Institutes of Health, Bethesda, MD, USA) 38 . We focused on the displacement along the direction of melanosome motion (Fig. 2a). Data were collected from 62 melanosomes from two individual melanophores in the absence of ciliobrevin; 59 melanosomes from two melanophores with 10 μM ciliobrevin; 61 melanosomes from three melanophores with 20 μM ciliobrevin; 51 melanosomes from two melanophores with 30 μM ciliobrevin; and 48 melanosomes from two melanophores with 40 μM ciliobrevin. The accuracy of the position measurement was verified using 300-nm latex beads. The mean ± standard deviation of the position of the beads securely attached to the glass surface was σ = 8.3 ± 1.2 nm (n = 4 beads). theoretical background. Microscopically, the motion of a cargo (x: position) at a constant velocity, for example, is described by the model other m B where γk T 2 B is thermal noise acting on the cargo (k B : Boltzmann constant, T: temperature of the environment), ξ is Gaussian noise with 〈ξ (t) ξ (t′)〉 = δ (t−t′), where < > denotes the time average over the time course, U m is an interaction between cargo and motors, and U other is the other interactions acting on the cargo. As assumed in the Stokes formula, γ = 6πηr (r: radius of the cargo, η: viscosity of the cytosol). Note that when the environment of a cargo is viscoelastic, it is suggested that a viscoelastic memory for the friction term is added to equation (1) [39][40][41][42] , instead of U other : www.nature.com/scientificreports www.nature.com/scientificreports/ The coarse-grained model of equation (1), which represents the motion of the cargo (X) for the time scale long enough to be compared with the relaxation time (∼20 ms), is phenomenologically derived from equation (1) based on a previous report 43 : The microscopic effect of U m and U other contributes to the effective viscosity Γ, effective noise Γk T 2 B eff , and the drag force F. The statistical properties of Gaussian force noise in equation (3) were determined from the power spectrum densities of X for a small f (the frequency) ( Supplementary Fig. S1) and the Gaussian distribution of P(ΔX), where ΔX = X(t + Δt) − X(t). Subsequently, the fluctuation unit χ was considered for the macroscopic model (equation (3)).
Calculation of χ based on the fluctuation theorem. The fluctuation unit χ FT , which was introduced in our previous studies [11][12][13] , is defined as FT from the distribution P(ΔX) of the displacement ΔX = X (t + t) − X (t) (e.g. Fig. 2b, inset). When P(ΔX) is fitted by a Gaussian function, neglecting non-typical fluctuation (extreme events): where the fitting parameters a and b correspond to the variance and mean of the distribution, respectively (e.g. Fig. 2c), χ FT is practically calculated as b a 2 / (6) χ = Because ΔX < 0 could not be observed for a large Δt because of the limitations of our experiment, P(ΔX) for ΔX > 0 was fitted with a Gaussian function to calculate the ratio P(ΔX)/P(−ΔX) in equation (4). χ is compared with χ FT for small Δt (Supplementary Fig. S7).
χ is calculated for each P(ΔX) for various intervals Δt from 1.25 to 37.5 ms (e.g. Fig. 2d). The relaxation time of χ in Fig. 2d represents the time scale by which the motion of X is described with the coarse-grained model (equation (3)). The microenvironment around the vesicle, especially its viscoelastic nature, affects the relaxation time as well as the enzymatic cycle time. Note that the Δt dependence of χ resulted from the non-linear behaviour of a for small Δt (<20 ms) (Fig. S8) rather than the linear behaviour of b for all Δt. For small Δt (<20 ms), a might not go below σ 2 where σ is the magnitude of the intrinsic noise of the experimental system (see "Melanophore observation by microscopy" for σ).
The converged value (χ* = χ at Δt = 36.25 ms) (Fig. 2d) of χ was related to the drag force (F) acting on a cargo according to the equation B eff ⁎ where k B is the Boltzmann constant, and T eff is the effective temperature, which is a generalised temperature in a non-equilibrium system 37,44,45 . A prior study 46 is helpful for derivation of the functional form of equation (7) from equation (3). Previous experiments 13 suggested that Note that equation (6) can be interpreted as χ = v/D, where v and D are the velocity and diffusion coefficient of a cargo. Then, χ is considered the inverse of the randomness parameter introduced in the reference 47 . Smoothing and affinity propagation. A smoothing filter was applied to the values of χ to reduce variation in the raw data for χ as a function of Δt (Fig. 2d). We used the averaging filter t t t t ( ) ( ( 1 25ms) ( ) ( 125 ms))/3, The data for χ ([ciliobrevin] = 0 μM) before applying the filter are shown in Supplementary Fig. S9. Affinity propagation 48,49 , an exemplar-based clustering method that does not require the number of clusters, was then adopted to cluster the smoothing-filtered two-dimensional data (χ*, χ m ), where χ m is the mean value of χ from Δt = 1.25-36.25 ms. The method was applied using the ' APCluster' function of R software 49 . The clustering was stable for the wide range of values for the parameter (q) (Fig. 3c). power spectrum. The power spectrum (S(f)) of the position X of a melanosome moving at a constant velocity was calculated as www.nature.com/scientificreports www.nature.com/scientificreports/ where f is the frequency, and τs = N w /800 s; 800 fps is the recording rate for the pause motion, and N w (=256) is the window size (plotted in Supplementary Fig. S1). S(f) for small f showed Brownian noise behaviour as S(f) ∝ f −2 .

Data Availability
Data supporting the findings of this study are available within the article and its Supplementary Information files and from the corresponding author on reasonable request.