Abstract

   Primary and secondary growth of the tree stem are responsible for
   corresponding increases in trunk height and diameter. However, our
   molecular understanding of the biological processes that underlie these
   two types of growth is incomplete. In this study, we used single-cell
   RNA sequencing and spatial transcriptome sequencing to reveal the
   transcriptional landscapes of primary and secondary growth tissues in
   the Populus stem. Comparison between the cell atlas and differentiation
   trajectory of primary and secondary growth revealed different
   regulatory networks involved in cell differentiation from cambium to
   xylem precursors and phloem precursors. These regulatory networks may
   be controlled by auxin accumulation and distribution. Analysis of cell
   differentiation trajectories suggested that vessel and fiber
   development followed a sequential pattern of progressive
   transcriptional regulation. This research provides new insights into
   the processes of cell identity and differentiation that occur
   throughout primary and secondary growth of tree stems, increasing our
   understanding of the cellular differentiation dynamics that occur
   during stem growth in trees.

   Key words: scRNA-seq, spatial transcriptome, primary growth, secondary
   growth, secondary cell wall, SCW, Populus
     __________________________________________________________________

   This study reveals cell atlas and dynamic molecular maps involved in
   cambium differentiation during the primary and secondary stem growth of
   Populus by using single-cell RNA sequencing and spatial transcriptome
   sequencing. These findings shed light on the dynamics of cell identity
   and differentiation during stem growth in woody plants.

Introduction

   Primary growth occurs at the shoot apical meristem through divisions
   that begin in the procambium and give rise to primary vascular bundles.
   Secondary growth depends on periclinal divisions in the vascular
   cambium that produce secondary xylem and phloem ([29]Larson, 2012;
   [30]Zhu et al., 2018; [31]Luo and Li, 2022). Although both cambiums
   share some biological properties, cell division and differentiation in
   the vascular cambium during secondary growth differ from those in the
   apical meristem ([32]Baucher et al., 2007; [33]Carlquist, 2013). The
   secondary xylem, often known as wood, is a key pool for the
   assimilation of the majority of photosynthetic carbon products, which
   are renewable and required by humans ([34]Pan et al., 2011).
   Understanding wood formation has thus attracted significant attention,
   as wood is critical for carbon storage and production of sustainable
   green materials.

   The cellular structure and organization of wood derived from
   proliferation of the vascular cambium have been described in the
   literature ([35]Larson, 2012). The vascular cambium is thought to be a
   zone of dividing cells that contains cambial initials (stem cells) and
   phloem or xylem precursors ([36]Esau, 1977; [37]Bossinger and
   Spokevicius, 2018). Cambium zone cells are typically distinguished by
   their undifferentiated state and smaller size compared with neighboring
   cells. Cambium cell division and subsequent cell differentiation have
   been investigated in the model plant Arabidopsis ([38]Nieminen et al.,
   2015). Studies have revealed that a local maximum of the phytohormone
   auxin in the Arabidopsis root contributes to specification of a stem
   cell organizer for cambium meristem activity ([39]Zhang et al., 2019).
   By increasing cytokinin biosynthesis, overexpression of AHL15 changes
   herbaceous Arabidopsis into a plant with secondary growth. This finding
   suggests that constant cytokinin biosynthesis may be needed for
   secondary growth based on cambium division ([40]Rahimi et al., 2022).
   However, the identity and differentiation pathways of cambial cells in
   woody plants have not been fully established owing to a lack of
   sufficient molecular markers. Wood formation is linked to complex
   transcriptional programs that regulate cambium cell activity, xylem
   differentiation, and secondary cell wall (SCW) biosynthesis
   ([41]Baucher et al., 2007; [42]Luo and Li, 2022). However, owing to
   technical difficulties in monitoring molecular activities at the
   cellular level during tree stem growth, studies of wood formation have
   primarily concentrated on tissue-level transcriptional activities and
   reverse genetics. It remains to be seen how transcriptional activities
   in tree secondary growth can be resolved dynamically and sequentially
   at the cellular level. Single-cell RNA sequencing (scRNA-seq) has
   recently revolutionized our understanding of gene expression in
   individual cells. The use of scRNA-seq methods to analyze wood-forming
   tissue at a single-cell resolution led to identification of
   heterogeneous cell transcriptomes that could be classified into
   different cell clusters ([43]Chen et al., 2021; [44]Li et al., 2021;
   [45]Tung et al., 2023). scRNA-seq of protoplasts prepared from
   differentiating xylem of Populus revealed that the xylem tissue was
   divided into 12 cell clusters ([46]Li et al., 2021). Another study
   identified 20 cell clusters using RNA-seq of protoplasts from highly
   lignified poplar stems ([47]Chen et al., 2021). Because of the various
   parameters used to identify cell clusters, the same tissue can be
   divided into different numbers of clusters. Cell-cluster identification
   is based primarily on the reduction of linear dimensionality
   via principal-component analysis (PCA) of cell transcriptomes obtained
   from single-cell samples. Tissues are organized from various cell types
   with varying developmental statuses, and the dimensional complexity of
   cell transcriptomes is influenced by both cell type and developmental
   state. It is therefore challenging to precisely map cell types and
   dynamic cell status through scRNA-seq data processing. In this study,
   we constructed a cell atlas of primary growth tissue (PGT) and
   secondary growth tissue (SGT) from the Populus stem. Because these
   tissues show both similarities and differences in many aspects of cell
   identification and differentiation, comparison of their cell atlases
   and differentiation trajectories provided cross-verification of
   cell-cluster identification and cell development. The combination
   of scRNA-seq analysis and spatial transcriptome information aided in
   validation of the cell atlas and spatial differentiation. The results
   of this combined approach provide new insights into the cellular atlas
   and differentiation of tree stem growth.

Results

