Plant material and growth conditions
All wild type and transgenic lines used here were Arabidopsis thaliana ecotype Columbia-0 (Col-0). Seeds were surface-sterilized and sown on 0.8% agar plates containing Murashige and Skoog medium with 0.5% sucrose or liquid media as described previously14. Plants were grown under 12L/12D (12 h light and 12 h dark, 84 µmol m-2 s-1) conditions at 22°C. For induction of ectopic vascular cell differentiation, plants were entrained by 12L/12D conditions for 7 days, although the original protocols14 called for growth of plants under continuous light (LL) conditions. Bikinin, 2,4-dichlorophenoxyacetic acid (2,4-D), and kinetin were added to the 8-day-old plants at ZT0. cca1-1 lhy-11 toc1-2 and cca1-1 lhy-11 were provided by Takafumi Yamashino47 (Nagoya University), lux nox and its parental CAB2::LUC were provided by Dmitri A. Nusinow15 (Donald Danforth Plant Science Center), LUX::LUX-GFP/lux-4 was provided by Philip A. Wigge48 (University of Cambridge).
Measurement of cell differentiation, cell proliferation, and cell cycles.
For lignin staining, plants were fixed in acetic ethanol fixative (75% glacial acetic acid and 25% ethanol) for 1 day, stained with 20% phloroglucinol in 99.5% ethanol and concentrated HCl (1:19, v/v) for 1 h, cleared with chloral hydrate/glycerol/H2O mixture (8 g of chloral hydrate in 1 mL of glycerol and 2 mL of H2O) for 2 h, and observed under a light microscope. For quantification of the vascular cell induction ratio, areas for xylem cells in the cotyledons were calculated using ImageJ.
For measurement of stomatal index, 10-day-old plants grown under 12L/12D conditions were stained with 50 mg/mL of propidium iodide (PI) and observed using a confocal laser scanning microscope, FV1000 (Olympus) in three square areas of 0.48 mm2 per cotyledon from 5 cotyledons of 5 independent plants. The stomatal index was calculated as previously described49. Error bars, representing standard errors, were calculated from the results of 5 independent cotyledons.
For root elongation measurements, plants were grown vertically under 12L/12D conditions. Root length was measured at ZT0. For quantification of root meristem size, 7-day-old plants grown under 12L/12D conditions were stained with 20 mg/mL PI and observed using an FV1000 microscope as described above. Meristematic cell numbers were determined from observations of the cortical cells, using confocal microscopy images.
For GUS staining, plants were fixed in 90%(v/v) acetone for 15 min on ice, vacuum-infiltrated and incubated at 37 °C for 2 h (overnight for IRX3::GUS) in the GUS assay solution containing 100 mM sodium phosphate buffer (pH 7.2), 1 mM potassium-ferrocyanide (5 mM for TDR::GUS), 1 mM potassium-ferricyanide (5 mM for TDR::GUS), 0.1%(v/v) Triton X-100 and 0.5 mg mL−1 5-bromo-4-chloro-3-indolyl-β-D-glucuronic acid (X-Gluc). Chlorophylls in the tissue were removed by incubation in 70%(v/v) ethanol.
qRT-PCR
Total RNA was extracted using an RNeasy Plant Mini Kit (QIAGEN) and reverse-transcribed using a Transcriptor First Strand cDNA Synthesis Kit (Roche) according to the manufacturer's instructions.
Real-time gene expression was analyzed with a CFX96 Real-Time PCR Detection System (Bio-Rad). UBQ14 were used as an internal control for VISUAL experiments14, respectively. Specific sequences for each primer pair were:
CAB3-RT-F, 5'- ACCCAGAGGCTTTCGCGGAGT;
CAB3-RT-R, 5'- TGCGAAGGCCCATGCGTTGT;
LHCB2.1_Fw 5'- ACCCGGAGACATTCGCTAAGAAC;
LHCB2.1_Rv 5'- TGGATCAAGTTAGGGTTTCCGAGG;
RGF6 Fw 5'- CACAAGGTTGCTTCTTTGGACG;
RGF6 Rv 5'- CACCAACATGGTGAATGCATATGC;
HAM2 Fw 5'- CAAACGGCGGAGATAACAAT;
HAM2 Rv 5'- CTGTGGAACGGAGGTTTAGG;
TDR-RT-F, 5'- TGGTGGAAGTTACTTTGAAGGAG;
TDR-RT-R, 5'- TTCAATCTCTGTAAACCACCGTAA;
ANT_fw, 5'- TCAATACCGAGGCGTTACAAGAC;
ANT_rv, 5'- TCGAGCAGCTTTCTCCTCCATATC;
IRX3-RT-F, 5'- CCTCGGCCACAGCGGAGGAT;
IRX3-RT-R, 5'- CGCCTGCCACTCGAACCAGG;
VND7_fw, 5'- TTCGAAACGCAGTCGTATAATCC;
VND7_rv, 5'- ATTAGCTTCGACCTCATTATAGCTTTG;
CCA1-RT-F, 5'- GAGGCTTTATGGTAGAGCATGGCA;
CCA1-RT-R, 5'- TCAGCCTCTTTCTCTACCTTGGAGA;
LHY-RT-F, 5'- ACGAAACAGGTAAGTGGCGACA;
LHY-RT-R, 5'- TGCGTGGAAATGCCAAGGGT;
CYCD3;1-qPCR-Fw, 5'- CGAAGAATTCGTCAGGCTCT;
CYCD3;1-qPCR-Rv, 5'- ACTTCCACAACCGGCATATC;
E2Fc-qPCR-Fw, 5'-; GAGTCTCC-CACGGTTTCAG;
E2Fc -qPCR-Rv, 5'-; TCACCATCCGGTACTGTTGC;
UBQ14-qPCR-Fw, 5'- TCCGGATCAGGAGAGGTT; and
UBQ14-qPCR-Rv, 5'- TCTGGATGTTGTAGTCAGCAAGA.
The following thermal cycling profile was used,
CAB3, 95˚C for 10 s, ~40 cycles of 95˚C for 10 s, 62˚C for 15 s, and 72˚C for 15 s;
CCA1, 95˚C for 60 s, ~40 cycles of 95˚C for 10 s, 60˚C for 15 s, and 72˚C for 7 s;
IRX3, 95˚C for 60 s, ~40 cycles of 95˚C for 10 s, 64.5˚C for 15 s, and 72˚C for 10 s;
ANT, VND7, RGF6, and HAM2, 95˚C for 30 s, ~40 cycles of 95˚C for 5 s and 60˚C for 30 s; and
TDR, LHY, CYCD3;1, E2Fc, and UBQ14, 95˚C for 60 s, ~40 cycles of 95˚C for 10 s, 60˚C for 15 s, and 72˚C for 15 s.
Each sample was run in technical triplicate to reduce experimental errors. Error bars, representing standard errors, were calculated from the results of biological triplicates.
scRNA-seq and cpRNA-seq
For scRNA-seq, the process closely followed the method described by Kubo et al24. A cotyledon before and after treatment with VISUAL for the indicated time periods was placed adaxial side down on a glass slide, and fixed with an adhesive tape, e.g., cellophane tape. Then the center of cotyledon was cut using a razor blade, and with the aid of a microscope, the contents of a single cell were collected using a glass capillary. Samples were subjected to UMI-tagged sequencing using a NextSeq 500 system (Illumina). Detailed methods are described in Supplementary information. For cpRNA-seq, total RNA was extracted using an RNeasy Plant Mini Kit and subjected to UMI-tagged sequencing, as for scRNA-seq, except that 10 cycles of the PCR amplification step were required. Sequence reads were mapped against to the TAIR10 Arabidopsis cDNA sequence by Bowtie50 with the parameter “-a --best --strata”. Gene expression level was quantified by counting the absolute number of UMIs mapped to each gene. scRNA-seq data was normalized together with cpRNA-seq using the R/Bioconductor DESeq51 package (http://bioconductor.org/packages/release/bioc/html/DESeq.html; v.1.34.1).
Wishbone – Seurat – PeakMatch pipeline
To identify scRNA-seq data related to the xylem cell lineage, Wishbone was performed as previously described11. Samples after induction in VISUAL were subjected to Wishbone. The first and third components of principal component analysis (PCA) were used for t-SNE. The xylem lineage was selected according to the expression of xylem cell marker genes on the t-SNE plots.
Clustering of cells was performed using the Seurat R package12. In brief, digital gene expression matrices were column-normalized and log-transformed. To obtain a landmark gene set for Seurat, we divided all genes in cpRNA-seq data into 17 groups according to the peak expression time of each gene. The 17 genes showing the highest correlation coefficient with scRNA-seq data in each respective group were selected as landmark genes. In addition to the 17 genes, cell-type-specific markers (CAB3, LHCB2.1, TDR, AtHB8, IRX3, IRX8, and SEOR1) were also selected as landmark genes. In total, 24 genes were used as a landmark gene set for Seurat.
Finally, we selected the genes whose correlation coefficient between scRNA-seq and cpRNA-seq was more than 0.5. Among the selected genes, 2,217 genes whose sum total of expression levels in the scRNA-seq data were higher than the value equivalent to ten times the cell numbers were subjected to PeakMatch with the following parameters: T = 2, last = 0, intv = 1, inter = 7.
Overview of the PeakMatch algorithm
Let be the set of whole genes under consideration. Discretizing pseudo and actual times into integers for simplicity, we denote by and the sets of available pseudo and actual times, respectively. Suppose that, for each gene , we are given pseudo time-series based scRNA-seq data and actual time-series based cpRNA-seq data , where and represent the gene z's expression levels at a pseudo time p in the scRNA-seq data, and at an actual time a in the cpRNA-seq data, respectively.
To estimate the actual times of gene expressions in the scRNA-seq data, we would like to find pairs of pseudo and actual times so that the expression levels and are likely to be “comparable” for many genes . Once such pairs are found, we may estimate the actual time of by that of .
The point is that, among the observed gene expression levels, “peaks” are the most important phenomena. Then it is desired that a peak in and a peak in should be matched. It is also required that the pseudo time order should be preserved in the time pairs. To be more precise, whenever a pseudo time p is matched to an actual time a, any pseudo time p′ > p should be matched to an actual time a′ > a.
We formulated the problem of finding such time pairs as the maximum weighted non-crossing matching (MWNCM) problem for a bipartite graph. The problem is polynomially solvable52, meaning that it is efficiently solvable from the standpoint of the theory of computational complexity.
We took the bipartite graph so that one vertex subset was the pseudo time set and the other vertex subset was the actual time set . For the edge set, we considered all possible pairs , where we determined the weight of an edge heuristically by how the pseudo time p and the actual time a were comparable peaks. We determined the weight of an edge as follows. For each gene , we decided whether or not the value was within a “peak area” in . We considered that was within a peak area if it was significantly larger than a general trend of , which was estimated by an exponential moving average. We set the weight of to a larger value if both and are among peak areas for more genes.
Given the scRNA-seq and cpRNA-seq data, the algorithm PeakMatch constructed the bipartite graph, derived an MWNCM for it, and estimated the actual times of all pseudo times in based on the derived MWNCM.
For more details and python-based programs, see https://github.com/endo-lab/PeakMatch.
Plasmid construction
For IRX3::GUS, the promoter of IRX3 was amplified from Col-0 genomic DNA using the following primers:
IRX3p-Fw (XhoI), 5′-cgcggccgcactcgaTCGAGAGCCCGA; and
IRX3p-Rv (XhoI), 5′-gtctagatatctcgaAGGGACGGCCGGAGATTAGCAGCGA.
The amplified fragment was cloned into the XhoI site of pENTR1A (no ccdB). After sequence verification, this plasmid was recombined with pGWB353, and introduced into Col-0 plants.
For BES1::BES1-GFP, the genomic DNA fragment of BES1 containing the 2 kb promoter sequence and the 1 kb downstream sequence from the stop codon was amplified from Col-0 genomic DNA using the following primers:
CACC-pBES1_2k_F, 5′- CACCTCTCAACCTGCTCGGT; and
gBES1-R, 5′- CTCTGATGTGGAGTCAATG.
The amplified fragment was cloned into the pENTR/D-TOPO (Thermo Fisher Scientific). After sequence verification, this plasmid was recombined with pGWB153. The resulting pGWB1-gBES1 was linearized by PCR using the following primers:
gBES1-C_inverse-F, 5′- TGAGATGAAGTATACATGAACCTG; and
bes1_R_nonstop, 5′- ACTATGAGCTTTACCATTTCCAAGCG.
Then, sGFP fragment was amplified using the following primers:
sGFP_for-gBES1-C_F, 5′- GGTAAAGCTCATAGTATGGTGAGCAAGGGCG; and
sGFP_for-gBES1-C_R, 5′- GTATACTTCATCTCACTTGTACAGCTCGTCCATG.
The sGFP fragment was inserted just before the stop codon of BES1 by fusion of the two fragments using In-Fusion HD Cloning Kit. After sequence verification, this plasmid was introduced into the bes1-3 mutant.
For transactivation assay constructs, the NOS terminator was amplified by PCR using the following primers:
nosT-F, 5′-gccgcactcgagatatctagaATCGTTCAAACATTTGGCAA; and
nosT-R, 5′-tacaagaaagctgggtctagaGATCTAGTAACATAGATGAC.
The amplified fragment was cloned into the XbaI site of pENTR1A (no ccdB) using an In-Fusion HD Cloning Kit (TaKaRa). In the same way, the coding sequence of LUC+ was amplified and cloned into the XhoI-EcoRV site of the plasmid using the following primers:
LUC-XhoI-F, 5′-ttcgcggccgcactcgagATGGAAGACG; and
LUC-EcoRV-R, 5′-tgaacgattctagatatcTTACACGGCGATCTTTCCGC.
The resulting pENTR1A (LUC-nosT) was used for LHY::LUC. The promoter of LHY was amplified from Col-0 genomic DNA using the following primers:
LHY-pro-F (XhoI), 5′-cgcggccgcactcgaTTTTGGAATAATTTCGGTTATTTC; and
LHY-pro-R (XhoI), 5′-atccgcggatctcgaAACAGGACCGGTGCA.
The amplified fragment was cloned into the XhoI site of pENTR1A (LUC-nosT). The coding sequences of bes1-D(L) from cDNA of the bes1-D mutant and RLUC from the vector pRL vector (Promega) were amplified by PCR using the following primers:
bes1-D(L)-F, 5′- GGACTCTAGAGGATCATGAAAAGATTCTTCTATAATTCC;
bes1-D(L)-R, 5′- CGGTACCCGGGGATCTCAACTATGAGCTTTACCATTTCC;
RLUC-F, 5′- GGACTCTAGAGGATCATGACTTCGAAAGTTTATGATCC; and
RLUC-R, 5′-CGGTACCCGGGGATCTTATTGTTCATTTTTGAGAAC.
The amplified fragments were cloned into the BamHI site of pPZP211/NP/35S-nosT54 using an In-Fusion HD Cloning Kit.
Western blotting
Plants expressing BES1::BES1-GFP were harvested every 8 h from 24 h before induction in VISUAL, up to 72 h after induction. Approximately 50 mg of seedlings were ground into a fine powder in liquid nitrogen with a mortar and pestle, mixed with equal volumes of 2´Laemmli sample buffer (100 mM Tris-HCl(pH 6.8), 4%(w/v) SDS, 10%(v/v) 2-mercaptoethanol, and 20%(v/v) glycerol) and boiled at 95°C for 5 min. Samples were separated by SDS-PAGE on a 7.5% acrylamide gel, and transferred onto polyvinylidene fluoride membranes (Bio-Rad Laboratories). For the primary antibody, polyclonal anti-GFP (MBL-598) was diluted 1:2,000. For the secondary antibody, ECL Rabbit IgG, HRP-linked whole Ab (GE Healthcare) was diluted 1:10,000. Blots were visualized with ECL Prime reagent (GE Healthcare) and ImageQuant LAS 4010 (GE Healthcare).
ChIP-qPCR and ChIP-seq
Chromatin immunoprecipitation assays using BES1::BES1-GFP/bes1-3 and LUX::LUX-GFP/lux-4 were performed as described55 with modifications. Briefly, 600 mg of seedlings at 8 h after induction in VISUAL were fixed in PBS containing 1% paraformaldehyde for 10 min at room temperature, and nuclei and chromatin were isolated. The isolated chromatin was sheared with a Covaris S220 sonicator under these parameters: 4-6°C, 175 W peak power, 5% duty factor, 200 cycles/burst, for 50 s of treatment. To immunoprecipitate chromatin, 10 µL of anti-GFP antibody (MBL-598) and 50 µL of Dynabeads Protein G (Thermo Fisher Scientific) were used. The precipitated samples were subjected to qPCR or library preparation for ChIP-seq. MYB3056 was used as a positive control for ChIP with BES1. For ChIP-qPCR, specific sequences for each primer pair were:
LHYp-ChIP-1F, 5'- GATTCGGGTAGTTCAGTTCTTCG;
LHYp-ChIP-1R, 5'- GGTTAGTTCGGTTCGGTTCTAGG;
LHYp-ChIP-2F, 5'- CACCGTACCCACTTGTTTAGTCG;
LHYp-ChIP-2R, 5'- CGAGCCAGAAGCTTCAATGTG;
LHYp-ChIP-3F, 5'- GGCTCGTAGAGAAGCAACTTGAG;
LHYp-ChIP-3R, 5'- AGTCATCGCAGATCGACACG;
LHYp-ChIP-4F, 5'- GTGGATTCGTTTGGGTGAGG;
LHYp-ChIP-4R, 5'- AACAGTCGCTGCTTCTCCAG.
MYB30-ChIP-F, 5'-AGGTATTTTACGCTGGAAAATGTGT;
MYB30-ChIP-R, 5'- GAATCATCATAATAAGTATGGAGGTG;
ACT2-ChIP-F, 5'- CGTTTCGCTTTCCTTAGTGTTAGCT; and
ACT2-ChIP-R, 5'- AGCGAACGGATCTAGAGACTCACCTTG.
The following thermal cycling profile was used for all the primers: 95˚C for 60 s, ~40 cycles of 95˚C for 10 s, 60˚C for 45 s.
For ChIP-seq, the sequence libraries were prepared using a TruSeq ChIP Sample Preparation Kit v2 (Illumina), and sequenced using an Illumina NextSeq 500 system with a 75-nt single-end sequencing protocol. The sequence reads were mapped to the TAIR10 Arabidopsis genome sequence by HISAT257 with default parameters. Peaks were identified by MACS258, using the matching INPUT control with the genome size parameter “-g 1.3e8”.
Transactivation assay
Agrobacterium cultures carrying plasmids for the transactivation assay were grown overnight at 28°C, collected by centrifugation, and adjusted to an OD600 of 0.4 with infiltration buffer (10 mM MES(pH 5.6), 10 mM MgCl2, 150 μM acetosyringone, and 0.02% Silwet-L77). Cells were kept at 28°C in the dark for 1 to 2 h and then infiltrated into the abaxial air spaces of 4-week-old N. benthamiana plants grown under 12L/12D conditions at ZT0. After infiltration, plants were kept under 12L/12D conditions for 36 h and harvested at ZT12. The transactivation assay was performed with a Dual-Luciferase Reporter Assay System (Promega) according to the manufacturer’s instructions.
Data Availability
Code for PeakMatch is available online at https://github.com/endo-lab/PeakMatch. Other related, relevant lines and data supporting the findings of this study are available from the corresponding authors upon reasonable request.
Sequence data from this article can be found in The Arabidopsis Information Resource (TAIR) databases (https://www.arabidopsis.org) under the following accession numbers: UBQ14 (At4g02890), ACT2 (At3g18780), LHCB2.1 (At2g05100), CAB2 (At1g29920), CAB3 (At1g29910), AtHB8 (At4g32880), TDR (At5g61480), ANT (At4g37750), IRX3 (At5g17420), IRX8 (At5g54690), VND7 (At1g71930), SEOR1 (At3g01680), COR15A (At2g42540), ADH1 (At1g77120), RD29A (At5g52310), HAM2 (At3g60630), RGF6 (At4g16515), MYB30 (At3g28910), CCA1 (At2g46830), LHY (At1g01060), PRR3 (At5g60100), PRR5 (At5g24470), PRR7 (At5g02810), PRR9 (At2g46790), TOC1 (At5g61380), ELF3 (At2g25930), ELF4 (At2g40080), LUX/PCL1 (At3g46640), BOA/NOX (At5g59570), BES1 (At1g19350), CYCD3;1 (At4g34160), E2Fc (At1g47870), and RBR (At3g12280).