Collective cell migration during human mammary gland organoid morphogenesis

Organ morphogenesis is driven by cellular migration patterns, which become accessible for observation in organoid cultures. We demonstrate here that mammary gland organoids cultured from human primary cells, exhibit oscillatory and collective migration patterns during their development into highly branched structures, as well as persistent rotational motion within the developed alveoli. Using high-resolution live-cell imaging, we observed cellular movement over the course of several days and subsequently characterized the underlying migration pattern by means of optical flow algorithms. Confined by the surrounding collagen matrix, characteristic correlated back-and-forth movements emerge due to a mismatch between branch invasion and cell migration speeds throughout the branch invasion phase. In contrast, alveolar cells exhibit continuous movement in the same direction. By modulating cell–cell adhesions, we identified collective migration as a prerequisite for sustaining these migration patterns both during the branching elongation process and after alveolus maturation.


INTRODUCTION
Complex multicellular morphogenetic processes, such as embryogenesis, organogenesis, and tumorigenesis, are enabled by collective cellular behavior [1][2][3] and orchestrated by biochemical 4,5 and mechanical 6,7 cues.Throughout a wide variety of processes, such as wound healing, branching morphogenesis, and tumorigenesis, collective migration is a key factor in enabling the formation of functional structures. 8In this case, long-ranged cellular patterns can emerge if correlated cell behavior is promoted by coupling of neighboring cells via cell-cell adhesions. 9his coupling leads to an integration of the individual traction forces and, thus, to a long-ranged stress field within the epithelial layer giving rise to high-order phenomena within the tissue. 10n order to recreate collective cell migration in vitro, cells are commonly cultivated on spatially confining 2D biopolymer patterns 9,[11][12][13][14] or inside simplified artificial 3D topologies. 15,16Within such confinement, collective cell migration is steered by the geometrical boundary conditions.It has been shown that cells on strip-like patterns exhibit a directed outgrowth as long as the strip diameter is narrower than the specific correlation length of the cells. 9In comparison, on circular and ring-shaped patterns with dimensions below, the correlation length cells perform rotational migration as long as they are coupled by cell-cell adhesions. 17Bridging the gap toward more complex systems, studies on multicellular spheroids have revealed complex rotational motion induced by the curvature 18 that resembles the coherent angular motion previously found for small numbers of breast epithelial cells, 19 which has been linked to basement membrane assembly, 20 while cells under quasi-3D tubular confinement are reported to invade into the tube as a result of a rotational but directed migration. 15In organoid systems, this transition of linear to rotational motion has been shown to result from the isotropification of the tension field by the migrating cells. 21uch collective migration persists until the cell density increases above a critical threshold or the cell-cell adhesions maturate, which ultimately leads to a jamming of the cells and an inhibition of cell migration. 6,22Within biological tissues in 3D migrating cells interact strongly with the surrounding extracellular matrix (ECM), thus leading to remodeling and invasion of the matrix. 3,23During branching morphogenesis of human mammary gland organoids, the mechanical interaction of the cells with their extracellular matrix ultimately leads to the formation of a mechanically stable collagen cage surrounding the extending organoids. 24,25Such cell-ECM interactions and matrix remodeling play a central role in organ development and rely on the detailed nature of the collective cell migration patterns and their spatial-temporal coordination.
Herein, we report on collective cell migration during the morphogenesis of primary human mammary gland organoids.We used long-term confocal live-cell imaging to observe the fluorescently labeled nuclei of individual cells over time.Analysis of the movement by means of automatic and manual single cell tracking as well as the calculation of optical flow fields revealed different dynamics in various stages of development.In contrast to persistently rotating alveoli, the internal flow field within the spatially confined and elongated branches exhibits a collective back-and-forth migration of the cells, while the tip cells in the front of the branches largely show uncoupled migration behavior.By inhibiting cell-cell adhesion, matrix degradation, and cell contractility, we demonstrate that the collective cell patterns arise as a result of the relation between tip invasion speed into the ECM, the individual cell migration velocities in the branch and the subsequent localized jamming events.

