Microtubule associated protein WAVE DAMPENED2-LIKE (WDL) controls microtubule bundling and the stability of the site of tip-growth in Marchantia polymorpha rhizoids

Tip-growth is a mode of polarized cell expansion where incorporation of new membrane and wall is stably restricted to a single, small domain of the cell surface resulting in the formation of a tubular projection that extends away from the body of the cell. The organization of the microtubule cytoskeleton is conserved among tip-growing cells of land plants: bundles of microtubules run longitudinally along the non-growing shank and a network of fine microtubules grow into the apical dome where growth occurs. Together, these microtubule networks control the stable positioning of the growth site at the cell surface. This conserved dynamic organization is required for the spatial stability of tip-growth, as demonstrated by the formation of sinuous tip-growing cells upon treatment with microtubule-stabilizing or microtubule-destabilizing drugs. Microtubule associated proteins (MAPs) that either stabilize or destabilize microtubule networks are required for the maintenance of stable tip-growth in root hairs of flowering plants. NIMA RELATED KINASE (NEK) is a MAP that destabilizes microtubule growing ends in the apical dome of tip-growing rhizoid cells in the liverwort Marchantia polymorpha. We hypothesized that both microtubule stabilizing and destabilizing MAPs are required for the maintenance of the stable tip-growth in liverworts. To identify genes encoding microtubule-stabilizing and microtubule-destabilizing activities we generated 120,000 UV-B mutagenized and 336,000 T-DNA transformed Marchantia polymorpha plants and screened for defective rhizoid phenotypes. We identified 119 mutants and retained 30 mutants in which the sinuous rhizoid phenotype was inherited. The 30 mutants were classified into at least 4 linkage groups. Characterisation of two of the linkage groups showed that MAP genes–WAVE DAMPENED2-LIKE (WDL) and NIMA-RELATED KINASE (NEK)–are required to stabilize the site of tip growth in elongating rhizoids. Furthermore, we show that MpWDL is required for the formation of a bundled array of parallel and longitudinally orientated microtubules in the non-growing shank of rhizoids where MpWDL-YFP localizes to microtubule bundles. We propose a model where the opposite functions of MpWDL and MpNEK on microtubule bundling are spatially separated and promote tip-growth spatial stability.


Introduction
Filamentous cells, such as root hairs of vascular plant sporophytes and rhizoids on vascular and non-vascular plant gametophytes, form at the interface between plants and soil. They carry out rooting functions, such as anchorage, water and nutrient uptake, and interact with microorganisms [1]. Their tubular shape is key to this function because anchorage is defective in Marchantia polymorpha mutants with defective rhizoid morphology (S1 Fig) and in Arabidopsis thaliana mutants with defective root hairs [2]. Filamentous rooting cells elongate by tipgrowth, a mechanism where growth is stably restricted to a small domain of the cell surface from which the tubular projection grows. Root hairs and rhizoids form as straight cylinders when growing in the air. However, when growing through soil substrates their growth direction continually changes as the tip manoeuvres around objects in the soil.
Tip growth involves the delivery to the apical dome of secretory vesicles containing cell wall material in land plants [3], including liverworts (S3A Fig). The stability of the position of the growing tip is controlled by external factors, such as soil particles, and internal factors such as microtubules [4][5][6]. A common feature of the microtubule array in tip-growing cells of all land plants is its organization in the shank of the cell, proximal to apical dome. In the non-growing shank, bundled cortical microtubules and endoplasmic microtubules are oriented longitudinally and grow towards the tip [7][8][9][10][11][12][13]. The conservation of this organization of microtubules suggests that parallel, longitudinal microtubules growing from a basal position in the nongrowing shank to the apical dome is required for tip-growth.
A function of the microtubule cytoskeleton in stabilizing tip-growth has been inferred from observations of plant tip-growing cells treated with drugs that inhibit microtubule assembly and microtubule organization. Tip-growing cells exposed to drugs that either stabilize or destabilize microtubules display wavy tubular morphologies, suggesting that the growth region is unstable and shifts laterally from an apical position to a subapical position [6,11]. Branching root hairs are also observed following the application of microtubule stabilizing drugs, suggesting that dynamic microtubules are required to maintain a single growth point. These observations suggest that the regulation of microtubule dynamics is required to stabilize the position of the growing region in tip-growing cells. Consistent with this hypothesis, loss of function mutations in genes coding for microtubule associated proteins (MAPs) that promote microtubule growth or destabilize microtubules can lead to loss of tip-growth spatial stability. In A. thaliana, MICROTUBULE ORGANIZATION1 (MOR1) promotes microtubule growth and bundling, and the temperature-sensitive mutant mor1-1 forms wavy root hairs at the restrictive temperature [14]. By contrast, ARMADILLO-REPEAT KINESIN1 (ARK1) destabilizes microtubules by promoting the transition of microtubules from the polymerizing to the depolymerizing state, and root hairs are wavy in the ark1 loss of function mutant [15]. In the liverwort Marchantia polymorpha, the microtubule destabilizing MAP NIMA-RELATED PROTEIN (MpNEK) is required to maintain a stable position of the growing region in rhizoids [16]. However, no microtubule stabilizing MAP has been reported to function in stabilizing tipgrowth in the liverwort model to date.
We hypothesized that microtubule stabilizing MAPs may also be required in tip-growing cells of liverworts to stabilize the position of the growth point. To identify microtubule stabilizing MAPs that stabilize tip-growth in liverworts, we screened T-DNA-and UV-mutagenized Marchantia polymorpha populations for mutants with wavy rhizoids that are morphologically similar to oryzalin-treated or Taxol-treated rhizoids. We isolated two nek mutants, demonstrating that mutants with defective microtubule organization were identified in the screen [16]. Moreover, we identified four mutant lines with T-DNA insertions in the gene encoding the MAP, WAVE DAMPENED2 LIKE (WDL). We show that MpWDL localizes preferentially to microtubules in the non-growing shank, where it promotes the bundling and parallel organization of longitudinal microtubules. Together, these results indicate that MAPs with opposite effects on microtubule bundling-MpNEK depolymerizes while MpWDL bundles-are required to regulate the organization of the microtubule array that stabilizes the site of tipgrowth in liverwort rhizoid.