Single-cell transcriptome analysis of Populus tissues undergoing primary and
secondary growth

   Primary growth of the tree stem depends on the shoot apical meristem,
   whereas secondary growth results from proliferation of the vascular
   cambium ([48]Prassinos et al., 2005; [49]Zhu et al., 2018). We
   harvested the first to third internodes of the Populus stem tip, which
   contain PGT ([50]supplemental Figure 1A–1C), and the fifth to tenth
   internodes of the stem, which contain SGT ([51]supplemental Figure 1A
   and 1E–1G), in order to examine transcriptomic activities of the two
   growth processes at single-cell resolution. Protoplasts were prepared
   according to a an optimized procedure ([52]Lin et al., 2014), and
   high-purity, viable protoplasts were separated from the stem tissues.
   The 10x Genomics Next GEM Single Cell Kit was used to create scRNA-seq
   libraries, which were then sequenced on the Illumina NovaSeq 6000
   platform. After cell- and gene-level filtration, we recovered 7416 PGT
   cells and 11 769 SGT cells, which showed expression of about 27 000 and
   29 317 genes, respectively ([53]supplemental Table 1). After linear
   dimensional reduction, we used the Clustree package to examine dynamic
   changes in clusters under various resolution parameters. We looked at
   the total number of clusters, the number of clusters that merged into
   one, and the number of clusters that split into different clusters to
   determine the best resolution parameter for cluster identification.
   Optimized resolution parameters were used to ensure that the cell types
   were correctly clustered together and visualized by Uniform Manifold
   Approximation and Projection (UMAP) ([54]Becht et al., 2018). The
   analysis revealed 12 cell clusters in each of the two tissues
   ([55]Figure 1A and 1B). Genes that were specifically expressed in one
   or two clusters were identified as cluster-enriched genes
   ([56]supplemental Tables 2 and [57]3). We examined the expression of
   marker genes in distinct clusters to identify the detected clusters
   with the correct cell type ([58]Figure 1C and 1D). The marker genes
   included the cambium cell markers WUSCHEL-RELATED HOMEOBOX 4 (WOX4) and
   PHLOEM INTERCALATED WITH XYLEM (PXY) ([59]Suer et al., 2011); the xylem
   precursor markers ACAULIS 5 (ACL5), LONELY GUY 3 (LOG3), AUXIN RESPONSE
   FACTOR 5 (MP), and PHABULOSA (PHB) ([60]Clay and Nelson, 2005;
   [61]Prigge et al., 2005; [62]Vidaurre et al., 2007; [63]Kuroha et al.,
   2009); the phloem precursor markers AINTEGUMENTA (ANT) and BETA-AMYLASE
   3 (BAM3) ([64]Rodriguez-Villalon et al., 2014); the xylem fiber markers
   SECONDARY WALL-ASSOCIATED NAC DOMAIN 1 (SND1) and MYB DOMAIN PROTEIN
   103 (MYB103) ([65]Zhong et al., 2006, [66]2008); the vessel markers
   XYLEM CYSTEINE PEPTIDASE 1 (XCP1), XCP2, ENDO-BETA-MANNASE 6 (MAN6),
   and VASCULAR-RELATED NAC-DOMAIN PROTEIN 1 (VND1) ([67]Funk et al.,
   2002; [68]Avci et al., 2008; [69]Zhao et al., 2013; [70]Zhou et al.,
   2014); the sieve markers CALLOSE SYNTHASE 7 (Cals7), ALTERED PHLOEM
   DEVELOPMENT (APL), and SIEVE-ELEMENT-OCCLUSION-RELATED 1 (SEOR1)
   ([71]Bonke et al., 2003; [72]Pelissier et al., 2008; [73]Barratt
   et al., 2011; [74]Xie et al., 2011); and the phloem parenchyma markers
   USUALLY MULTIPLE ACIDS MOVE IN AND OUT TRANSPORTERS 12 (UMAMIT12),
   UMAMIT19, and UMAMIT22 ([75]Kim et al., 2021). In addition, spatial
   transcriptomics analysis revealed that expression of 3-KETOACYL-COA
   SYNTHASE 10 (KCS10) was mapped to the epidermis and that of CCR-LIKE
   (CCL) and PHOTOSYSTEM II SUBUNIT X (PSBX) was mapped to the cortex and
   pith ([76]supplemental Table 4). Cell-type assignment for cluster
   identification was studied further using GO pathway analysis.
   Expression of genes in the secondary cell wall biogenesis pathway was
   significantly concentrated in fiber and vessel cell clusters, whereas
   phloem development genes were expressed in sieve companion cells
   ([77]supplemental Figure 1H). The analysis led to identification of 12
   cell clusters in PGT. Nine clusters were identified as cambium, xylem
   precursor, phloem precursor, sieve, vessel, epidermis, ray parenchyma,
   cortex, and pith cells. Cluster 5 contained two cell types: cortex and
   pith. Both clusters 4 and 8 contained vessel cells, but cluster 4
   represented early-stage cells and cluster 8 represented late-stage
   cells ([78]Figure 1A and 1B). However, clusters 0, 6, and 10 in PGT did
   not have a particular cell type assigned to them ([79]Figure 1A).
   Twelve cell clusters were also identified in SGT. Eight were assigned
   to specific cell types: cambium, xylem precursor, phloem precursor,
   sieve, phloem parenchyma, vessel, fiber, and ray parenchyma
   ([80]Figure 1B). Genes related to mitosis and the cell cycle were
   highly expressed in clusters 6 and 9, indicating that these clusters
   contained a high density of dividing cells. Clusters 1, 6, 9, and 11
   were not assigned to a specific cell type ([81]Figure 1B).

Figure 1.

   [82]Figure 1
   [83]Open in a new tab

   Identification of cell clusters from the primary growth stem and
   secondary growth stem.

   (A) UMAP representation of 12 PGT clusters. Each dot corresponds to a
   single cell. Cell clusters are labeled with numbers and different
   colors. PGT stands for primary growth tissue.

   (B) UMAP representation of 12 SGT clusters. Each dot corresponds to a
   single cell. Cell clusters are labeled with numbers and different
   colors. SGT stands for secondary growth tissue.

   (C) Expression profiles of representative marker genes in each PGT
   cluster. The size of the dot indicates the percentage of cluster cells
   that express a given gene. Color indicates the average expression
   level.

   (D) Expression profiles of representative marker genes in each SGT
   cluster. The size of the dot indicates the percentage of cluster cells
   that express a given gene. Color indicates the average expression
   level.

Spatial transcriptome sequencing of the Populus stem

   We performed spatial transcriptome sequencing of the Populus stem to
   characterize cell transcriptomes in cross-sections of growing stem
   tissue. The stem was frozen and sectioned for spatial sequencing on the
   10x Genomics Visium platform ([84]Figure 2A). After filtration, we
   obtained 2406 spots, with a median of 1962 genes sequenced for each
   spot, covering transcription of a total of 23 886 genes
   ([85]supplemental Table 1). The spots were divided into six clusters
   ([86]Figure 2B) after linear regression and UMAP visualization of the
   spatial transcriptomes. The sequenced spots were mapped to the
   corresponding regions of the stem cross-section ([87]Figure 2C) and
   were assigned to the pith, parenchyma, xylem and cambium, phloem,
   cortex, and epidermis ([88]Figure 2B–2D). Expression of cell-type
   marker genes in the six clusters corresponded to their locations in the
   cross-section ([89]Figure 2E). Specific cell types were considerably
   enriched in their corresponding marker genes: the vessel markers XCP1
   and XCP2, the xylem precursor marker ACL5, the cambium marker WOX4, and
   the phloem markers SEOR1 and CLAVATA3/ESR-RELATED 41 (CLE41)
   ([90]Figure 2F). GO pathway analysis revealed that the majority of
   lipid transport genes, lignin biosynthesis genes, and phloem
   development genes were expressed in the epidermis, xylem, and phloem,
   respectively ([91]Figure 2G). The results showed a high degree of
   consistency among the anatomical regions of gene expression, the cell
   type identifications, and the single-cell transcriptomes.