Cells are highly motile within human mammary gland organoids
Human mammary gland organoids were prepared by seeding singularized human primary cells in floating collagen gels.Over the course of two weeks, single cells grew into highly branched and multicellular organoids with a diameter of up to 2 mm [Figs.1(a) and Movie S1].The growth of the organoids was finalized by the formation of spherical alveoli at the end of elongated branches, which arise due to an isotropification of the internal stress field within the branches. 21uring branch elongation, their expansion is accompanied by a longrange and anisotropic deformation field in the surrounding collagen network, which is a result of the forces exerted by the organoids by means of internal cell migration.These plastic deformations ultimately lead to the formation of a mechanically stable collagen cage surrounding and stabilizing the organoids. 24This creates as spatial confinement for the cells within the organoid thus restricting and defining the cellular migration within the branches.
To analyze the resulting modes of cellular migration, we fluorescently labeled the cell nuclei using SiR-DNA stain and conducted live cell imaging of the organoids throughout all developmental stages.Due to heterogeneities in organoid morphology and age, the branch lengths ranged from 79 to 671 lm, and the observed branch widths were between 15 and 63 lm.The analyzed alveoli had diameters between 78 and 124 lm and were all hollow in the center.The cells within organoids appeared to be highly motile throughout all developmental stages from day 7 until day 14 [Fig.S1(a)].The average velocity across the analyzed cylindrical branches was around 0.14 lm/min, reaching maximal values of up to 0.7 lm/min.No statistically significant differences in average cell velocities were found between different organoid ages.Furthermore, the average and maximal velocities appeared to be independent of branch length and width [Fig.S1(b)].We observed that the 3D cell migration within the branches was highly correlated with branch geometry.Cells within cylindrical branches moved mainly parallel to the branch axis and only exhibited a negligible migration orthogonal to the branch axis [Fig.1(b)].In comparison, the cells comprising the alveoli exhibited a rotational migration perpendicular to the main branch axis with speeds of up to 1 lm/min [Figs.1(b) and 1(c) and Movie S2].The cumulative displacement of single cell tracks demonstrated a linear behavior, thus indicating continuous nonstop motion at the single cell level [Fig.

1(d)]
. There was no prevailing chirality in rotation, with 45% of the analyzed alveoli (n ¼ 20) rotating clockwise and 50% rotating counterclockwise.Only in 5% of the analyzed events a reversal of the rotation direction was observed.Therefore, once a rotational migration within the alveoli had been established, the cells persistently rotated in one direction.The observed migration modes in the elongated branches and spherical alveoli resembled those found on 2D-patterned surfaces having corresponding geometries. 9,11However, in contrast to these previous results in 2D systems, the cells in the branches did not exhibit a persistent migration in one direction but rather changed their migration direction over time.To reveal the temporal migration patterns, the optical flow fields obtained by analysis of the fluorescence live cell imaging data were projected and temporally resolved [Fig.1(e)]: The resulting histograms of the parallel velocities for each frame were color-coded and displayed as a single line.Subsequent stacking of the individual lines illustrates the temporal development of the parallel velocity distribution within a branch.It, thereby, became apparent that the cells were not continuously migrating in only one direction, but rather exhibiting several phases, during which most cells within a branch collectively migrated either inwards or outward [Fig.1(f)].These back-and-forth migration patterns took place over several hours and were interrupted by comparatively short phases without motion and the subsequent migration reversals.The 52 analyzed turnaround events showed an average time between direction reversals of 108 min [Fig.1(g)].