Mutant screen identified 30 mutants with putative defects in microtubule dynamics
To identify proteins that are active in microtubule-mediated stabilization of the apex during tip-growth, we screened mutagenized Marchantia polymorpha plants for mutants with defective rhizoids. To define the expected rhizoid phenotypes of mutants with defective microtubule organization, we first defined the morphology of wild type rhizoids with impaired microtubule dynamics caused by oryzalin-treatment, which inhibits microtubule polymerisation at the plus-end, and Taxol-treatment, which stabilizes microtubules by inhibiting depolymerisation at both plus and minus ends. Rhizoids were grown in the presence of the microtubule-depolymerizing drug oryzalin (Fig 1) and observed with a confocal microscope. At lower oryzalin concentrations (0.1 μM-0.5 μM), rhizoids developed a wavy phenotype compared to rhizoids grown in control conditions which were straight. At higher oryzalin concentrations rhizoid sinuosity increased and some rhizoids branched. Similarly, growth of rhizoids in the presence of 3.3 μM Taxol resulted in a wavy rhizoid morphology and branched rhizoids. Taken together, these data indicate that that microtubules are not required for rhizoid growth per se.
However, it suggests that a dynamic network of microtubules is required to stabilize the position of the apical domain of the tip growing rhizoid.
To identify genes that control the stability of the growing apex in rhizoids we carried out forward genetic screens for mutants that have rhizoids resembling oryzalin-and Taxol-treated rhizoids. 120,000 UV-B mutagenized plants and 336,000 T-DNA transformed lines were generated and screened for wavy rhizoid phenotypes. We initially selected 97 independent mutants with wavy rhizoids on the mature gametophyte from the 120,000 UV-B mutagenized plants in the initial screens. We also selected 22 wavy rhizoid mutants from the 336,000 T-DNA transformed lines [17]. Of these 119 mutants, 30 (22 UV-B-induced wavy mutants and 8 T-DNA-induced wavy mutants) were retained for further analysis. The remaining 89 mutants were not characterized further in this study, either because their mutant phenotype could not be observed in subsequent vegetative generations or because their growth was too severely stunted.
To quantify the defective phenotypes of rhizoids in these mutant lines, rhizoids that formed on the upper side of gemmae, which are vegetative propagules genetically identical to the plant from which they form, were imaged two days after plating. Two parameters were measured: rhizoid diameter (thickness) and rhizoid sinuosity (waviness) ( Table 1). An Anova test indicated that wild type and mutants could be discriminated by rhizoid sinuosity and rhizoid diameter (p kruskal-wallis < 2.10 −6 ). Therefore, these two phenotypic parameters were used to categorize the mutant lines.
A first group of 7 mutants (Group 1) developed rhizoids whose diameter was greater than those of wild type and sinuosity was equal or greater than in wild type. The thicker rhizoid phenotype suggests that the domain of cell growth at the apex may be larger in rhizoids of these mutants than in wild type. Furthermore, an additional phenotype was observed in each mutant with thicker rhizoids and greater sinuosity: they developed gaping air pores on the dorsal epidermis in the mature gametophyte; openings were elongated instead of disc-shaped in wild type (S2 Fig). These phenotypic similarities shared by these mutants suggested that each line might harbour mutations in the same gene, i.e. they might be allelic. These mutants were not characterized further in this study because the increase in rhizoid diameter was not observed in wildtype rhizoids treated with microtubule stabilizing or microtubule destabilizing drugs. Table 1. Rhizoid phenotypes of T-DNA and UV mutant lines fall in two categories. The mean values of rhizoid sinuosity and rhizoid diameter were measured in 10 mature rhizoids per genotype and are given +/-SD. Single or double stars indicate a Benjamini-Hochberg adjusted p-value from non-parametric Dunn test lower than 0.05 and 0.01 respectively. Phenotypes significantly different from wild type were highlighted to help visualize the independence of the large rhizoid and wavy rhizoid phenotypes. A second group of 13 mutants (Group 2) developed rhizoids whose diameter was identical to wild type but whose sinuosity was greater than in wild type. The rhizoid phenotype of this mutant class-the same diameter as wild type rhizoids, but more sinuous-suggests that these mutants polarize growth similarly to wild type but fail to stabilize the growth site to a constant position in the apical dome. We hypothesized that the mutations causing the Group 2 phenotype defined several linkage groups.