Figure 2.

   [92]Figure 2
   [93]Open in a new tab

   Reconstruction of the cellular organization of the Populus stem based
   on spatial transcriptomics.

   (A) Frozen sections of Populus stem were stained with hematoxylin and
   eosin.

   (B) UMAP display of cell clusters identified by spatial
   transcriptomics. Different colors represent cell clusters.

   (C) Cell clusters mapped over images of sections stained with
   hematoxylin and eosin. Cell clusters are colored as in (B).

   (D) Paraffin-embedded sections of Populus stem stained with toluidine
   blue.

   (E) Spatial expression pattern of marker genes. The color bar indicates
   the expression level.

   (F) Expression heatmap of significantly enriched genes in each cluster.
   Representative marker genes are listed on the right. Colors reflect the
   expression level.

   (G) GO enrichment of highly enriched genes in each cluster. The size of
   the dot represents the number of genes, and the color denotes the
   adjusted P value in the related pathway.

Cambium identity

   Cambium is a type of cell that is essential for the growth and
   development of tree stems, and WOX4 and PXY are considered to be
   cambium cell markers. PXY expression was found in four PGT clusters
   (xylem precursors, early vessel cells, cluster 3, and cluster 7)
   ([94]Figure 1C), and WOX4 expression was found in clusters 3 and 7. PXY
   and WOX4 expression was also found in five SGT clusters (cell cycle
   clusters 6 and 9, xylem precursor cells, and clusters 0 and 3)
   ([95]Figure 1D). The cambium markers WOX4 and PXY were thus expressed
   in more than one cluster. Although the resolution parameter has an
   effect on the number of clusters created, clustering with different
   resolution settings had no effect on the expression of PXY and WOX4 in
   multiple cell clusters ([96]Figure 3A). To validate this result, we
   examined PXY and WOX4 expression in other studies and found that they
   were likewise expressed in various clusters ([97]supplemental
   Figure 2A). The highly enriched genes in the four clusters expressing
   PXY and WOX4 were then used for GO pathway analysis. The clusters PC3
   and SC0 were enriched in the phloem and xylem histogenesis pathway,
   whereas PC7 and SC3 were enriched in three stress-response pathways:
   response to heat, response to hydrogen peroxide, and response to
   reactive oxygen species ([98]Figure 3B). In situ hybridization was
   performed to identify the cell types that expressed the stress-response
   pathway genes ([99]Figure 3C–3E). Specific genes from the PC7 and SC3
   clusters were expressed in early ray parenchyma cells, whereas specific
   genes from PC3 and SC0 were expressed in cambium ([100]Figure 3E). We
   also examined how the cluster genes were expressed by comparing them
   with microdissection gene expression data from xylem, phloem, cambium,
   and other cell types ([101]Dai et al., 2023; [102]Tung et al., 2023).
   The PC3 and SC0 cluster genes were matched to cambium
   ([103]supplemental Figure 2B and 2C), and the genes in the PC7 and SC3
   clusters were found to be associated with fiber, vessel, and ray cells
   ([104]supplemental Figure 2D and 2E). WOX4 and PXY were also matched to
   ray cells ([105]supplemental Figure 2F). Thus, PC3 and SC0 were assumed
   to be cambium initials for xylem and phloem cells, but PC7 and SC3 were
   assumed to be ray initial cells. Because ray cells express WOX4 and
   PXY, as well as stress-response pathway genes, it would be worthwhile
   to determine the role played by these cells during stem secondary
   growth.

Figure 3.

   [106]Figure 3
   [107]Open in a new tab

   Identification of cambium.

   (A) Clustree of clustering in Seurat at various resolution settings.
   The size of the circle represents the number of cells in each cluster.
   The color of the circle shows the resolution setting. The transparency
   of the arrows represents the change in proportion between clusters with
   different resolution settings. In this study, a resolution setting of
   0.3 was selected. PGT is on the left and SGT on the right.

   (B) GO enrichment of the highly enriched genes in four clusters (PC3,
   cluster 3 in PGT; PC7, cluster 7 in PGT; SC0, cluster 0 in SGT; SC3,
   cluster 3 in SGT). The size of the dot represents the number of genes,
   and the color denotes the adjusted P value in the related pathway.

   (C) Violin plots showing expression levels of PC7 cluster genes
   (Potri.008G023866, Potri.008G181200, ATAF1a, and MEF7c) and PC3 cluster
   genes (CBSX5c, PLAT1a, CEL3, and Potri.018G140500).

   (D) Violin plots showing the expression levels of SC0 cluster genes
   (Potri.017G128800, Potri.002G113000, ERF16b, and ERF109b) and SC3
   cluster genes (NPY2a, PIN6b, Potri.006G215000, and Potri.016G051600).

   (E) RNA in situ hybridization. Scale bars correspond to 50 μm. ATAF1a
   expression in PC7, ERF16b expression in SC3, CEL3 expression in PC3,
   and NPY2a expression in SC0 are shown.