Cells migrate collectively in space and time
A spatiotemporal analysis of the velocity fields was performed in order to further investigate the collective properties of the observed migration patterns.Given velocity histograms for a given branch exhibiting multiple distinct peaks at single time points, coloring pixels contributing to the individual peaks revealed clusters of collectively migrating cells [Fig.2(a)].By manually tracking single cells over multiple back-and-forth phases, it was confirmed that cells were confined to individual branches and thus, the organoid branches did not act as a large, interconnected flow system [Fig.2(b)].Instead, the migration patterns appeared to emerge locally and independently of other branches, leading to a non-synchronized movement of branches belonging to the same organoid [Fig. 2

(c) and Movie S3].
In order to visualize the spatiotemporal development of the velocity field within a branch, parallel cell velocities were color-coded and subsequently averaged along the orthogonal branch axis, thus taking advantage of the cylindrical symmetry of the branches and the neglectable orthogonal components.Combining multiple of these one-dimensional representations of the branch velocity profile for different timesteps yielded velocity kymographs [Fig.2(d)].During phases of stable inwards or outward movement, large sections of the branches were moving in unison, which confirmed the temporal stability of the previously observed clusters, whereas uncorrelated movement was observed in the comparatively short turnaround phases [Figs.Given that antibody staining of f-actin indicated strong coupling between the cells in a branch, antibody staining of E-Cadherin (which is a transmembrane protein known for mediating cell-cell junctions) was performed to further understand the role of collective cell migration in branching morphogenesis [Fig.3(a)].The branch cells expressed high levels of localized E-Cadherin, which enables cellular crosstalk between neighboring cells and has been reported to be crucial to collective migration.We hypothesized that, in the densely packed branches, collective migration plays a key role in sustaining any cell movement.To test this hypothesis, the E-Cadherin antibody HECD1

Mismatch between cell velocity and tip invasion induces back-and-forth migration
To further understand these findings, the structural properties of the surrounding ECM must be considered.Notably, the organoids' branches were surrounded by a cage-like collagen structure as previously reported. 24This mechanically stable encasing, which had pores much smaller than average cell sizes created a barrier for the cells.As a result, cell migration within human mammary gland organoids should be understood as a form of cell migration under geometrical confinement.However, the spatial confinement of the branches in this system is not fixed.Instead, it is constantly being changed by the dynamic invasion process into the ECM.This process is facilitated by the tip cells, whose dynamics have been shown to be decoupled from the stalk cell migration.To investigate the effect of dynamic spatial confinement on cell migration patterns, we quantified the branch tip velocity and relative branch extension during organoid growth [Fig.4(a)].When comparing the determined invasion velocities of the branch tip with the migration velocity of the cells in the branch stalk, a clear mismatch became apparent.Throughout every inwards or outward directed migration phase, there were a substantial number of cells within a branch, which exceeded the tip invasion velocity.As a result, local This region of high density brought cell migration to a halt and led to a direct reversal of the stalk cells, thus promoting the establishment of a new inward-directed migration phase during which the tip cells continue to pull on the surrounding collagen matrix.Importantly, by analyzing the tip/branch thickness ratio a statistically significant increase in relative tip thickness became apparent, further supporting the claim of increased cell crowding in the tip region at the time of migration reversal [Fig.4(e)].
However, the ECM invasion process does not solely rely on mechanical remodeling; it also requires enzymatic matrix digestion. 26Cells secrete matrix metalloproteinases (MMPs) that degrade the surrounding collagen fibers, thus allowing the branch to increase in length.In order to demonstrate the impact of the observed velocity mismatch between invasion and migration velocities on branch contractions, drug treatments with Marimastat were performed on the organoids so as to inhibit matrix metalloproteinases.Cells were subsequently imaged for more than 70 h immediately after the drug was added in order to investigate the effect of decreased invasion speed on internal cell migration (see the supplementary material Movie S6).
During the initial hours of treatment, the cells still exhibited the characteristic back-and-forth patterns observed in the control.However, even though the cells continued to proliferate, the branch no longer significantly increased in length and cell velocities simultaneously decreased after 15-20 h [Fig.4(c)].Subsequently, cell migration came to a halt and, after 50 h of treatment, the branch tip slowly started to retract, indicating decreasing tension within the organoid branch.

DISCUSSION
It has been previously reported that cells on 2D-patterned substrates exhibit migration patterns that are dictated by the geometrical boundary conditions.While cells seeded on strip-like patterns migrate mainly in parallel with respect to the main axis, circular patterns lead to the emergence of rotational migration. 11We found that similar migration modes can be observed during the growth of human mammary gland organoids in a three-dimensional matrix, in which the geometrical confinement is created by the dynamic formation of a self-built collagen cage around the organoids.As the branches span at most 5-7 cells in diameter during the branching morphogenesis phase in early organoid development and show a high degree of rotational symmetry with neglectable orthogonal velocities, we were able to quantify the most important aspects of the large-scale cell migration by analyzing the 2D projections of confocal image stacks.Here, the correlation length of cells migrating in extending branches was comparable to those found in 2D.In contrast to the previously described 2D and 3D experiments, 11,15 the spatial confinement of extending organoid branches constantly changes over time and is remodeled by the tip cells themselves. 24Nevertheless, the matrix remodeling process at the invasion site imposes a limit on the maximum possible invasion speed of the branch tip.As a result, there is a significant mismatch between stalk migration velocities in the proximal part of the branch and the tip extension.Consequently, fast outward migrating cells will periodically jam at the branch tip, thus slowing down and ultimately changing their migration direction.Over time, this effect gives rise to an oscillatory back-and-forth migration mode, which is at the heart of branching morphogenesis in basal human mammary gland organoids.The dynamic coupling of the motion and force exertion of the branches to the ECM results in the mechano-plastic remodeling of the ECMand leads to branch extension.This bidirectional interplay is not possible in Matrigel, due to the lack of mechanical plasticity, which might also explain different branching modes observed.While in pure Matrigel murine, mammary epithelial fragments of mouse cell lines have shown to form simple cyst-liked protrusions without a clear mechanical interaction of the cells with the matrix, 27,28 already enriching the Matrigel with 30% prepolymerized collagen leads to branched structures of higher complexity. 29The mechanical plasticity of these networks is unknown and it remains to be determined what kind of cellular migration dynamics drive the branching morphogenesis in this model systems.To this end, the observed cellular migration patterns in the morphogenesis of human mammary gland organoids presented here may, thus, serve as a benchmark to further understand the morphogenesis of other mammary gland organoid systems as well as those of other organs with possibly even more complicated geometries.

MATERIALS AND METHODS Organoid cultivation
Primary human mammary gland organoids were prepared following a previously published protocol (Linnemann).Freshly isolated human mammary gland epithelial cells from healthy women undergoing reduction mammoplasty were embedded in freely floating collagen I gels (Cornings).Branched organoids were formed by EpCAMþ/ CD49fþ/CD10þ basal cells.The collagen concentration was set to 1.3 mg/ml.For the first five days, cells were cultivated in mammary epithelial growth medium (promocell MECGM) mixed with 3 lM Y-27632 (Biomol), 10 lM Forskolin (Biomol) and 0.5% FBS.After five days, the MECGM was enriched with 10 lM Forskolin only.Organoids were grown at 37 C at an atmosphere of 5% CO 2 and 3% O 2 .Finally, organoids were imaged between days 6 and 12 to analyze the branching morphogenesis phase, while older organoids with fully formed alveoli were used to describe the rotational motion in the alveoli.

Confocal live-cell imaging
Long-term live-cell imaging was conducted with a Leica SP8 lightning confocal microscope.Data were collected with Leica Application Suite X v. 3.57.The CO 2 level (5%) and both temperature (37 C) and humidity (80%) were controlled using a gas incubation system (Ibidi).In order to visualize the nuclei, cells were incubated with 10 lM sirDNA (Spirochrome AG) 3 h prior to the measurements.Labeled nuclei were excited at a wavelength of 633 nm.The resulting signal was collected at an emission maximum of around 674 nm through a HCX PL APO 10Â/0.40CS dry objective.Dependent on the experiment the observation time varied between 12 and 72 h with time steps between each frame ranging from 5 to 30 min.

Inhibitor treatments
Inhibitors were added at various developmental stages to the cell culture media.In particular, the media was enriched with the inhibitors directly prior to the experiment, which was then present throughout the whole measurement.Matrix degradation was inhibited by 10 lM Marimastat (Merck), and cell-cell adhesion was lowered by blocking E-Cadherin via HECD1 (Abcam) in a dilution of 1:25.

Immunofluorescence
Collagen gels were first washed with PBS solution, subsequently fixed with 4% paraformaldehyde for 15 min and finally rinsed twice with PBS solution again.The cell membranes were permeabilized with 0.2% Triton X-100 and blocked with 10% goat or donkey serum in 0.1% BSA.Regarding immunofluorescence stainings, the samples were incubated with primary and secondary antibodies diluted in 0.1% BSA.DAPI was used to visualize the cell nuclei, HECD1 (Abcam) was used to label E-cadherin and Phalloidin conjugated with Atto 647 (Sigma) was used to visualize the actin network.Alexa 488 (Life Technologies) was used as the secondary antibody for the E-Cadherin.

Optical flow
Optical flow analysis was performed using the Gunnar-Farneback optical flow implementation in the Python 3 OpenCV library.Initial image processing and data analysis were performed using Python's scikitimage, SciPy, and OpenCV as well as ImageJ.The pixel displacements in between frames were determined by applying the Gunnar-Farneback optical flow implementation of the Python 3 OpenCV library to 2D projections of confocal microscopy stacks, thus enabling the calculation of a velocity vector for every pixel.First, noise reduction was performed on the images using a median filter.Second, regions of interest (ROIs) were defined together with a branch axis for every region.Third, the image sequence was analyzed using the optical flow algorithm and the resulting velocity vectors were split into parallel and orthogonal velocity components with respect to the defined branch axis.

Single cell tracking
Single cell tracking was performed manually for 2D projected branches using FIJI.3D-tracking of nuclei in alveoli was achieved using Leica AIVIA Version 10.0.9/2021.

Calculation of correlation length and time
Correlation functions of the mean centered kymographs were computed using the numpy correlate function.For temporal correlation, the 1D correlation was calculated for every position along the branch axis and all obtained curves for a given ROI were subsequently averaged.Similarly, for the spatial correlation, single curves were computed for every timepoint and then averaged across all positions.

Statistics and reproducibility
For statistical testing, p-values were calculated using two-tailed wilcoxon ranksum analysis as implemented Python's SciPy library.All experiments were repeated for at least three independently prepared organoid batches.The exact numbers of observed phenomena are reported in the figure legends.

SUPPLEMENTARY MATERIAL
See the supplementary material for additional data as referenced in the article as well as live microscopy data and analysis movies.
2(f) and S2].Temporal and spatial autocorrelation functions of the mean-centered kymographs were calculated to further examine the collectiveness of the underlying cell migration [Fig.2(e)].Examining the first zero-crossings of the resulting curves revealed average correlation times of 2 (61) h and correlation lengths of between 65 and 105 lm [Fig.2(g)].