Multiple mutant alleles define MpNIMA-RELATED PROTEIN KINASE (MpNEK) as a regulator of stability of the site of tip growth
To test the hypothesis that Group 2 mutants defined multiple linkage groups, all 14 Group 2 mutants were crossed with each other. We analysed the segregation of rhizoid phenotypes in their progeny. The phenotype of progeny of F1 crosses between T-DNA and UV-B of Group 2 mutant lines-the same diameter as wild type rhizoids, but more sinuous-indicated that 3 of these 14 mutations constitute a single linkage group. Crosses between UV3.4, ST47-6 and VJ72 resulted in 100% mutant progeny ( Table 2). The absence of wild type progeny is consistent with the mutations being in the same gene or closely linked on the same chromosome. To identify the causative mutation in UV3.4 we sequenced the genomes of UV3.4 and independent non-allelic UV-B mutants. Then, we compared the mismatches between UV3.4 and the independent non-allelic UV-B mutants. Of the 190,843 mismatches in the initial UV3.4 pool of mismatches, six candidate mutations remained after filtering steps were carried out ( Table 3). Of those, only one mutation is a nonsense mutation and it is in the coding sequence of NIMA-RELATED PROTEIN KINASE (MpNEK) ( Table 3). To verify that the mutation in MpNEK causes the wavy rhizoid phenotype of UV3.4 we tested if there was a mutation in the MpNEK gene of the VJ72 mutant which is in the same linkage group as UV3.4. We first tested if the T-DNA insertion caused the wavy rhizoid phenotype. We scored hygromycin resistance among progeny of the VJ72 line (harbouring a hygromycin resistance gene on its T-DNA) backcrossed to wild type ( Table 4). Half of the progeny was hygromycin resistant, consistent with the segregation of a single T-DNA in this population. Of these progeny, all hygromycin resistant plants formed wavy rhizoids, indicating that the T-DNA co-segregated with the wavy mutant phenotype. These data are consistent with the hypothesis that a T-DNA insertion causes the wavy phenotype in mutant plants. To identify the gene mutated by the T-DNA insertion, sequences flanking the T-DNA from VJ72 were isolated by TAIL PCR (S1 Data). The flanking sequences indicated that the T-DNA was inserted in the promoter region of MpNEK 1760 bp 5' of the predicted start codon (Fig 2). Taken together these data indicate that the nonsense mutation in the in MpNEK gene of UV3.4 and a T-DNA insertion in the promoter of MpNEK in VJ72 cause a defective wavy rhizoid phenotype.  There is a single NEK protein encoded in the M. polymorpha genome, P. patens, three in Selaginella moellendorfii, six O. sativa and seven in A. thaliana. To determine the phylogenetic relationships between these NEK proteins we generated gene trees using maximum likelihood statistics ( Fig 2I). MpNEK is a member of a monophyletic group that contains land plant NEK proteins (bootstrap 0.995). Within this group MpNEK is a member of a monophyletic group that includes the single copy P. patens NEK protein and the three S. moellendorffii proteins. However, the bootstrap support for this monophyletic group is not strong. The similarity between the morphology of oryzalin and Taxol-treated rhizoids with the nek loss of function, is consistent with a role of MpNEK in modulating microtubule dynamics (Figs 1 and 2) required for stabilization of the apex during tip-growth as proposed by [16].
The identification of MpNEK in the screen validates our forward genetic approach to identify proteins that regulate microtubule dynamics and organization during tip-growth. We carried out no further characterization of MpNEK because MpNEK is already known to be involved in tip-growth stability and destabilize microtubules in growing rhizoids of M. polymorpha (16).

