Computational Tools for Serial Block Electron Microscopy Reveal Plasmodesmata Distributions and Wall Environments

Plasmodesmata are small channels that connect plant cells. While recent technological advances have facilitated analysis of the ultrastructure of these channels, there are limitations to ef ﬁ ciently addressing their presence over an entire cellular interface. Here, we highlight the value of serial block electron microscopy for this purpose. We developed a computational pipeline to study plasmodesmata distributions and detect the presence/absence of plasmodesmata clusters, or pit ﬁ elds, at the phloem unloading interfaces of Arabidopsis ( Arabidopsis thaliana ) roots. Pit ﬁ elds were visualized and quanti ﬁ ed. As the wall environment of plasmodesmata is highly specialized, we also designed a tool to extract the thickness of the extracellular matrix at and outside of plasmodesmata positions. We detected and quanti ﬁ ed clear wall thinning around plasmodesmata with differences between genotypes, including the recently published plm-2 sphingolipid mutant. Our tools open avenues for quantitative approaches in the analysis of symplastic traf ﬁ cking.

The cellular units of complex organisms have an intrinsic need for communication. In plants, effective signal exchange is enabled by plasmodesmata (PD), small channels connecting neighboring plant cells (for review, see Nicolas et al., 2017a). While research has largely focused on the structure and biological regulation of the aperture of the PD, recent insights point to the importance of PD spatial arrangements and their cell wall environments for the flow of materials through them (Deinum et al., 2019). Novel methods to get comprehensive information on PD are therefore now required.
The nanometer size of PD pores poses a challenge for their study. A tradeoff exists between resolving the detailed structure of the channels and capturing their overall distribution. Electron microscopy (EM) resolved the structure of these channels, identifying a continuous plasma membrane and a constricted form of the endoplasmic reticulum (ER), the desmotubule, running across the pore between the two cells (Lopez-Saez et al., 1966;Robards, 1971). More recently, with the application of electron tomography, variable apposition of plasma membrane-ER membranes was shown (Nicolas et al., 2017b;Yan et al., 2019). Classical EM can also be used to study PD occurrence. An inventory of PD densities along the Arabidopsis (Arabidopsis thaliana) root highlighted interesting variation between cellular interfaces, which might underpin qualitative or quantitative differences in PD-mediated communication between cells (Zhu et al., 1998). However, EM approaches, when looking at single or separate slices, largely lose information about the positions of PD relative to each other and only capture approximate densities. This is problematic because the distribution of PD across an interface is predicted to have a significant impact on flow properties (Deinum et al., 2019). Limited alternatives to comprehensively address the presence of PD have since emerged. Confocal microscopy was applied in leaves, using specific PD markers, to show that the development and distribution of particular PD morphologies in the epidermis were strongly increased by treatments eliciting nutrient, osmotic, and pathogen stresses (Fitzgibbon et al., 2013). Fluorescent approaches are, however, limited to relatively accessible cell-cell interfaces and often cannot resolve the signal from individual PD. Faulkner et al. (2008) used freeze-fractured trichomes and EM to analyze PD distributions across the entire fractured surface. They observed that new PD (not generated during cell division) seemed to insert themselves in close proximity to existing PD, suggesting the use of the latter as nucleation centers. The process ultimately results in clusters of PD pit fields. A method to obtain similar interface-level estimates of PD densities, in this case in the mesophyll layer of leaves, was introduced by Danila et al., (2016), combining 3D immunolocalization, to determine the area of the pit field relative to that of the interface, and scanning EM, to assess the number of PD per pit field. However, both immunochemistry-scanning EM and freeze-fracture approaches remain confined to tissues that are readily accessible to such sample processing.
Serial block face electron microscopy (SB-EM; Denk and Horstmann, 2004) can overcome these limitations, offering the opportunity to look at interfaces deep in tissues. A block of fixed and embedded tissue is mounted inside a scanning EM device and the upper face of the block is cut away using an internal microtome. After each slice, the newly exposed block surface is imaged. The process is repeated, ultimately generating a stack of images along a z-axis with the z-resolution defined by the thickness of the slices. Importantly, the positions of cellular objects are retained relative to one another, and the data sets are good starting points for 3D reconstruction (for review, see Kittelmann et al., 2016). SB-EM technology has been successfully employed to study PD, demonstrating defects in sieve pore (a modified form of PD) structure and distribution (Dettmer et al., 2014) and allowing the quantification of PD densities at the interfaces of the sieve element (SE; Ross-Elliott et al., 2017) and at the endodermal (EN) face of phloem pole pericycle (PPP) cells (Yan et al., 2019). Both SE and PPP cells are key players in the process of phloem unloading, largely mediated by PD (for review, see Truernit, 2017). These data sets are, however, underexploited in part due to limitations in the technology to extract such information from them. Consequently, important parameters, such as the specific distributions of PD and the cell wall environment of the pores, despite being contained in these data sets, have so far been ignored.
Here, we address these two biological aspects. Dense clustering of PD into pit fields is often assumed as a general feature of these structures (Sager and Lee, 2014). However, while this is certainly the case at some interfaces (Faulkner et al., 2008;Danila et al., 2016), additional evidence is needed to support a generalization. Recent modeling efforts have highlighted how the arrangement of PD in clusters might actually reduce flow between cells (compared with a random arrangement; Deinum et al., 2019). Having detailed information on distributions in actual cells would greatly inform these models. The local wall environment in which PD reside is also of relevance for flow. The thickness of the wall at PD defines the length the path substances have to travel before entering the neighboring cell. Thinning at PD is often assumed, but the evidence is not comprehensive and quantifications are not available.
Correlations between wall thicknesses and different PD ultrastructures have been reported (Nicolas et al., 2017b), and this is now being integrated into models, with predicted effects on flow (Deinum et al., 2019). We also know that the PD environment is peculiar in terms of wall polysaccharides, with an enrichment in callose and pectins and a concomitant reduction in cellulose (for review, see Knox and Benitez-Alfonso, 2014). Overall, the properties generated by wall components have not been extensively explored in planta, partly due to the difficulty of efficiently imaging phenotypic effects.
To extract the relevant information from SB-EM data sets, we developed novel computational and visualization tools dedicated to PD analysis. We deployed the SB-EM data sets from Yan et al. (2019) as a study case. We first address the spatial distribution of PD. We detected clusters of PD at the PPP-EN interface, while we did not see signs of clustering at the SE-PPP interface. We quantified the number and size of the clusters. We quantified specific wall thinning at PD positions and detected changes in the wall environment in the plm-2 Arabidopsis mutant.