FIG. 1 .
FIG. 1. Experimental setup.(a) Bright-field and SiR-DNA stained nucleus signal of an exemplary mammary gland organoid (b) optical Flow reveals direction reversals in branches and rotations in alveoli.(c) Alveoli persistently rotating at nearly constant velocities without changing the rotation direction.(d) Single alveolar cells showing linear cumulative displacement, indicating persistent continuous migration at the single cell level.(e) By splitting the calculated velocities with respect to the branch axis and stacking the resulting histogram for all time steps the temporal change in the velocity distributions can be visualized.(f) The time-resolved histogram for a cylindrical branch reveals a distinct back-and-forth migration pattern parallel to the branch axis.In comparison, the orthogonal velocities were negligible.Scale bars, 100 lm.(g) Histogram of the distribution of turnaround times for n ¼ 52 turnaround events in n b ¼ 8 independent branches.The average turnaround time of the analyzed events is t ¼ 108 min.

FIG. 2 .
FIG. 2. Spatio-temporal correlation analysis.(a) Peaks in the velocity distribution (green, magenta) correspond to spatially connected areas within the branch.(b) Individual single cell tracks show that single cells commuted within single branches and did not migrate into neighboring branches over time.(c) Branches of an organoid exhibit an oscillatory mean velocity.Their frequency and velocity are not correlated but rather independent off another.(d) For migration analysis in branches, the parallel component was averaged perpendicular to the branch axis.The resulting velocity line profiles were stacked for all timepoints revealing the spatiotemporal development of the oscillations.(e) Representative curves for the correlation analysis in space and time.The first zero crossing of correlation functions computed on mean-centered velocity kymographs was used to determine the correlation lengths and times.(f) Representative velocity kymograph of an extending branch.In red, outward directed phases appeared repeatedly after those of inwards oriented migration colored in blue.Large sections of the branch move in unison.(g) Distribution of correlation lengths and times in n l ¼ 31 and n b ¼ 28 branches.All scale bars 100 lm.