Multiple mutant alleles define WAVE DAMPENED-LIKE (WDL) as a regulator of stability of the site of tip growth
The phenotype of progeny from F1 crosses between T-DNA and UV-B of Group 2 mutant lines-with higher sinuosity than wild type but same diameter as wild type-indicated that 4 of these 14 mutations comprised a second linkage group. CR1, CR2, ST45-1 and ST33-5 T-DNA mutants were crossed to each other. The rhizoids of all (100%) of the F1 progeny were sinuous compared to the straight rhizoid of wild type ( Table 5). The absence of wild type progeny is consistent with the mutations being in the same gene or closely linked on the same chromosome. This suggested that each mutant line harbours an independent mutation in the same gene (Table 5). DNA sequences flanking the T-DNAs were identified by TAIL PCR. The flanking sequences demonstrated that the T-DNA of CR1, CR2 and ST45-1 were inserted into a gene encoding a protein that is homologous to the A. thaliana WAVE DAMPENED2 LIKE [17]. These mutant alleles were designated wdl-1, wdl-2 and wdl-3 respectively. TAIL PCR failed to amplify sequences flanking the T-DNA in the ST33-5 mutant. To identify the causative mutation in ST33-5, we sequenced the entire genome of ST33-5 and identified a T-DNA insertion site (S1 Data) by aligning sequencing reads against the reference genome and the T-DNA sequence. A T-DNA in ST33-5 was inserted in the promoter region of MpWDL 3436 bp upstream of the start codon. (Fig 3A). The presence of an insertion into the MpWDL promoter region of ST33-5 and the absence of wild type in the next generation when crossed to wdl-1, wdl-2 and wdl-3 mutants indicates that the sinuous rhizoid phenotype of ST33-5 results from a mutation in the MpWDL gene.
T-DNAs are inserted into the 7th exon and 7th intron of the MpWDL gene in wdl-2 and wdl-3 respectively (Fig 3A). Insertions in introns and exons often block the production of a stable transcript. To test if wdl-2 and wdl-3 produce a functional MpWDL transcript we measured steady state levels of MpWDL mRNA in wdl-2 and wdl-3 mutants. MpWDL transcript was not detectable in either wdl-2 or wdl-3 ( Fig 3B). This suggests that wdl-2 and wdl-3 are complete loss of function mutants. Furthermore, the location of the T-DNA insertions in wdl-2 and wdl- 3 at the predicted microtubule-binding domain TPX2 suggests that even if a truncated protein were produced, it would not be functional. To test the hypothesis that wdl-2 and wdl-3 are loss of function alleles, we transformed these mutants with wild type MpWDL genomic sequence that included 5 kb upstream of the transcriptional start site and the coding region (proMpWDL:MpWDL-YFP). The wild type rhizoid phenotype was restored by expressing proMpWDL:MpWDL-YFP in the wdl-2 and wdl-3 background (Fig 3J-3L). We conclude that the wavy rhizoid phenotype is caused by the loss of WDL function in wdl-2 and wdl-3.
A T-DNA is inserted in the promoter of MpWDL of the ST33-5 wdl mutant, with the right border of the T-DNA oriented 5' to the MpWDL locus. Insertions into the 5' regulatory regions where the right border of the T-DNA is 5' to the coding sequence have been reported to be associated with the overexpression of the sequences 3' of the T-DNA [17]. To test the hypothesis that MpWDL is expressed at higher levels in ST33-5 mutants than in wild type, we measured the steady state levels of MpWDL transcript. MpWDL transcript levels are 4-times higher in ST33-5 than in Tak-1 wild type (Fig 3B). Furthermore, there are no mutations in the coding sequence of MpWDL in ST33-5. These data suggest that the ST33-5 line harbours a gain of function MpWDL mutation in. ST33-5 was named wdl GOF-1 .
The fact that both loss of function and gain of function wdl mutants develop sinuous rhizoids suggests that loss of MpWDL activity and extra MpWDL activity had a similar effect on the stability of the apex during tip growth. To determine if the loss and gain of function mutants were phenotypically distinguishable, we measured waviness of the rhizoids by calculating their sinuosity. We could not distinguish between the sinuosity of the gain of function mutant wdl GOF-1 and the loss of function mutant wdl-3 (Fig 3C, 3F and 3I). Because the phenotypes of some wdl loss and wdl gain of function mutants are indistinguishable-like Taxol-and oryzalin-treated rhizoids at certain concentrations (Fig 1)-these data suggest that the MpWDL protein regulates microtubule dynamics required for the stabilisation of the apex during tipgrowth.