Trajectory of xylem and phloem differentiation from cambium

   Differentiation of the procambium during primary growth and the
   vascular cambium during secondary growth into xylem and phloem results
   in development of the vascular structure. Pseudo-time analysis was
   performed to arrange cells with comparable developmental transcriptomes
   in order and thus examine the trajectory of cambium differentiation in
   PGT and SGT. To construct pseudo-time trajectories, we first collected
   cambium, xylem precursor, and phloem precursor cell clusters from PGT
   and SGT. As shown in [108]Figure 4A, cambium cells, xylem precursors,
   and phloem precursors were built sequentially, with cambium serving as
   the starting point and xylem and phloem precursors differentiating into
   two branches in PGT. The trajectories were mapped by computing gene
   changes to arrange the cells in an ordered sequence. Branch-dependent
   genes whose expression was related to cell differentiation into xylem
   or phloem precursors in PGT were mapped ([109]Figure 4B). The
   expression of branch-dependent genes showed three patterns: increased
   expression in phloem precursors, increased expression in xylem
   precursors, or high expression in cambium followed by decreasing
   expression during differentiation. Enriched GO pathways in these genes
   are listed in [110]Figure 4C. The xylem branch-dependent genes included
   xylem precursor marker genes. For example, expression of ACL5, which
   controls Arabidopsis xylem specification, and MP, which encodes a major
   auxin response transcription factor in xylem ([111]Muñiz et al., 2008),
   increased steadily toward the xylem precursor branch ([112]Figure 4D).
   Homologs of LOB DOMAIN-CONTAINING PROTEIN 4 (LBD4), NAKED PINS IN YUC
   MUTANTS 2 (NPY2), and BRASSINOSTEROID-SIGNALING KINASE 6 (BSK6) were
   also strongly expressed in the xylem precursor branch ([113]Figure 4D).
   Another study ([114]supplemental Figure 3A) also showed that these
   genes were involved in xylem cell development ([115]Chen et al., 2021).
   Homologs of ARABIDOPSIS RESPONSE REGULATOR 4 (ARR4), a negative
   regulator of cytokinin signaling ([116]To et al., 2007), and LACCASE 6
   (LAC6) were discovered in the xylem precursor branch of PGT
   ([117]Figure 4D). By contrast, these genes were rarely expressed in the
   SGT trajectory ([118]supplemental Figure 3B). Tissue transcriptome
   sequencing revealed that these genes were considerably more abundant in
   PGT than in SGT ([119]supplemental Figure 3C and [120]supplemental
   Table 5). This suggested that these genes have a role in PGT xylem
   differentiation. As cells moved toward the phloem precursor branch, the
   phloem precursor markers CLE41, which encodes the TDIF peptide for
   control of cambium proliferation ([121]To et al., 2007), and BAM3, a
   receptor kinase ([122]Depuydt et al., 2013), gradually increased in
   expression ([123]Figure 4E). KNOTTED-LIKE HOMEOBOX OF ARABIDOPSIS
   THALIANA 1 (KNAT1), PA2, MEE32, GH9B13, and PER73 were also phloem
   precursor branch-dependent genes ([124]Figure 4E). Surprisingly, these
   genes showed almost no expression in SGT ([125]supplemental Figure 3D),
   but tissue RNA-seq showed that they were highly expressed in PGT
   ([126]supplemental Figure 3E). The trajectory of cambium
   differentiation in SGT revealed the differentiation path of xylem and
   phloem precursor cells from cambium ([127]Figure 4F). [128]Figure 4G
   shows the mapping of the branch-dependent genes whose expression was
   linked to cell differentiation into xylem or phloem precursors in SGT.
   As in PGT, the expression of branch-dependent genes in SGT showed three
   patterns: increased expression in phloem precursors, increased
   expression in xylem precursors, and high expression in the cambium
   followed by decreasing expression during differentiation. GO analysis
   suggested the likely biological processes involved ([129]Figure 4H). We
   next compared xylem precursor branch genes in PGT with those in SGT and
   discovered that more than half of them overlapped, including ACL5, MP,
   LBD4, NPY2, and PSBK6 ([130]Figure 4I and 4J). However, homologs of
   PLASMA MEMBRANE INTRINSIC PROTEIN 2;3 (PIP2;3) and several genes of
   unknown function were significantly expressed in SGT but not in PGT
   ([131]Figure 4J; [132]supplemental Figure 3F and 3G). CLE41, BAM3, and
   KNAT1 were strongly expressed in phloem precursor branches of both PGT
   and SGT ([133]Figure 4K and 4L). Expression of REMORIN GROUP 4.2
   (REM4.2), SENESCENCE-ASSOCIATED GENE 20 (SAG20), DROUGHT-INDUCED 21
   (DI21), and MORE LATERAL GROWTH1 (MOL1) was high in SGT but low in PGT
   ([134]Figure 4L; [135]supplemental Figure 3H). Tissue transcriptome
   sequencing confirmed this disparity ([136]supplemental Figure 3I).
   These findings imply that these genes have a role in the
   differentiation of phloem precursor cells in SGT. We also compared the
   expression of branch-dependent genes with previously published data
   ([137]Chen et al., 2021). The appropriate cell type was assigned to 153
   xylem precursor branch genes and 125 phloem precursor branch genes
   ([138]supplemental Tables 6 and [139]7). For example, HB8, PHB, and
   LBD4 and the auxin-related genes ACL5, MP, LIKE AUXIN-RESISTANT 2
   (LAX2), INDOLE-3-ACETIC ACID INDUCIBLE 9 (IAA9), IAA11, IAA13, and
   IAA19, which reflect the auxin activation signaling pathway, were
   enriched in xylem precursors ([140]Figure 4C and 4H; [141]supplemental
   Table 6), whereas ANT, SMAX1-LIKE 5 (SMXL5), and AINTEGUMENTA-LIKE 5
   (AIL5) were expressed in phloem precursors in previously published data
   ([142]supplemental Table 7) ([143]Chen et al., 2021). In summary, xylem
   and phloem differentiation in PGT and SGT involves both shared
   molecular networks and PGT- and SGT-specific molecular modules. Further
   studies of their functions could increase our understanding of cambial
   differentiation during tree growth.

Figure 4.

   [144]Figure 4
   [145]Open in a new tab

   Trajectories of cambium differentiation toward xylem and phloem
   precursors in PGT and SGT.

   (A) Pseudo-time trajectory of cambium differentiation in PGT. The
   pseudo-time of the trajectory (above) and cell-type differentiation
   (Ca, cambium; XP, xylem precursor; PP, phloem precursor) in the
   trajectory (below) are shown. Each dot represents a cell.

   (B) Heatmap depicting the expression dynamics of branch-dependent genes
   in PGT along the differentiation trajectory toward xylem and phloem
   precursors on two sides. The branch-dependent genes show three
   differentiation patterns.

   (C) Bar chart of GO enrichment of branch-dependent genes. The top 5
   pathways were identified for the three differentiation patterns in (B).

   (D) Expression of xylem precursor branch genes along the
   differentiation trajectory in PGT. Each dot represents a cell. The
   color bar indicates the expression level.

   (E) Expression of phloem precursor branch genes along the
   differentiation trajectory in PGT. Each dot represents a cell. The
   color bar indicates the expression level.

   (F) Pseudo-time trajectory of cambium differentiation in SGT. The
   pseudo-time of the trajectory (above) and cell-type differentiation
   (Ca, cambium; XP, xylem precursor; PP, phloem precursor) in the
   trajectory (below) are shown. Each dot represents a cell.

   (G) Heatmap depicting the expression dynamics of branch-dependent genes
   in SGT along the differentiation trajectory toward xylem and phloem
   precursors on two sides. The branch-dependent genes show three
   differentiation patterns.

   (H) Bar chart of GO enrichment of branch-dependent genes. The top 5
   pathways were identified for the three differentiation patterns in (G).

   (I) Venn diagram showing numbers of shared and phase-specific expressed
   genes for xylem precursor branch genes. Green indicates genes in PGT,
   and purple indicates genes in SGT.

   (J) Expression of xylem precursor branch genes along the
   differentiation trajectory in SGT. Each dot represents a cell. The
   color bar indicates the expression level.

   (K) Venn diagram showing numbers of shared and phase-specific expressed
   genes for phloem precursor branch genes. Green indicates genes in PGT,
   and purple indicates genes in SGT.

   (L) Expression of phloem precursor branch genes along the
   differentiation trajectory in SGT. Each dot represents a cell. The
   color bar indicates the expression level.