SB-EM Allows Spatial Positioning of PD over Wall Interfaces
SB-EM data sets can cover large portions of a tissue. The data sets from Yan et al. (2019) employed here cover an area encompassing the cells around the protophloem of Arabidopsis roots. The data sets can be visualized either in a longitudinal orientation (Fig. 1A) or in an axial one (Fig. 1B). In the latter, PD at various radial interfaces are more easily detectable due to better xy resolution of the SB-EM technique (Fig. 1,C and D,respectively). The image resolution of the specific data sets employed here is good enough to identify individual PD within those areas but not to distinguish the detailed morphology of the PD. A unique aspect of SB-EM is that such annotated PD positions (relative to the cell surface) can be addressed globally within the full length of cells. While density calculation approaches have taken advantage of this (Yan et al., 2019), the spatial component, namely the 3D distribution of PD, has been largely neglected. Traditional bidimensional visualizations fail to convey the distribution of PD over the interfaces. Here, we show that identified PD can be exported (as clouds of dots) alongside the segmented wall, generating effective 3D spatial representations that capture the distribution. We show this both at the SE-PPP and PPP-EN interfaces (Fig. 1, E and F). The rendering can also be stored as movies (Supplemental Video S1).

Rendering of PD in the Cellular Context
SB-EM data sets also contain data on various organelles within the cells, putting PD in a wider and more realistic context. We can generate highly structured and dense cellular models. By segmenting an area of the data sets (represented in the overall data set in Fig. 1H, inset) and color coding the different organelles ( Fig. 1G), the 3D model can be eventually exported in visualization programs (Fig. 1H). We can show the ER strands of PD crossing into the wall (Fig. 1I) and then merging on either side to the wider ER network system. Various animations can then be applied to the 3D model (Supplemental Video S2). Models such as these highlight how symplastic transport needs to navigate a dense cytoplasm before reaching the PD.

Signs of Clustering of PD at the PPP-EN Interface in the Root
Taking advantage of the spatial information on PD contained in the SB-EM data sets, we studied their distribution at selected interfaces. To describe the distribution of PD on the cell wall, we calculated pairwise Euclidean distances between each of them, using the x, y, and z coordinates available in the data sets. This revealed a multimodal distribution of distances (red lines in Fig. 2A, where we show two examples of cells). To provide a meaningful comparison, an equal number of points with a uniform distribution was generated on the same surface using the SpatialControlPoints tool described in "Materials and Methods." One thousand simulations were generated for each cell. We calculated Euclidean distances for each of the simulations and observed that in each case the distributions of distances approximated a normal distributions (yellow lines in Fig. 2A, for the two example cells). The surfaces are not perfectly flat, resulting in deviations from full normality. This immediately suggests some overall differences compared with the distances between real points. Next to each plot for the distribution of distances we also show the original 3D distribution of points, to give an  . The red bar highlights the 0 value, representing identity between real and simulated distributions. The left dark gray curve represents data from the top cell in A, and the right dark gray curve represents data from the bottom cell in A. E, Comparison of KS test values at the appreciation of the biases in point distribution between real data and simulation (Fig. 2B, showing the same two example cells). At the PPP-EN interface, in addition, an excess of short distances between PD points was always visually detectable compared with the distributions of distances of the simulated points. This suggests some form of clustering ( Fig. 2A). Note that PD present in kinked areas of the wall (Fig. 2B, PPP to the right), which represent a minority of cases, would result in some diagonal Euclidean distances outside of the wall surface. However, simulated points on that same wall would experience similar effects, making them comparable. For each of the eight cells tested for Col-0 or the five cells tested for plm-2 (the mutant available in the data sets of Yan et al. [2019]), pairwise Kolmogorov-Smirnov (KS) tests of the distribution of real points against each one of the 1,000 simulated point sets were performed. The distribution of P values is shown in Figure 2C (in pink for the PPP-EN interface). All cells fell below a P value of 0.05 (black vertical line), suggesting a nonuniform distribution of PD. To give a quantitative appreciation of variation across cells in spatial distributions, the distribution of the KS test result values can be plotted (Fig. 2D). Overall, test values for the single comparisons ranged from 0.025 to 0.347. Figure 2, A and B, actually display the two most extreme cases among cells of the Col-0 genotype: we took the cell with the distribution most shifted to lower KS values (PPP1-ENa) and that shifted to higher values (PPP2a-EN), and these are shaded in a darker gray in the figure. The spatial plots in Figure 2B match well with the expectations. Comparisons between genotypes can be performed by plotting mean KS test values for each cell of both genotypes (Fig. 2E). Summary values per cell remove the otherwise present problem of interdependencies of points, which would complicate statistical comparisons. The plm-2 genotype did not show appreciable shifts compared with the wild type (medians of 0.09 and 0.08, respectively).