MpWDL-YFP localizes to microtubules and promotes the formation of a longitudinal array of parallel-arranged bundled microtubules in the shank of growing rhizoids
There is a single MpWDL encoding gene in M. polymorpha, 13 in P. patens, 6 in O. sativa and 8 in A. thaliana. MpWDL proteins are members of the TPX2 domain-containing microtubule binding proteins that contain the KLEEK motif (at position 352 in MpWDL) (Fig 3M). In A. thaliana, MpWDL proteins promote microtubule growth, and bind to and bundle microtubules in vitro [18]. We hypothesize that MpWDL modulates microtubule organization during the growth of the M. polymorpha rhizoid.
If MpWDL protein regulates microtubules in stabilizing the apex during tip growth, we hypothesized that MpWDL protein would bind to microtubules. To test this hypothesis, we transformed wild type M. polymorpha with a gene construct in which the MpWDL-YFP protein fusion was placed under the control of the endogenous MpWDL promoter (proMpWDL-MpWDL:YFP). We imaged YFP fluorescence in transformed lines. MpWDL-YFP localizes to interphase, spindle and phragmoplast microtubules in all gemma epidermal cells investigated (Fig 4A and 4B). Moreover, the MpWDL-YFP fluorescence distributes homogeneously along the length of microtubules, suggesting that MpWDL does not preferentially decorate the plusends of microtubules.
Microtubules are arranged in two distinct domains in the growing rhizoid. Dense bundles of microtubules run longitudinally along the shank of the rhizoid and extend into the apical dome where microtubules are relatively less bundled (S3B- S3E Fig). Microtubules converge on a region just behind the tip of the apex to form a microtubule focus. The distribution of MpWDL-YFP is characteristic microtubule-localized signal and resembles the localisation of GFP-MpTUB1 in growing wild type rhizoids (Fig 4C and 4D). MpWDL-YFP decorates the bundled microtubules along the shank and the microtubule network in the apical dome. This suggests that MpWDL and MpTUB1 co-localize. However, MpWDL-YFP fluorescence is stronger in microtubules along the shank than in the dome (Fig 4D). By contrast, GFP-MpTUB1 fluorescence is stronger in the apical dome and at the position of the microtubule focus than it is along the shank (Fig 4D). Consistent with this observation is the demonstration that the ratio between MpWDL-YFP and GFP-MpTUB1 is higher in the shank of growing rhizoids than in the apical dome ( Fig 4E). These data suggest that MpWDL-YFP localizes preferentially to microtubules along the shank in comparison to microtubules in the apical dome. The fact that MpWDL-YFP is more abundant in the region of the cell where microtubules are highly bundled than in regions where microtubules are less bundled suggests that MpWDL may be involved in microtubule bundling.

MpWDL-YFP promotes formation of a longitudinal array of parallelarranged bundled microtubules in the shank of growing rhizoids
Since the MpWDL microtubule-binding protein is more abundant in the rhizoid shank where microtubules are arranged in longitudinal, parallel bundles, we hypothesized that MpWDL was responsible for bundling the shank microtubules. To test this hypothesis, we compared the organization of microtubules decorated with MpGFP-MpTUB1 in growing wild type and wdl-2 rhizoids (Fig 5A and 5B). Cortical microtubules in the shank of wild type rhizoids are bundled (Fig 5A). Bundles are predominantly parallel. They run longitudinally from base to apex and converge at the tip where microtubule plus ends are located (Figs 5A and S3F and S3G). Cortical microtubules in the shank of growing wdl-2 rhizoids were significantly less parallel than in wild type and form a "net-like" array ( Fig 5B and 5D). Furthermore, cortical microtubules appeared less bundled in wld-2 rhizoids than wild type. To test the hypothesis that microtubules were less bundled in wld-2 mutants than wild type, we estimated relative bundling in mutant and wild type rhizoids. The relative bundling of microtubules can be estimated by the skewness in the fluorescence intensity distribution of microtubule pixels (see [19] for an application of skewness as an estimator of actin bundling). The skewness values were significantly lower (Fig 5B and 5C) in the shank of wdl-2 than in the shank of wild type growing rhizoids (Fig 5C and 5D). This suggests that cortical microtubules in the shank of growing rhizoids in wdl-2 were less bundled than in the shank growing rhizoids in wild type. These data are consistent with the hypothesis that MpWDL activity is required to bundle microtubules which promotes the formation of the parallel, longitudinal organization of the microtubule cytoskeleton in the shank of growing rhizoids. The wavy morphology and the less parallel, less bundled microtubles in the shank of wdl-2 loss of function mutants suggest that MpWDL promotes the organization of an array of bundled, parallel and longitudinally orientated microtubules in the shank that participate in stabilisizing the apex where the rhizoid extends by tip-growth.