Differential expression of genes related to auxin distribution during
vascular development

   Auxin plays a critical role in regulating vascular development. Auxin
   accumulation and distribution are affected by genes encoding auxin
   transporters from the AUXIN1/LIKE-AUX1 (AUX/LAX), PIN-FORMED (PIN), and
   P-GLYCOPROTEIN/ATP-BINDING CASSETTE B4 (PGP/ABCB) families
   ([146]Geisler et al., 2005; [147]Petrásek et al., 2006; [148]Cho
   et al., 2007). We examined the expression of various genes from these
   families in the scRNA-seq data to investigate the connection between
   auxin transporters and cambium differentiation and division in the
   Populus stem. Both the cell differentiation trajectories in PGT and SGT
   and the expression of auxin transporter genes, including AUX1a/b/c,
   LAX2a, LAX3a/b, PIN6b, and PIN3a/b, were examined ([149]Figure 5;
   [150]supplemental Figure 2). LAX3a was significantly expressed in the
   cambium throughout the pseudo-time differentiation trajectory, and its
   expression decreased as the cambium differentiated into xylem or phloem
   precursors in PGT and SGT ([151]Figure 5). PGT and SGT also showed
   enhanced LAX2a and AUX1b expression as xylem identity was established
   ([152]Figure 5). As xylem precursors differentiated, AUX1c expression
   increased in SGT but remained unchanged in PGT ([153]Figure 5). When
   phloem identity was established, expression of PIN3a/b and LAX3b
   increased in PGT ([154]Figure 5A) but not in SGT ([155]Figure 5B).
   Differentiation did not cause dynamic changes in PIN6 or AUX1a in PGT
   or SGT ([156]supplemental Figure 2A). Spatial transcriptome
   localization revealed that auxin transporters were expressed in various
   regions of the stem cross-section ([157]supplemental Figure 2A). PIN6,
   PIN3, AUX1, LAX2, and LAX3 were all expressed in the cambium region
   (cambium, phloem, and xylem precursor cells), whereas LAX2a, AUX1b, and
   AUX1c were expressed mainly at xylem cell sites ([158]supplemental
   Figure 2B). These findings suggest that auxin accumulation and
   distribution during cambium differentiation may be controlled by
   several auxin transporters, some of which may differ between PGT and
   SGT.

Figure 5.

   [159]Figure 5
   [160]Open in a new tab

   Dynamics of auxin transporter expression in PGT and SGT through the
   cambium differentiation trajectory.

   (A) Auxin transporter expression dynamics over the differentiation
   trajectory in PGT. A dot represents a cell, and the colors denote the
   cell-type cluster.

   (B) Auxin transporter expression dynamics over the differentiation
   trajectory in SGT. A dot represents a cell, and the colors denote the
   cell-type cluster.

   The solid line indicates the phloem precursor branch, and the dotted
   line indicates the xylem precursor branch.

Molecular components of SCW thickening in fibers and vessels

   SCW thickening is the basic process of wood formation. Different
   transcriptional networks are thought to control SCW biosynthesis in
   various cell types of wood tissue ([161]Zhu and Li, 2021; [162]Luo and
   Li, 2022). However, the specific transcriptional activity that controls
   SCW thickening in different cell types has yet to be established.
   Analysis of scRNA-seq data from developing wood revealed a number of
   potential markers for SCW formation in fibers and vessels
   ([163]supplemental Figure 5A and 5B). Several genes encoding SCW
   biosynthesis enzymes exhibited a cell-type preference or specificity.
   For example, IRREGULAR XYLEM 14 (IRX14a) xylan biosynthesis genes were
   found only in fibers ([164]supplemental Figure 5A), whereas
   PHENYLALANINE AMMONIA LYASE 1 (PAL1a) lignin biosynthesis-related genes
   were found only in vessels ([165]supplemental Figure 5A). LAC17 has
   four homologs in Populus; LAC17a and LAC17d were expressed primarily in
   vessel cells, and LAC17c was expressed only in fiber cells
   ([166]supplemental Figure 5A and 5B). Numerous SCW synthesis enzymes
   were expressed in both fibers and vessels without regard for cell type
   ([167]supplemental Figure 5C). These data imply that SCW biosynthesis
   may be regulated differently in different cell types. ROPGAP3, PLANT
   GLYCOGENIN-LIKE STARCH INITIATION PROTEIN 5 (PGSIP1), and
   PLASMODESMATA-LOCATED PROTEIN 6 (PDLP6) were specifically expressed in
   fibers ([168]supplemental Figure 5A), whereas GERMIN-LIKE PROTEIN 7
   (GLP7) and PIP2 were specifically expressed in vessels
   ([169]supplemental Figure 5B). The particular cell-type preference of
   these genes has also been observed in other Populus scRNA-seq studies
   ([170]supplemental Figure 6A–6F) ([171]Li et al., 2021; [172]Tung
   et al., 2023). Although the cell-type specificity of genes must be
   confirmed through further cellular studies, identification of new
   cell-type markers would be useful for obtaining novel insights into
   specific cell-wall biosynthesis in xylem and phloem.

   We next analyzed the cell subpopulations in this study. We began by
   extracting vessel and fiber clusters. A combination of mutual nearest
   neighbors ([173]Haghverdi et al., 2018) and canonical correlation
   analysis ([174]Butler et al., 2018) methods in the Seurat package
   (v4.1.3) ([175]Satija et al., 2015) was used for cell integration. The
   Clustree package was used to examine cluster changes under various
   resolution options. Highly enriched genes in subgroups were examined at
   several resolutions, and the resolution with the clearest association
   between gene expression and recognized developmental processes was
   chosen for subgroup analysis. There were four subgroups in the cell
   cluster designated as vessels (V1, V2, V3, and V4) ([176]Figure 6A),
   and substantially enriched genes in each subgroup corresponded to the
   four successive stages of vessel development ([177]Figure 6B;
   [178]supplemental Table 8). Expression of WOX4, ACL5, TRANSLOCASE OF
   THE OUTER MITOCHONDRIAL MEMBRANE 6 (TOM6), and other markers for
   cambium and xylem precursor cells revealed that the V1 subgroup
   retained these traits ([179]Figure 6B). In V2, the xylem precursor
   genes ACL5 and AUX1 were expressed, but the cambium gene WOX4 was not
   ([180]Figure 6B). At the same time, expression of the transcription
   factor MYB46, a crucial regulator for SCW formation, started the
   thickening of SCWs in the vessels ([181]Figure 6B). Prior to SCW
   thickening, the EXPANSIN (EXPA) and XYLOGLUCAN ENDOTRANSGLUCOSYLASE
   (XTH) family of genes, which is also expressed in V2, is needed to
   relax the cell wall ([182]Figure 6B). GO analysis revealed that
   enriched genes at this stage had functional annotations for cell wall
   organization and cellular glucan metabolic pathway ([183]Figure 6B and
   6C). In V3, the SCW biosynthesis regulators MYB42c and MYB52e and genes
   encoding SCW biosynthesis enzymes, such as SCW-associated CELLULOSE
   SYNTHASE 8 (CESA8a) and the lignin biosynthesis-related genes PAL1c and
   4-COUMARATE:COA LIGASE 1 (4CL1a), were expressed ([184]Figure 6B).
   Expression of these genes was reflected in GO functional annotations
   such as plant-type SCW biogenesis and lignin catabolic process
   ([185]Figure 6C). In V4, SCW formation genes were continuously
   expressed, together with genes related to normal programmed cell death,
   such as XCP1a and XCP2a/b ([186]Figure 6B). In addition to genes with
   known functions in vessel development, genes from the BASIC
   HELIX-LOOP-HELIX (bHLH), BEL1-LIKE HOMEOBOX (BLH), ETHYLENE RESPONSE
   FACTOR (ERF), ARABIDOPSIS NAC DOMAIN-CONTAINING PROTEIN (ANAC), LBD,
   NITRATE TRANSPORTER (NPF), and short peptide RAPID ALKALINIZATION
   FACTOR (RALF) families were also discovered ([187]Figure 6B) and may
   represent new vessel regulators. BHLH106b, ERF105b, and NRT1.1 were
   highly enriched in the V1 subgroup ([188]Figure 6B); BHLH145b, ERF003b,
   and NPF4.4a in the V2 subgroup ([189]Figure 6B); BLH1b, ERF109c, and
   NPF5.2a in the V3 subgroup ([190]Figure 6B); and BHLH128b, ERF110b, and
   NPF5.2a in the V4 subgroup ([191]Figure 6B). In addition, V3 and V4
   expressed ANAC087b, which is involved in regulating the onset of cell
   death ([192]Huysmans et al., 2018). This expression was corroborated by
   the spatial transcription location ([193]Figure 6B and 6D).