FIG. 3 .
FIG. 3. Cell-cell adhesion and contractility.(a) Close-up of an organoid branch stained for cell nuclei (DAPI), F-Actin (Phalloidin), and E-Cadherin (HECD1).Actin bundles that span several cells and expression of E-Cadherin across all cell membranes indicate high cellular coupling.(b) Velocity kymograph after E-Cadherin inhibition with HECD1.Initially, cells still showed back-and-forth like behavior, which comes to a standstill after approximately 15 h.(c) Mean cell velocity, tip velocity and tip extension over time after HECD1 treatment.As the mean cell velocities decreased, there was no further tip extension in the branches.Notably, the 25-75 percentiles also converged to zero indicating that virtually all cells in a branch stop migration.All scale bars, 100 lm.

FIG. 4 .
FIG. 4. Mismatch between cell invasion and cell migration.(a) Mean cell velocity, tip velocity and relative tip extension over time.The tip velocity was predominantly within the 25-75 percentiles of the cell velocities.(b) Cell nuclei with optical flow field and bright-field images of an organoid branch.Cells jammed at the front and subsequently reversed their direction (c) metalloproteinase inhibition upon treatment with Marimastat led to a stagnation in matrix invasion and, ultimately, highly reduced cell velocities.(d) Schematic representation of the proposed mechanism for the emergence of back-and-forth-oscillations.Scale bar, 100 lm (b).(e) Comparison of tip/branch thickness ratio 1 h before and at the point of direction reversal for ten turnaround events in five branches (colors) of four different organoids (markers).The tip thickness increased relative to the branch statistically significant with a p-value of 0.006 91.