Discussion
We demonstrate that proteins that are required for the formation of an array of bundled, parallel and longitudinally orientated microtubules in the non-growing shank stabilize the tip of rhizoid cells during growth. We also confirm that proteins that repress microtubule bundling in the apical dome of rhizoids stabilize tip growth. Taken together, our data indicate that MpWDL and MpNEK proteins act oppositely on microtubule bundling and both are required to stabilize tip-growth. We propose a model where tip-growth stability requires the microtubule bundling-promoting function of MpWDL and the microtubule destabilizing function of MpNEK to be spatially separated, thus maintaining a population of microtubule bundles in the non-growing shank of rhizoids and a population of dynamic microtubules growing toward the apical dome.
We report here that MpNEK is required to maintain the stable position of the growing point in tip-growing cells of Marchantia polymorpha. We isolated nek mutants from UV-B and T-DNA mutant screens and observed that an early stop codon in MpNEK causes rhizoids to grow wavy instead of straight. Our data are consistent with published data that demonstrate that nek loss of function mutant develop wavy rhizoids [16]. MpNEK had not been reported to control tip-growth stability in other land plant models. This may be because there are larger numbers of genes in the NEK gene family, with potentially redundant functions in tip growth, in the genomes of other land plant models. This means that loss of function mutations in these genes would not result in defective rhizoid morphology. Because MAP gene families are smaller in Marchantia polymorpha than in other land plant models, the probability of observing a defective phenotype is higher in single mutants than in other models like A. thaliana. Therefore, it is formally possible that identifying the causative mutations in other Marchantia polymorpha mutants with wavy rhizoids may reveal new regulators of microtubule-mediated tip-growth stability.
We report that MpWDL is a MAP that regulates tip-growth stability in M. polymorpha. WDL proteins form a land plant-specific MAP family that promotes microtubule bundling in vitro [18]. WDL homologs promote the longitudinal orientation of the cortical microtubule array and repress the anisotropic elongation of root epidermal cells and light-grown epidermal cells in A. thaliana [18,20,21]. However, they have not been shown to be required for the spatial stability of the site of tip-growth in A. thaliana although ectopic overexpression results in defective root hair development. In M. polymorpha, the loss of function wdl mutants form wavy rhizoids and microtubules are misaligned in the shank of growing rhizoids: while microtubules are oriented longitudinally in the shank of wild type rhizoids, they are more randomly oriented in wdl-2. This suggests that the function of WDL proteins in promoting the longitudinal orientation of microtubules may be conserved in cylinder-shaped cells of land plants that elongate by tip growth.
Taken together, our results suggest that tip-growth stability in M. polymorpha requires the function of a microtubule destabilizing MAP, MpNEK, and a MAP that promotes microtubule bundling, MpWDL. This is consistent with previous studies in A. thaliana where both the microtubule destabilizing protein ARK1 [9] and the microtubule bundling protein MOR1 [14] are required to stabilize tip-growth. Microtubules are more bundled in the root hairs of ark1 mutants in A. thaliana and the rhizoids of nek mutants in M. polymorpha compared to their respective wild type root hairs and rhizoids [9,16]. Therefore, MpWDL and MpNEK act oppositely on microtubule bundling and both are required to stabilize tip-growth.
We propose that tip-growth stability requires a spatial separation of MAP-mediated microtubule bundling and microtubule destabilization. MAP-mediated microtubule bundling may promote microtubule longitudinal orientation in the shank. In turn, longitudinally oriented microtubules grow toward the apical dome, where plus-end targeted microtubule destabilizing MAPs prevent excessive microtubule stability that would limit the pool of free tubulin available for the polymerisation of a self-renewing population of microtubules that target polarity markers to the apical plasma membrane.

Plant growth
Male and female M. polymorpha accessions Takaragaike-1 (Tak-1) and Takaragaike-2 (Tak-2) respectively provided by K. Ishizaki (Kobe University, Japan) were used as wild-type. T-DNA mutant lines were selected from a T-DNA mutant collection reported in [17]. Plants were grown horizontally and under continuous white light (~60 PPFD) on a Johnson's medium modified as detailed in [17]. To induce the formation of reproductive organs, one month-old gemmalings were transferred to a 1:3 mixture of medium vermiculite: John Innes No. 2 compost and grown under 150 PPFD white light supplemented with GreenPower LED research modules FR (8727900 809336 00, Phillips).