Figure 6.

   [194]Figure 6
   [195]Open in a new tab

   Vessel development analysis.

   (A) UMAP analysis showing that the vessel cell cluster is separated
   into four subgroups. The subgroups are color coded.

   (B) Heatmap depicting the highly enriched genes in each vessel
   subgroup. Two columns list representative genes. Genes with known and
   unknown functions in vessel development are listed separately. The
   color of the dots indicates the expression level.

   (C) GO enrichment of the most enriched genes in each vessel subgroup.

   (D) Spatial expression pattern of vessel-enriched genes in the stem
   cross-section. The color of the dots indicates the expression level.

   Four subgroups were identified in the fiber cell cluster, and enriched
   genes in each subgroup corresponded to four successive developmental
   stages ([196]Figure 7A; [197]supplemental Table 9). Three (F1, F2, and
   F4) appeared to correspond with the progression of fiber development,
   but subgroup F3 did not clearly demonstrate a link with developmental
   phases ([198]Figure 7B). Cells in the F1 subgroup expressed the xylem
   precursor genes PHBa and KNAT7, as well as the SCW formation regulators
   MYB46b and MYB83, indicating that SCW development had begun in this
   subgroup ([199]Figure 7B). Expression of the SCW
   biosynthesis-regulating transcription factors MYB52c/e and the enzyme
   genes CESA7b, LAC17c, and 4CL1a was enriched in the F2 subgroup
   ([200]Figure 7B). More SCW synthesis enzyme genes, particularly those
   involved in lignin production such as PAL1b, 4CLa, and LAC17b, and more
   transcription factors regulating SCW synthesis, including SND2b/d,
   MYB52b, and MYB85, were significantly enriched in the F4 subgroup
   ([201]Figure 7B). GO functional annotations for plant-type SCW
   biogenesis also indicated the status of F2 and F4 ([202]Figure 7C).
   These findings suggested that fiber-related genes act sequentially in
   fiber formation. In addition, we discovered several transcription
   factors and genes involved in signal transduction that may be
   recognized as new candidate regulators in fiber formation. These
   included RALFL33a/b/c, BLH3, BHLH62d, and ROPGEF7b in F1; BLH7a,
   BHLH147a, and ROPGAP3b in F2; and BSK1a and TORTIFOLIA 1 (TOR1b) in
   F4 ([203]Figure 7C). Fibers selectively expressed TOR1b, which
   regulates microtubule-dependent cell elongation ([204]Furutani et al.,
   2000) ([205]Figure 7B). BLH3, which controls the timing of the
   vegetative-to-reproductive phase transition ([206]Zhang et al., 2016),
   was found in the F1 subgroup ([207]Figure 7B). The ROP activator gene
   ROP (RHO OF PLANTS) GUANINE NUCLEOTIDE EXCHANGE FACTOR 7 (ROPGEF7b),
   which is involved in determining the metaxylem vessel pit pattern
   ([208]Nagashima et al., 2018), was particularly upregulated in the F1
   subgroup ([209]Figure 7B). Spatial transcriptomics revealed where these
   genes were expressed on the stem cross-section ([210]Figure 7D). By
   contrast, F3 lacked expression of any known xylem marker genes, and its
   enriched GO pathways did not involve cell wall organization or SCW
   biogenesis ([211]Figure 7C). However, the F3 subgroup had a greater
   number of genes that were associated with reactions to heat, stress,
   and reactive oxygen ([212]Figure 7C), and this subgroup overlapped with
   the other three subgroups in the UMAP graphic ([213]Figure 7A). The F3
   subgroup cells were grouped together within the fiber cells, but they
   exhibited numerous fiber-cell-specific traits. To understand how the
   subgroup cells are incorporated in subsequent SCW thickening, more
   research would be highly worthwhile.

Figure 7.

   [214]Figure 7
   [215]Open in a new tab

   Fiber development analysis.

   (A) UMAP analysis showing that the fiber cell cluster is separated into
   four subgroups. The subgroups are color coded.

   (B) Heatmap depicting the highly enriched genes in each fiber subgroup.
   Two columns list representative genes. Genes with known and unknown
   functions in fiber development are listed separately. The color of the
   dots indicates the expression level.

   (C) GO enrichment of the most enriched genes in each fiber subgroup.

   (D) Spatial expression pattern of fiber-enriched genes in the stem
   cross-section. The dot color indicates the expression level.

Discussion

   The majority of the carbon fixed during photosynthesis in trees is
   converted and stored in their secondary xylem to produce wood, an
   essential and renewable biomaterial for human needs. Despite thorough
   descriptions of wood cell organization and structure in trees
   ([216]Larson, 2012), we lack a thorough understanding of how cells
   differentiate to form wood. Here, we attempted to describe cellular
   atlases and regulatory networks of wood formation by building and
   comparing cell developmental trajectories in wood-forming tissues of
   Populus. The findings provide fresh insight into the processes of
   cellular differentiation during tree stem growth.