Lack of PD Clustering at the SE-PPP Interface in the Root
At the SE-PPP interface in Col-0, conversely, there was no evidence to reject the null hypothesis of a uniform distribution for the PD. The P values for the KS test comparison between real and simulated points were shifted toward or above 0.05 in Figure 1C (light blue curves). Mean P values for each cell were above 0.05. We show the distribution of Euclidean distances for one of the cells. The distribution of distances for the real points appeared less diverse relative to the simulated points, compared with those observed at the PPP-EN interface. In addition, no visual excesses of shorter distances could be detected (Supplemental Fig. 1A). The 3D distributions of real and simulated points are shown in Supplemental Figure 1B. The distributions of the KS test values, while being at times higher than those at the PPP-EN interface in terms of absolute values, were much shallower (Supplemental Fig. 1C; the dark shaded cell in this image was used as the example in Supplemental Fig. 1, A and B).
Because the SE-PPP interface has a lower number of PD compared with the PPP-EN interface (Supplemental Fig. 1D), we tried to assess if the high P values at the SE-PPP interface were just due to lower statistical test power or were an indication of actual lack of clustering. We sampled a lower number of PD (and simulated points) at the PPP-EN to achieve the same PD density as that seen at the SE-PPP interface. The number of new points was calculated by multiplying the density of PD at the SE-PPP interface by the surface area at the PPP-EN interface (Supplemental Fig. 1D). We then tested if a difference between the distribution of Euclidean distances of real points and simulated ones could still be detected. While the P values did indeed on average shift toward higher values, in six of eight cells of Col-0 the mean P value was still below 0.05 (red vertical line; Supplemental Fig. 1E). In only two cells (purple ones in the figure), the PD distribution could no longer be robustly differentiated from a uniform one. Overall, this suggests that at the SE-PPP interface there are indeed no obvious signs of PD clustering, and this highlights differences between this interface and the PPP-EN interface.