Drug treatments
Gemmae were placed on modified Johnson's medium [17] supplemented with oryzalin, Taxol or DMSO and mature rhizoids were imaged two days later. Oryzalin was dissolved in pure DMSO to produce 0.1 mM and 0.75 mM 1000X stock solutions. Taxol was dissolved in pure DMSO to produce 3.3 mM and 5 mM 1000X stock solutions. Treatment plates were compared to 0.1% DMSO control plates.
The open reading frame of MpWDL was amplified from genomic DNA using Phusion High-Fidelity DNA polymerase (NEB) and gene-specific primers (WDL_infusion _F1: ACCA ATTCAGTCGACATGAGTGAAGCCGGA and WDL_infusion_GS1_R1: GCCCTTGCTCA CCATGGATCCTGATGCAACAGCCA) designed to introduce a 15bp sequence homologous to HB39 and to replace the stop codon by a glycine-serine C-terminal linker. The YFP C-term entry vector HB39 was linearized using the restriction enzymes NcoI and SalI (NEB). The PCR product was gel purified and recombined with double digested HB39 using the In-Fusion kit (#638910, Clonetech). Finally, the MpWDL fusion entry vector was recombined by LR reaction (#12538-200, ThermoFisher) with the destination vector HB444 to produce a constitutive expression construct.
A 5 kb region upstream of the start codon of MpWDL was defined as the promoter region and amplified from genomic DNA using Phusion High-Fidelity DNA polymerase (NEB) and gene-specific primers (pWDL_F1: AAACGACGCCCAGTGCCACCCAATGCTTCAAAG TTTGAAG) designed to introduce a 15 bp sequence homologous to MPGWB207 and a SmaI restriction site downstream of the 5'UTR. Destination vector MPGWB207 was linearized using the restriction enzymes HindIII and SacI (NEB) to remove the pro MpEF1a cassette. The promoter amplicon was gel purified and recombined with double digested MPGWB207 using the In-Fusion kit (#638910, Clonetech). The resulting construct was linearized with the restriction enzyme SmaI and ligated with the Gateway Cassette frame A (Invitrogren) to produce a destination vector. The resulting pro MpWDL destination vector was finally recombined by LR reaction (#12538-200, ThermoFisher) with the MpWDL-GS-YFP entry vector to produce the pro MpWDL:MpWDL-GS-YFP transformation vector.
The expression clones were transformed into M. polymorpha spores obtained from a cross between Tak-1 and Tak-2 by co-cultivation with Agrobacterium tumefaciens GV3101 using a method from [17]. Co-cultivated sporelings were plated on modified Johnson's medium supplemented with 100 ppm cefotaxime plus 10 ppm hygromycin, 150 ppm gentamycin or 0,15 ppm chlorsulfuron depending on the constructs.

Reporter lines
Reporter lines for cytoskeleton and endomembrane compartments were generated by transformation of spores obtained from a cross between Tak-1 and Tak-2. Transformants were screened for adequate signal strength and absence of aberrant morphological phenotype. Selected reporter lines were then crossed to mutant backgrounds when required. Because the level of expression of certain reporter constructs increased or decreased after sexual reproduction, the F1 progeny with wavy rhizoids was compared with the F1 progeny with straight rhizoids, rather than with the parental reporter line.

Confocal imaging
To image growing rhizoids, one-day-old gemmalings were used. To image growth-arrested rhizoids, two-day-old gemmalings were used. Whether rhizoids were growing was tested by imaging rhizoids over the course of several minutes.
Imaging chambers were designed as described in [23], using cavity slides to allow rhizoids to grow without entering in contact with the cover slip. Gemmae were plated on the agar slice before adding perfluorodecalin and transferring to normal growth conditions. Imaging of rhizoids was finally carried out using a 20X or a 60X water immersion lens.
For the purpose of staining with fluorescent dyes, gemmae were grown on Petri dishes with 10 mL of modified Johnson's medium, rather than in imaging chambers. Gemmalings were immersed in 10 μM propidium iodide for 4 minutes and washed once by removing the staining solution and adding 20 mL ddH 2 O. Imaging was then conducted using a 40X water dipping lens.
All imaging was carried out using a TCS SP5 confocal microscope (Leica). Excitation and emission settings were determined empirically to minimize chloroplast and cell wall autofluorescence bleed-through. The same settings were used for GFP and Venus: excitation with 488 nm and emission between 500 nm and 570 nm. Brightness was adjusted in raw images to reach the point of saturation and the same modifications were applied to all samples to be compared.

Quantification of microtubule phenotype in growing rhizoids
We imaged 10 growing rhizoids from independent gemmae of either wild type or wdl-2 pro M-pEF1a:GFP-MpTUB1 originating from the same transformation event. The four cortical-most optical slices (395 nm) for 10 rhizoids from independent gemmae were maximum intensity projected and analysed using the LPX package for ImageJ, as described in [19]. A Mann-Whitney U test (or "two sample Wilcoxon" test) was performed to test the similarity between the wild type and wdl-2 skewness of the distribution of intensity signal in pixels corresponding to microtubules and the wild type and wdl-2 angles between lines of the skeletonized images.

Quantification of the ratio between MpWDL-YFP and MpTUB1-GFP signal
Fluorescence was imaged in growing rhizoids of 10 pro MpEF1a:GFP-MpTUB1 transformed wild type and 13 pro MpWDL:MpWDL-YFP transformed wild type. Signal intensity values were extracted along a 35 μm longitudinal line starting from the apical plasma membrane. The ratio between the average signal intensity values for the two reporter proteins was computed at each position to represent the relative MpWDL-YFP abundance.

Quantification of mature rhizoid morphological phenotype
Rhizoid sinuosity and rhizoid diameter were measured from the apical most region of growtharrested rhizoids that would fit in a 200 x 100 μm field of view. The quantification of rhizoid diameter was carried out by measuring manually the shortest diameter along the shank of rhizoids every 5 μm. The quantification of rhizoid 3D sinuosity was carried out using a python script. Software to analyse the waviness of rhizoids in 3D was developed in the Python programming language [24] and uses Python Imaging Library [25] (3.3.1), NumPy [26] (1.11.2) and scikit-image [27] (0.12. 3) to open and read images. Contours of cells and its inner pixels were segmented in each slice by applying a pixel intensity threshold with OpenCV (2.8.4) [28]. The noise of the image is removed by selecting contours larger than 1% of the total area of the image. The centres of gravity (CoGs) of these contours were calculated and used to sort the contours to represent the rhizoid in 3D. This order, referred to as the path, is determined by defining the shortest possible route between all central points of the contours using a Simulated Annealing algorithm [29] modified to deal with a three dimensional distance calculation for a non-circular path with a fixed starting point. After defining the path, the angle to the xaxis between each CoG and its successor was calculated to define the local direction of the path. The path is split when a significant change in direction occurs. The contours of the resulting fragments of the path are used for a non-parametric regression analysis using Locally Weighted Scatterplot Smoothing (LOWESS) [30], with the Statsmodels [31] (0.6.1) module. Parameterized with frac. set as 0.1. After the regression analysis of all parts, the regressions were assembled to one path. The pipeline from angle calculation to assembly of the regression was performed for both the x-y and x-z axis in order to maintain the 3D dataset. Finally, this assembly is used to calculate sinuosity as defined by the ratio between the total path length after regression by the shortest distance between base and tip of the rhizoid.

Expression studies
Three-week-old gemmalings were ground in liquid nitrogen using pestle and mortar. cDNA synthesis and qPCR were carried out as described in [32]. Absolute levels of expression were normalized using MpEF1a as a reference gene. Primers used for qPCR were as follows. MpE-F1a_forward: CCGAGATCCTGACCAAGG. MpEF1a_reverse: GAGGTGGGTACTCAGCG AAG. MpWDL_forward: GTTGCCTGTCCTCACGATCA. MpWDL_reverse: TCATGAC GCTTGGGCAGTAG.

Genomic DNA extraction
Genomic DNA was isolated from one-month-old plants grown in axenic conditions. Whole tissue was ground in liquid nitrogen using pestle and mortar. Genomic DNA was extracted with 2% CTAB buffer as described in [33]. To assess the suitability downstream applications, the quantity of extracted DNA was measured with a Qubit 2.0 Fluorometer (Invitrogen) and the purity was estimated from absorbance ratios measured with a ND-1000 spectrophotometer (Nanodrop).

Thermal asymmetric interlaced PCR
Thermal asymmetric interlaced PCR were carried out as detailed in [17].

UV-B mutagenesis
Fresh sporangia obtained from a cross between Tak-1 and Tak-2 were harvested and surface sterilized with 1% sodium dichloroisocyanurate (Sigma) solution for 4 minutes and subsequently washed three times in ddH 2 0 before plating on square plates with modified Johnson's medium at a density of 2000 spores per plate. Excess of water on the plates was allowed to dry out in the flow hood before opening the plates in a Benchtop 2UV Transilluminator (UVP) set at 302 nm and equipped with 6 G8T5E UV-B fluorescent lamps (Ushio). Spores were placed facing the UV-B light source with no plastic or medium in between and exposed for a duration of time sufficient to yield 50% kill, which was 60 s in our experimental conditions. Then, irradiated plates were closed and plates kept in the dark at room temperature overnight before being transferred to normal growth conditions. Finally, irradiated sporelings were assessed for survival and screened for mutant phenotypes 14 days after plating.

Non-allelism-based mutation discovery pipeline
The sequencing of UV-B mutants was carried out with Illumina's HiSeq-2000 platform by the Oxford Welcome Trust Center for Human Genomics. Raw reads were quality trimmed using Trimmomatic-0.32 [34] and normalized using Khmer-0.7.1 [35] with a k-mer size of 31. Resulting reads were aligned against the genome assembly using bowtie2-2.1.0 set in-verysensitive-local mode. Alignments were position sorted and mismatches within reads with q quality higher than 35 were extracted using the function sort and mpileup from bio-samtools-2.0.5. Mismatches in regions with coverage exceeding 100X were excluded using the varFilter function from bcftools of the samtools-0.1.9 package. Then, mismatches were retained only if they were supported by more than seven reads and if they appeared sufficiently homozygous based on a negative FQ value or AF1 value higher than 0.5001. Finally, retained mismatches were considered candidate SNPs if they were absent from other sequenced non-allelic lines, Tak-1 and Tak-2, if they were G2A or C2T substitutions or INDELs consistently with the expected UV-B mutation signature, and if they caused a change in the amino acid sequence of a predicted gene.

S1 Fig. Rhizoid-mediated soil adhesion of wild type (left) and wdl-3 mutant (right). Four
week-old gametophytes were transferred on soil and grown for 8 weeks before being pulled up by lifting the gametophyte and its rhizoid system in a vertical motion. More soil was attached to wild type than of wdl-3 suggesting that straight rhizoids anchor the gametophyte to the substratum more effectively than wavy rhizoids.