Single-cell transcriptomics of Populus stem development

   The recent invention of scRNA-seq has revolutionized our understanding
   of transcriptional activity within tissues. However, few scRNA-seq
   experiments have been performed in plants, and the majority of these
   studies have used herbaceous Arabidopsis tissues. Spatial transcriptome
   sequencing has been used in animal studies, but it has only recently
   been applied to the study of plants ([217]Liu et al., 2022; [218]Xia
   et al., 2022). Perennial trees are distinguished by cambium
   proliferation and stem growth differentiation. Tree stem growth is
   supported by two cambium meristems, one in the shoot tip and the other
   in the vascular cambium ring. Similar and unique characteristics in
   cell identification and differentiation can be found in both meristem
   tissues ([219]Baucher et al., 2007; [220]Campbell and Turner, 2017).

   We acquired scRNA transcriptome datasets from primary and secondary
   growth tissues. Identification of cell clusters and cell development
   were cross-verified by comparing the two datasets and spatial
   transcriptome data. Twelve cell clusters were identified in both PGT
   and SGT by linear dimensional reduction and resolution parameter
   optimization. Cell-type markers for cambium, xylem precursor, phloem
   precursor, sieve, and vessel cells were found in cell clusters of both
   PGT and SGT. Markers for phloem parenchyma and fiber cells were found
   only in SGT cell clusters, and markers for epidermis, cortex, and pith
   cells were found only in PGT cell clusters. Through spatial
   transcriptome analysis, expression of marker genes could be localized
   to appropriate cells in stem cross-sections. Three cell clusters in PGT
   and four in SGT could not be classified into their relevant cell types
   owing to the absence of suitable markers. These results suggested that
   there were both similarities and differences in cell arrangement
   between primary and secondary growth. scRNA transcriptomes revealed the
   expression of approximately 27 000 genes in PGT and 29 317 genes in
   SGT; however, the expression of genes abundant in cell clusters varied
   between PGT and SGT. Comparison of scRNA transcriptome datasets from
   the two tissues could facilitate identification of essential regulators
   in the two developmental stages.

Transcriptional control of cambium differentiation

   In secondary growth, cell differentiation after cambium proliferation
   is essential for expansion of stem diameter. Studies have shown that a
   complex network of transcriptional programs and hormone signals
   regulates differentiation ([221]Ye and Zhong, 2015; [222]Luo and Li,
   2022). The majority of these studies have focused on tissue-level
   molecular dissection, and it remains unknown how differentiation is
   dynamically and sequentially regulated at the cellular level. Here,
   construction of pseudo-time trajectories revealed that
   branch-determining genes control cambium differentiation in two
   directions: toward xylem and phloem. The auxin transporters AUX1b,
   LAX2, IAA9, IAA11a, IAA13a/b, and IAA19a/b were dynamically expressed
   in the differentiation trajectory toward xylem in both PGT and SGT,
   suggesting that auxin signaling and distribution are crucial for xylem
   differentiation. Impaired xylem differentiation in aux1lax1lax2lax3
   mutants of Arabidopsis supports this conclusion ([223]Fabregas et al.,
   2015). Auxin distribution shows a peak in the vascular cambium of
   Populus ([224]Immanen et al., 2016), and the VCM1–PIN5 molecular module
   controls auxin concentration to control xylem differentiation
   ([225]Zheng et al., 2021). Dynamic regulation of auxin distribution
   during xylem development is implied by the involvement of distinct
   auxin transporters in differentiation. However, more research is needed
   to understand the genetic processes by which auxin transporters
   contribute to secondary growth in trees.

   In addition to auxin-associated genes, we discovered other
   transporters, including the SWEET sugar transporters SWEET1a, SWEET17a,
   and SWEET2a. A recent Arabidopsis study showed that the sweet16 mutant
   has fewer xylem cells and that SWEET16 is expressed at the
   procambium–xylem interface ([226]Aubry et al., 2022). Furthermore,
   overexpression of SWEET7 in Populus promotes xylem growth ([227]Zhang
   et al., 2021). All of these findings suggest that sugar transport
   influences cambium differentiation. The trajectory also revealed the
   presence of two small peptide genes, RALFL33a/b and RALFL34a/b. They
   appear to function as a signal in cambium differentiation. A large
   number of transcription factor genes also showed branch-dependent
   expression during differentiation. As potential new regulatory factors,
   these genes may participate in the regulation of cambial
   differentiation.

   The xylem is mainly composed of two cell types with thickened SCWs:
   vessel cells for mineral and water transport and fiber cells for
   mechanical support. These are also the main cell components for
   biomaterial production ([228]Ek et al., 2009; [229]Gorshkova et al.,
   2012). Although both cell types go through comparable developmental
   processes, such as cell expansion, secondary wall thickening, and
   programmed cell death, their secondary wall thickening may be
   controlled by different transcriptional regulatory networks ([230]Ye
   and Zhong, 2015; [231]Zhu and Li, 2021). In this work, we examined the
   developmental trajectories of fiber and vessel cells. Some functional
   genes were shared across fiber and vessel cells, but others were cell
   specific, such as four homologs of the laccase gene LAC17, which
   functions in lignin synthesis. Several transcription factors and
   signaling molecules were also expressed at different developmental
   stages of the two cell types, suggesting that they may have regulatory
   roles in the development of fiber and vessel cells. Overall, these
   findings contribute to a better understanding of the molecular
   regulatory networks involved in the formation of fiber and vessel cells
   and provide a basis for future wood modification via cell-specific
   mechanisms.

Understanding the growth of tree stems via cellular transcriptomes

   Tree stem growth is caused by proliferation of the vascular cambium,
   followed by cellular differentiation, and this is the primary pathway
   for storage of photosynthetically fixed carbon in lignocellulosic
   biomass in trees ([232]Pan et al., 2011). The use of scRNA-seq analysis
   allows for a cellular-level description of the molecular activity of a
   tissue, and analysis of single-cell transcriptomes in Populus stem
   tissues has recently been described ([233]Chen et al., 2021; [234]Li
   et al., 2021; [235]Tung et al., 2023). These scRNA-seq experiments
   focused on xylem tissue cells. However, because a limited number of
   marker genes are available for cell-type identification in growing stem
   tissue, cell clusters identified in scRNA-seq datasets could not be
   unambiguously assigned to specific cell types. This limitation results
   from a lack of understanding of how trees grow. To investigate cambium
   differentiation in various directions, we collected cells from both
   sides of the cambium in growing stem tissues. We were able to
   cross-check the cambium differentiation trajectories by comparing
   primary and secondary stem tissues at the same time. A large number of
   genes have been identified as participating in the process of cell
   differentiation during stem cell proliferation using scRNA-seq
   research. However, further experimental characterization is required to
   properly characterize their functional mechanisms in cambium
   differentiation and xylem development. Spatial transcriptome analysis,
   on the other hand, enables observation of cellular transcriptomes
   in situ. However, owing to its restricted resolution when applied to
   plant tissues, spatial transcriptome analysis in plants is still in its
   early stages. The combination of scRNA-seq with spatial transcriptome
   analysis therefore provides a systematic view of the molecular networks
   involved in wood formation at the cellular level.