Describing the Organization of PD in Pit Fields at the PPP-EN Interface
Upon establishing the presence of a nonuniform distribution of PD over the PPP-EN interface, we attempted to characterize the potential clusters. Namely, we tried to address the number of clusters, the number of PD per cluster, and the cluster sizes relative to the surface of the interface. To determine the number of clusters, we used two different clustering algorithms, a k-means-based method and a model based one, within the R environment. Variation was visible, as should be expected due to the relatively arbitrary computational classification, between the single cells and between clustering algorithms with median values of 11.5 (Col-0) and 10 (plm-2) using the mclust package and with median values of 11.5 (Col-0) and 14 (plm-2) in the silhouette approach (Fig. 3A). Differences between genotypes were not statistically significant, so, overall, a working range of 10 to 14 PD clusters can be suggested at the PPP-EN interface. As an example, we color-coded the PD of a cell according to the cluster they had been assigned with two methods (Fig. 3, B and C). In the image, the 3D coordinates had been reduced to 2D via PCA. Some of the strengths and pitfalls of these clustering methods are illustrated in this example, with the silhouette approach being possibly overly conservative while mclust assigned one PD to the wrong cluster (the olive green dot in the bright green cluster in Fig. 3B, highlighted with an asterisk). We strongly emphasize that cluster number values should be used as working ranges rather than absolute values. The number of PD in each cluster was similar between the two clustering methods, with a median of eight to 10 PD per cluster (Fig. 3D). Once again, no strong trends in the plm-2 mutant from the datasets of Yan et al. (2019) were detectable. Figure 3. Quantification of clustering parameters at the PPP-EN interface using SB-EM data sets. A, Number of clusters identified in the Col-0 (eight cells) and plm-2 (five cells) genotypes using the mclust or the silhouette approach. B and C, Visualizations of a principal component analysis (PCA) reduced interface (from a cell) with different cluster assignments. The surface is rendered in gray, while PD belonging to different clusters and the area they occupy are color-coded. D, Number of PD per cluster in Col-0 and plm-2 genotypes. The total clusters are 95 (mclust) and 96 (silhouette) for Col-0 and 52 (mclust) 55 (silhouette) for plm-2. E, Total percentage of the surface occupied by clusters. F, Percentage of the surface occupied by individual clusters. G, Absolute surface in micrometers squared of cells. H and I, 3D visualizations in the Imaris program of the same interface shown in B and C. The surface is rendered in gray, while PD belonging to different clusters are color-coded with the same scheme used previously. Note that D and F have logarithmic y axes. In the graphs, individual values are represented as dots, distributions as violin plots, and medians as horizontal bars. Cells from different roots are shown using different symbols. For each clustering approach, statistical comparisons between genotypes were made using the nonparametric Mann-Whitney U test for two samples; ns 5 P . 0.05.
To assess the conductive surfaces provided by these PD clusters, we calculated the area occupied by each of these clusters and that occupied in total by all the clusters on a cell interface. These areas were calculated and reported as percentages of the total surface of the interface on which they occur, as delimited by the most extreme cluster points in the PCA space (shaded gray area in Fig. 3, B and C). While this is an underestimate of the total surface, we feel it is more appropriate to calculate these scaling factors than attempting to use the cluster surfaces as absolute values (the data are indeed scaled and reduced by the PCA, so the units are no longer true mm 2 ). For the total conductive surface (i.e. the proportion occupied by the clusters relative to the overall surface), the mclust method suggests medians of 14% for both Col-0 and plm-2 while the silhouette approach provides median values of 15% and 18%, respectively (Fig. 3E). Each cluster accounts for a median surface of 1% in both genotypes using the mclust method or 0.7% (Col-0) and 0.8% (plm-2) with the silhouette approach (Fig. 3F). No differences between the two genotypes were highlighted by statistical testing for any of these parameters.
While these might be sufficient for some purposes, we also wanted the possibility to relate these percentages to the actual surface values in mm 2 . To do so, we had to employ the original image data within the MIB software, rather than the R-processed and PCAreduced ones. Using an updated version of the Surfa-ceArea3D plugin employed by Yan et al. (2019), we calculated the actual total surface of the PPP-EN interfaces in MIB (for details, see "Materials and Methods"). Median surface areas of 91.1 mm 2 for Col-0 and 119 mm 2 for plm-2 were determined (Fig. 3G). Given the variance in the data, this difference is not robust. Relating the median surface area values to the scaling factors described above, we can obtain the actual surfaces of the individual clusters. The surface can be exported to Imaris for visualization (Fig. 3, H and I, mirroring Fig. 3, B and C).

Extracting and Visualizing Wall Thickness at PD, Controls, and Every (Other) Position
The SB-EM data sets contain information on many components of a cell, and in the context of PD, a relevant one is the cell wall and its thickness. To do that, we developed the CellWallThickness plugin, which can extract the thickness of a given segmented wall at positions of interest (see "Materials and Methods"). As an example, we employ the plugin on one of the Col-0 cells available (the one used in Fig. 3, B and C). By using an equal number of random uniformly distributed points (median of 117 nm), one can accurately capture the thickness of the overall all-other wall (median of 123 nm) in a computationally effective manner. These values are not statistically significantly different. The data also show a clearly thinner wall in the proximity of PD (median of 46 nm; Fig. 4A). The thickness of the wall at the interface of interest can also be visualized graphically, by exporting the midline thickness map (generated by the plugin) into 3D rendering software such as Imaris (Fig. 4B). The wall color intensity matches the calculated thickness value at that position, brighter meaning thicker. The PD positions and those of the controls can also be exported as dots and their relative size made to match the wall thickness value. The thinning at PD positions is visually confirmed and shown to extend beyond the precise position of the channels, to the entire pit field PD are grouped into.
Cell wall thickness comparisons between genotypes are also possible with our tools, using mean thickness values for the three categories of points in each cell. At PD positions, the median thickness of the resulting values was 53 nm for Col-0 and 62 nm for plm-2 (Fig. 4C). The difference is supported by statistical testing. Conversely, no difference was supported for the random uniformly distributed points (132 nm in Col-0 versus 141 nm in plm-2) or for the all-otherpoints category (131 nm in Col-0 versus 142 nm in plm-2), despite a trend for increased thickness in plm-2 (Fig. 4C).
To test if the difference at PD positions between genotypes could be independently confirmed, we looked at another data set described by Yan et al. (2019). Electron tomography, a different technique, had been employed to study the ultrastructure of single PD at the PPP-EN interface in the two genotypes. When we redeployed those data, this time extracting the wall thicknesses in the immediate proximity of PD, we obtained median values of 66 nm in Col-0 and 85 nm for the plm-2 mutant (Fig. 4D). Statistical testing supported a difference. These values are remarkably close to those obtained with our plugin: an absolute match was unlikely considering the difference in scale of observation (individual PD compared with the entire tissue). Both techniques therefore agree in showing a trend of thicker walls at PD (and possibly across the entire wall) in the plm-2 mutant.
To further validate the reliability of the data generated, we assessed if known biological features could be detected in our data sets and if the values obtained matched those from different techniques. The wall of enucleated SEs is thicker compared with that of nucleated SEs or of surrounding cells. This reinforcement is likely necessary to withstand the pressure of sap flow (Furuta et al., 2014). The plugin output was able to effectively capture and quantify this difference using the all-points category. The median thickness of the SE-PPP wall, using averages per single cells, was 207 nm compared with 131 nm for the PPP (Fig. 4E). Note that here all points are used rather than all other points, as there is no need to exclude PD positions.