Methods

Tree growth conditions and protoplast isolation

   A phytotron was used to cultivate Populus euramericana cv. ‘Nanlin895’
   for 2 months under a 16-h light/8-h dark photoperiod. Intact
   protoplasts were isolated according to a previously described procedure
   with some modifications ([236]Lin et al., 2014). To obtain single-cell
   transcriptomes, the first to third internodes were extracted as PGT and
   the fifth to tenth internodes as SGT. Internodes from the sampled
   plants were divided into outer bark and inner woody tissue, then
   immersed in a solution for breaking down cell walls that contained 1.5%
   (wt/vol) cellulase R-10, 0.4% (wt/vol) macerozyme R-10, 0.6 M mannitol,
   20 mM KCl, 10 mM CaCl[2], and 0.1% (wt/vol) bovine serum albumin. After
   vacuum extraction and lysis in the dark for 2 h, protoplasts were
   collected by centrifuging twice at 300 g for 3 min. This step removed
   any clumped cells. Protoplasts were resuspended in a 0.6 M mannitol
   solution. After centrifuging at 700 g for 1 min, 1 mL of the
   resuspension was added to 5 mL of a 2 M sucrose solution. Trypan blue
   staining was used to check protoplast concentration and viability.
   Protoplasts with a viability rating of >75% were used immediately for
   single-cell sequencing analysis.

scRNA-seq library construction and sequencing

   As directed by the vendor, scRNA-seq libraries were created using the
   10x Chromium Single Cell 3′ Platform (10x Genomics, Pleasanton, CA). To
   create single-cell GEMs, protoplasts were briefly injected onto a
   chromium microfluidic chip (gel beads in emulsion). The Chromium Single
   Cell 3′ library kit v3.1 was used to create scRNA-seq libraries (10x
   Genomics). An Agilent 2100 Bioanalyzer was used to analyze the
   libraries, and the Illumina NovaSeq 6000 platform was used to sequence
   the data.

Spatial transcriptome library construction and sequencing

   The fourth internode was obtained from a 2-month-old Populus tree,
   embedded in optimal cutting temperature compound, and cryosectioned.
   Hematoxylin and eosin was used to stain the sections. A cDNA library
   was constructed on the 10x Visium platform according to the
   manufacturer’s instructions, and the DNBSEQ-T7 platform was used for
   sequencing.

Cell clustering and scRNA-seq data analysis

   The Populus trichocarpa v4.1 genome was used to align the raw scRNA-seq
   data. Cell Ranger algorithms produced expression matrices for every
   gene and every cell. Additional analysis was performed using Seurat
   software (v4.1.3) ([237]Satija et al., 2015). Doublets (two or more
   cells in a single oil droplet) were removed using DoubletFinder
   (v4.1.2) ([238]McGinnis et al., 2019). A gene was counted if it was
   expressed in at least three distinct cells. Cells that expressed more
   than 10 000 genes, fewer than 500 genes, or more than 5% mitochondrial
   genes were eliminated. In the data matrix, this method produced 7416
   and 11 769 cells from PGT and SGT, respectively. We used canonical
   correlation analysis and mutual closest neighbors approaches in Seurat
   software (v4.1.3) for integration of vascular cells. The remaining
   matrix was then normalized in Seurat (v4.1.3) via LogNormalize.
   FindVariableFeatures with the selection.method = "vst" option was used
   to identify highly variable genes, and the top 2000 genes were used for
   PCA to reduce the linear dimensionality. Clustree (v0.4.4) was used to
   determine the ideal resolution parameter for cluster identification.
   Using the Seurat function FindClusters with resolution = 0.3, the first
   20 PCs were used to find cell clusters, which were then visualized and
   investigated with UMAP ([239]Becht et al., 2018). FindAllMarkers was
   implemented in Seurat (v4.1.3) with the settings “min.pct = 0.25” and
   “logfc.threshold = 0.58” to identify cluster-enriched genes
   ([240]Satija et al., 2015).

Analysis of spatial transcriptome data and cell clustering

   The Populus trichocarpa v4.1 genome was used to align the raw scRNA-seq
   data. Expression matrices for each gene and each location were created
   using Visium Ranger algorithms. Further analysis was performed with
   Seurat software (v4.1.3) ([241]Satija et al., 2015). In Seurat
   (v4.1.3), the matrix was normalized using SCTransform. PCA was used to
   reduce linear dimensionality. Using the Seurat function FindClusters
   with “resolution = 0.5”, the first 30 PCs were used to find cell
   clusters, which were then visualized and investigated with UMAP
   ([242]Becht et al., 2018). FindAllMarkers was implemented in Seurat
   (v4.1.3) with the settings “min.pct = 0.25” and “logfc.threshold =
   0.58” to identify cluster-enriched genes.

Pseudo-time analysis

   A continuous developmental trajectory was created using the R program
   Monocle2 (v4.1.1) ([243]Qiu et al., 2017). In brief, the
   cluster-enriched genes were used to select and arrange the cells of
   interest. The data were split into two components using the DDRTree
   approach, and branch-dependent genes were identified using BEAM.

Pathway analysis

   The R package clusterProfiler (v4.1.2) was used for pathway enrichment
   analysis ([244]Yu et al., 2012). The pvalueCutoff, qvalueCutoff, and
   pAdjustMethod parameters were set to 0.01, 0.05, and BH, respectively.

Gene information

   [245]Supplemental Table 10 provides JGI numbers and other information
   for all genes mentioned in this study.

In situ hybridization

   In situ hybridization was performed according to a previously published
   method ([246]Dai et al., 2023). In brief, Populus stems were fixed with
   FAA (formalin-aceto-alcohol), dehydrated in an alcohol gradient,
   infiltrated with paraffin, and finally embedded in paraffin. A
   gene-specific segment of 100–150 bp was selected and amplified by PCR
   (primers used are listed in [247]supplemental Table 11), then assembled
   into the intermediate vector, transcribed in vitro, and labeled with
   digoxin according to the manufacturer’s instructions (Roche
   10999644001; Roche 3359247910). Sections were fixed on poly-lysine
   slides. After pretreatment, probes were hybridized with sections
   overnight at 48°C, after which antibodies were incubated for the final
   color reaction.

Funding

   This research was supported by the National Natural Science Foundation
   of China (32130072), the Chinese Academy of Sciences' Strategic
   Priority Research Program (XDB27020104), and the National Key Research
   and Development Program (2021YFD2200204).

Author contributions

   L.L. conceived the project. J.W. made suggestions for the project. R.L.
   and Z.W. performed the experiments. R.L. and L.L. analyzed the data and
   wrote the manuscript. All authors discussed the results and approved
   the final manuscript.

Acknowledgments