DISCUSSION
PD perform a key role in cell-cell transport across plant cells. We developed new computational tools to explore aspects of their distributions and of the wall they span. We performed this in part with the aim of informing future models of flow across PD with relevant experimental data. Published modeling approaches have so far studied PD transport in relation to the overall single-pore structure (Blake, 1978), phloem flow (Jensen et al., 2012), phloem-loading mechanisms (Comtet et al., 2017), and unloading flow type (Ross-Elliott et al., 2017). One modeling study tried to address some complexities of PD, integrating ultrastructure parameters for the cytoplasmic sleeve in their models (Liesche and Schulz, 2013).
The assumptions of that article, however, have been challenged by experimental data (Ding et al., 1992;Nicolas et al., 2017b), highlighting the difficulty of modeling flow across PD when limited experimental data are available. Only recently, the spatial distribution of the PD at interfaces, namely the assumption of clustering in pit fields, is starting to be included in models. Unequal distribution was shown to reduce the effective symplastic permeability of the interface (Deinum et al., 2019). However, detailed experimental data for distribution parameters, such as those we  Figure 3. The surface is rendered in shades of blue, on a scale matching the thickness of the wall. PD are shown as colored dots (red for real PD positions and yellow for those of a simulation with a random uniform distribution). The size of the dots relates to the wall thickness at that position. C, Comparison of wall thickness in Col-0 (eight cells) and plm-2 (five cells) genotypes. D, Comparison of wall thickness at PD positions in Col-0 and plm-2 genotypes using tomography data (n 5 30 for Col-0 and n 5 49 for plm-2). E, Comparison of overall wall thickness at the SE-PPP and PPP-EN interfaces in Col-0 (eight cells). In the graphs, average values for single cells are represented with symbols (different symbols for different roots), distributions as violin plots, and medians as horizontal bars. Note that A, C, and D have logarithmic y axes. Statistical comparisons between two samples (or two samples within a category) were performed using the nonparametric Mann-Whitney U test, and those for more than two samples were performed using the nonparametric Dunn's test. Supported differences are highlighted by an asterisks (*P , 0.05); ns 5 P . 0.05. present here, are lacking. Including PD spatial arrangements in models is a significant advancement from the use of the same spatial parameters (number of clusters and PD per cluster) to estimate the total conductive surface alone. Using traditional microscopy approaches, Kuo et al. (1974) had calculated assimilate fluxes in wheat (Triticum aestivum) leaves, and, more recently and more comprehensively, Danila et al. (2016) had used immunolabeling in accessible leaf tissues to compare fluxes in C3 and C4 monocot leaves. However, such studies might have underestimated the impediment to flow such PD arrangements may impose.
Using our first new plugin and its analysis pipeline, we show that at an interface important for postphloem unloading in the Arabidopsis root, that of PPP-EN (Yan et al., 2019), there is clear spatial clustering of PD (Fig. 2, A-C). We determined and visualized the numbers of pit fields at the PPP-EN interface, with a median of 10 to 14 clusters (Fig. 3A). Computer-driven clustering methods present their own limitations in terms of absolute value generation, possibly explaining part of the observed variation within and between clustering approaches. However, pit field numbers in different cells might also be affected by the surface and age of the actual cell. This in turn would explain the clustering metric variation observed in Figure 2A. We do not have enough cells spanning the z direction to empirically test such a hypothesis, but it remains plausible and interesting for future research. Overall, we recommend using these values as working ranges rather than absolute values. The remarkably stable median value for the number of PD per pit field, around 10 with all methods (Fig. 3D), might conversely hint at some biological process eventually constraining cells to create new clusters as they expand. There might be an upper boundary for secondary PD formation within a cluster. Lastly, in cells, the orientation of the PD and the overarching clusters (Figs. 1F and 3, H and I) often seem to follow the direction of vertical cell elongation.
We also calculated the median surface area of the clusters and related it to the total surface (Fig. 3, E and G). An approximate median conductive surface of 15% (each cluster accounting for a median of around 1%) is a remarkably low percentage and certainly suggests the possibility that flow in and out of a cell might not be uniform but rather resembles more a series of water currents in a cytoplasmic ocean. It will be extremely interesting to apply such a concept to modeling studies of flow between cells, as started in Deinum et al. (2019). Percentages can be related to actual surface areas using our second plugin.
In our data, we do not observe clustering at the SE-PPP interface ( Fig. 2C; Supplemental Fig. 1, B and C), which is fundamental for phloem unloading in the root (Ross-Elliott et al., 2017). A more uniform distribution of PD might reflect the enucleated nature of SEs (Furuta et al., 2014) and the impact this might have on secondary PD formation or a unique feature of the funnel PD at this interface, known to perform batch unloading (Ross-Elliott et al., 2017). It is interesting that in the study of Deinum et al. (2019), clustering of PD had the largest negative effect on parameters regulating flow at lower PD densities. The lower density of PD at the SE-PPP interface compared with neighboring tissues (Yan et al., 2019) might impose limits to clustering if this was to compromise the extensive flow that needs to take place at this interface. Overall, this result at least challenges the broad assumption that all PD might be grouped into pit fields. The distribution tools presented in this article could be used for more systematic studies within tissues.
In addition to the spatial distribution of PD, parameters describing the environment surrounding PD can also be of high value. We address wall thicknesses, affecting flow between cells in relation to the structure of PD in recent modeling studies (Deinum et al., 2019), with our third tool. The overall wall thickness is in the same order of magnitude as the estimation by Kramer et al. (2007). The value of around 200 6 30 nm was actually derived from a figure published by Andème-Onzighi et al. (2002), focused on root epidermal cells. We detect a thinner wall around PD at the PPP-EN interface in the root, by a factor of about 2.5 times (Fig. 4C), matching assumptions in the literature that PD lie in wall depressions. The agreement of wall thickness values at PD between SB-EM and electron tomography (a technique that focuses on the area of one PD) is extremely satisfactory, although we still caution on using these values as absolute. Thinner walls at PD might be the consequence of cell wall modifications required for PD de novo insertion (Faulkner et al., 2008) or a prerequisite for PD insertion at all. Regardless of the ontological reason of this wall thinning (or lack of thickening), it is likely that it carries functional relevance for conductivity (Thompson and Holbrook, 2003;Barratt et al., 2011). A few reports from other plant species mention that the sieve pores (highly modified forms of PD) in mature plates connecting SEs lie in wall depressions. This was correlated to callose deposition inhibiting wall thickening (Evert et al., 1966;Deshpande, 1974Deshpande, , 1975. Whether that is a shared mechanism to all PD, also rich in callose, is unknown. In modeling studies, a relative arbitrary value of 100 nm is employed as the wall thickness for PD (Liesche and Schulz, 2013). This is compatible with the averages for our PPP-EN cells (Fig. 4, C and D). However, our work flow might be highly valuable in future studies to inform cell-cell permeability models of differences between interfaces. For instance, we quantify the known biological difference of SE wall thickness (Furuta et al., 2014) compared with PPP cells (Fig. 4E).
Our tool could easily be applied to broad cell wall questions. For instance, to our knowledge, in the literature there are no tissue-specific studies of cell wall thickness in the Arabidopsis root. In addition, although SB-EM is not yet high throughput enough to allow mutant screens, targeted validation of mutants can be performed. We detected thicker walls around PD (and possibly globally) in the plm-2 mutant compared with the wild type (Fig. 4, C and D). The mutant is defective in the biosynthesis of very-long-chain fatty acid-containing sphingolipids (Yan et al., 2019). Glycosylinositol phosphorylceramide sphingolipids, known to be enriched at PD (Grison et al., 2015), are cross-linked via boron bridges with pectins (Voxeur and Fry, 2014), also likely enriched at PD (for review, see Knox and Benitez-Alfonso, 2014). It is therefore reasonable to speculate that there may be feedback effects on wall structure from lipid perturbations. The detected thicker wall emphasizes the importance of PD type change observed in this mutant for flow: it must provide a significant ease of trafficking to achieve the reported increase in communication (Yan et al., 2019). Alternatively, modeling studies suggest that different types of PD might be more suitable in different types of walls (Deinum et al., 2019).
Overall, interesting new lines of research might develop from a more systematic use of SB-EM and the associated analysis tools we provided in this article. We also envisage that future developments of the plugins will be able to expand the quantitative efforts to more aspects of the data sets.

Data Sets
For details on the equipment and settings for SB-EM image acquisition, we refer to the original article by Yan et al. (2019), for which the data sets were generated. Briefly, chemically fixed roots of 5-d-old Arabidopsis (Arabidopsis thaliana) plants were sectioned and imaged with cutting steps of 40 nm and xy resolutions of 7 to 10 nm. The collected images were assembled into a single calibrated, aligned, and contrast-normalized image stack. The resulting data sets are available at EMPIAR, the Electron Microscopy Public Image Archive (Iudin et al., 2016), with the accession code EMPIAR-10442. Images were loaded into the MIB software (Supplemental Fig. 2A; Belevich et al., 2016), downloadable at http://mib.helsinki.fi/ (last accessed March 2020), and they were filtered to reduce noise using deep neural network algorithms, which preserve the edges of the organelles (Supplemental Fig. 2B; Zhang et al., 2017). Tutorials on how to operate the software and its tools are available on the Web site. For the analysis presented here, we trimmed the original data sets, each relating to a root and containing multiple cells, into separate data sets for each cell (eight for Col-0 and five for plm-2). One of the original data sets for the plm-2 mutant had to be discarded, as the image quality was not sufficient for the specific purposes of this study, reducing the available cell number. Removal of the data set was performed before data analysis.

Annotations of PD and Cell Wall Segmentations
For PD annotations, we redeployed those from Yan et al. (2019). These are contained in the annotation layer of MIB, each annotation including the x, y, and z coordinates of the corresponding PD. For analysis, all coordinates were recalculated from pixels to physical units of the data set (micrometers) relative to the bounding box of each data set. In order to avoid duplicate counts, no new PD were annotated within 160 nm (12 or 22 slices from a central one) of an existing annotation, which provides a conservative estimate (Supplemental Fig. 2B). Within the MIB environment, we fully segmented the cell walls of interest by employing black-and-white thresholding within preselected masked areas (Supplemental Fig. 2C). Selection of the masked area encapsulating the cell wall was done using the brush tool and an interpolation process to infer the drawn areas on intermediate slices. Resulting models were smoothed and filtered so that the cell wall formed one continuous object in the 3D space. The final model was manually checked for any possible impurities. Small (less than five pixels in size) 2D profiles within the 3D model that might not be reliable were removed. High-quality segmentations and careful annotations of PD are the basis of any analysis employing the plugins we describe in this article.
They are the most time-consuming components of the pipeline, as they involve manual work from the user.

Plugins
All the computational tools employed in this article were written in Matlab language, but they do not require this proprietary software to operate. They are implemented as plugins for the freely available MIB software (Belevich et al., 2016). The plugins are included in the stand-alone version of the MIB software at http://mib.helsinki.fi/downloads.html (last accessed March 2020) or separately from https://github.com/AndreaPaterlini/Plasmodesmata_dist_wall (last accessed March 2020). In the latter case, they need to be saved in the Plugin folder of the MIB software. The plugins are provided with help sections. The research community can improve or tweak these plugins, according to their specific needs, by editing the source code. An overview of the type of file outputs generated by the plugins is provided in Supplemental Figure 2. The plugins require as inputs the initial manual steps described in previous sections.

The SpatialControlPoints Plugin
In order to ask questions relating to the distribution of PD, we felt that a comparable (in terms of points) simulated distribution had to be generated. The simulated distribution differs from the real PD one in that the points are placed with a spatially uniform pattern. To generate such a distribution, we developed a computational tool, the SpatialControlPoints plugin, capable of creating the control point distributions over the same surface as those present in the SB-EM data sets. A segmented wall and a list of annotated PD are fed into the plugin. The tool, in return, finds the midline of the wall by thinning the model to a single centerline without branches (Supplemental Fig. 2D). The thinning morphological operation (Lam et al., 1992) is applied to each slice, and then a function detects the longest available path within the thinned lines and removes all the others. The resulting single thin centerline is placed in the mask layer of the MIB interface. On this centerline surface, this tool generates an equal number of points to that of the PD, whose positions are sampled from a uniform distribution, with a randomly placed starting point. For reproducibility of results, in the user interface we provide an option to specify the random seed used by the sampling algorithm. The number of simulated distributions can be defined in the user interface; here, we employed 1,000 simulations. Matlab and csv file formats are available as outputs (Supplemental Fig. 2E).

The SurfaceArea3D Plugin
To calculate the surface of interfaces of interest in the SB-EM data sets, we employed an edited and improved version of the plugin used by Yan et al. (2019) for the same purpose. The plugin finds the midline of a supplied segmented wall on each slice of the model. This step is the same as that described for the SpatialControlPoints plugin (Supplemental Fig. 2D). This plugin then, additionally, connects such midlines across the slices, generating a surface (Supplemental Fig. 2F). The plugin employs the triangulateCurvePair function from the MatGeom toolbox for geometric computing with Matlab (https:// github.com/mattools/matGeom; last accessed March 2020). Matlab, Excel, and csv file formats are available for the numerical output of the surface. The surface itself can be exported as an object to Matlab, Amira, and Imaris programs.

The CellWallThickness Plugin
In order to explore the environment surrounding PD, namely the cell wall they span, we developed the CellWallThickness plugin to extract wall thickness from SB-EM data sets (Supplemental Fig. 2G). The plugin is fed a segmented wall and finds its centerline, as described for the SpatialControlPoints plugin (Supplemental Fig. 2D). A distance map, which assigns a value to each model pixel based on its distance to the closest edge of the model, is calculated at each slice using the Euclidean distance transformation algorithm (Maurer et al., 2003;Supplemental Fig. 2H). The values at each point of the masked centerline are then obtained by placing the centerline over the distance map image (Supplemental Fig. 2I). The values are expressed in pixels. Since the image is calibrated, the plugin then recomputes those numbers to actual physical thickness of the wall as thickness (in mm) 5 value (in pixels) 3 pixel size 3 2. The doubling factor is introduced to obtain wall thickness (and not just half thickness). The masked centerline, where each pixel encodes the rounded thickness of the cell wall at the corresponding point, can also be saved as an image file. Employing the annotated PD positions, the plugin looks for the closest position on the midline. A line over the wall to show one-half of the distance is displayed (Supplemental Fig. 2J). In addition to PD position, if requested, it generates a random uniform distribution over the same surface (employing the SpatialControlPoints plugin; Supplemental Fig. 2E) and samples an equal number of points to those of the PD. It can also extract the wall thickness at all points. Depending on the task, the values of real PD and randomly placed PD can be excluded from the list of all points using the corresponding option checkboxes. This ensures the independence of classes for statistical comparisons. Matlab, Excel, and csv file formats are available as outputs (Supplemental Fig. 2G).

R Scripts and Guided Pipeline Availability
The data obtained from MIB and its plugins were then imported and analyzed in R (R Core Team, 2017) to obtain a range of descriptive statistics. We stress that the data-analysis pipeline we employed here is one of many possible ones (from the same data outputs of the plugins). For instance, we calculated pairwise Euclidean distances between PD (or simulated points) to describe PD distributions. The approach was chosen because it is independent of the surface that PD sit in (and its boundaries). It rather focuses on the relationship between the individual points alone. Alternative spatial analyses, such as Ripley's K function, are also in principle possible. However, they will require specific implementations. The scatterplot3d package (Ligges and Mächler, 2003;version 0.3-41) was used to visualize PD in 3D space. KS tests were used to assess signs of clustering relative to the uniform distributions. We favored the broadly applicable KS test as its metric output is easy to interpret (higher KS test values relate to stronger differences between the real and simulated distributions), and this facilitates quantitative comparisons. Two clustering algorithms were employed to detect the number of clusters present at the PPP-EN interface in root cells. The first is a k-means method with a silhouette approach for estimating optimal cluster number (termed silhouette in the figures), which was implemented using the factoextra package (version 1.0.5; Kassambara and Mundt, 2017). The second one is a Bayesian Information Criterion for expectation-maximization, initialized by hierarchical clustering for parameterized Gaussian mixture models (termed mclust), which was implemented using the mclust package (version 5.4.2; Scrucca et al., 2016). In both cases, we arbitrarily defined the maximum numbers of clusters to 20, believing the 1:20 range to be biologically meaningful. In the case of the k-means algorithm, we additionally repeated the initial seed placing 100 times, in order to reduce the possibility of inaccurate clustering due to biases in initial seed placement. To determine the surface areas occupied by the identified clusters, we projected the 3D coordinates of the PD onto a 2D space using PCA. No significant loss of information in the distribution of the PD occurred (likely relating to the fact that the cell walls were mostly flat planes). The two first principal components of PCS captured more than 90% of the variance in the xyz coordinates of the original data in all cases reported here. The areas of the convex hulls delimited by the outer points of each cluster (or the outer points in general, in the case of the total surface) were extracted using the splancs package (version 2.01-40; Rowlingson and Diggle, 2017). Lastly, we estimated the cell wall thickness around PD. A guided tutorial with all the necessary code for this analysis is available at https://andreapaterlini. github.io/Plasmodesmata_dist_wall/ (last accessed March 2020). The Col-0 data sets used in this article, with corresponding models and annotations, are available on Figshare (https://doi.org/10.6084/m9.figshare.12488702.v1). They can be used as example data sets to test our pipeline. In addition to the specific packages listed above, we also employed the broader tidyverse environment (version 1.3.0; Wickham, 2017) and the data.table (version 1.12.0; Dowle and Srinivasan, 2019) and ggbeeswarm (version 0.6.0; Clarke and Sherrill-Mix, 2017) packages.

3D Visualizations
For 3D visualization, we employed both Imaris (version 8.4.2; Oxford Instruments) and Amira (version 2019.1; Thermo Scientific) imaging software. Export of features from the MIB environment is compatible with both visualization packages. For segmentations involving cellular organelles (ER, Golgi, and mitochondria), morphological features across the 3D stacks were used for organelle classification.

Accession Numbers
SB-EM image data sets are available at EMPIAR, the Electron Microscopy Public Image Archive, with the accession code EMPIAR-10442.

Supplemental Data
The following supplemental materials are available.
Supplemental Figure S1. Analysis of the SE-PPP interface in terms of spatial clustering using SB-EM data sets.
Supplemental Figure S2. Overview of the plugins developed for this article.
Supplemental Video S1. SB-EM data set with annotated PD on segmented walls.
Supplemental Video S2. Segmented cellular features in proximity of PD.