Light controls mesophyll-specific post-transcriptional splicing of photoregulatory genes by AtPRMT5

To explore the complex landscape of PTS events in response to light, we conducted a comprehensive study using Arabidopsis seedlings. We first grew seedlings in complete darkness for 5 d and then exposed these etiolated seedlings to white light for different durations: 1 h (DL1), or 6 h (DL6), using seedlings maintained in darkness as control (Dark condition). To identify PTS events, we extracted nuclear RNAs and constructed FLEP-seq2 libraries (Fig. 1A). In the nuclear fraction, we determined that about 20%, 15%, and 34% of all polyadenylated transcripts in the Dark, DL1, and DL6 libraries, respectively, retain unspliced introns (SI Appendix, Fig. S1).

To ensure accurate identification and quantification of PTS-associated introns, we also generated Illumina-based short-reads RNA-seq libraries from both nuclear and total RNA samples with three replicates, which showed strong reproducibility among replicates (SI Appendix, Table S1 and Fig. S2A). Importantly, we observed that the rate of IR in unspliced introns calculated by both Nanopore and Illumina RNA-seq is highly correlated (28) (SI Appendix, Fig. S2B).

We focused on PTS events during photomorphogenesis by comparing our sequencing data for the samples DL1 and DL6 to the Dark control samples, which identified 1,411 light-regulated PTS events (SI Appendix, Fig. S3) that can be classified into six clusters (Fig. 1 B and C). Importantly, these clusters followed different patterns of PTS behavior, including gradual increase (Cluster 3) (Dataset S1). Among these clusters, we noticed numerous PTS introns responding to light (Fig. 1D and SI Appendix, Fig. S4). For instance, the well-known regulator of flowering and stomatal opening SUPPRESSOR OF OVEREXPRESSION OF CO 1 (SOC1, also named AGAMOUS-LIKE 20 [AGL20]) (29–33) showed increased intron retention ratio within the nucleus upon exposure to light, indicating that SOC1 transcripts are post-transcriptionally spliced by light (Fig. 1E and SI Appendix, Fig. S5).

Interestingly, several genes exhibiting photodynamic PTS events have not previously been associated with photomorphogenesis. For example, the uncharacterized gene At5g14500 displayed a pattern opposite that of SOC1, with intron retention for its nuclear transcripts in the dark. After exposure to light, these intron-retaining transcripts were fully spliced and then exported to the cytoplasm (Fig. 1E and SI Appendix, Fig. S5). We speculate that these PTS events allow seedlings to rapidly respond to light signals during the transition from darkness to light.

To investigate whether light-induced PTS events exhibited a cell type–specific pattern, we used snRNA-seq, a method based on the 10 × Genomics Platform, to profile single-nucleus RNAs in Arabidopsis seedlings exposed to different light conditions. Following stringent filtering criteria (SI Appendix, Fig. S6 A and B), we obtained a dataset consisting of 17,947 (7,330 in Dark; 3,621 in DL1; and 6,996 in DL6) high-quality nuclei, containing 23,281 genes (SI Appendix, Table S2). To validate the effectiveness of our snRNA-seq, we performed a uniform manifold approximation and projection (UMAP) analysis (34–38) on the transcript count per nucleus (SI Appendix, Fig. S6 C–E), revealing that the majority of nuclei detected over 3,000 transcripts. We performed the UMAP analysis with top 30 principal components, estimated by JackStraw function (SI Appendix, Fig. S6F), to classify the nuclei into 17 distinct groups. With the help of marker genes for specific cell types, we annotated these 17 groups into 10 distinct cell types (Fig. 2 A and B and SI Appendix, Fig. S7) (39–50). The UMAP plot analysis demonstrated the presence of nuclei from all three conditions across all cell types following data integration. (Fig. 2C).

To uncover key regulators involved in the early light response, we identified 1,385 and 5,422 differentially expressed genes (DEGs) from comparisons between the Dark condition and the DL1 or DL6 condition, respectively (Dataset S2). Approximately 60% of these DEGs (810 out of 1,385 in DL1; and 3,193 out of 5,422 in DL6) were predominantly enriched in mesophyll cells relative to other cell types (Fig. 2D). Gene Ontology (GO) term enrichment analysis indicated that a significant proportion of these nuclear genes encode proteins that localize within chloroplast thylakoids and photosynthetic membranes (SI Appendix, Fig. S8), suggestive of the pivotal roles of light-responsive genes in mesophyll cells, particularly for chloroplast development and photosynthesis, following the transfer of etiolated seedlings from darkness into the light.

To assess the cell-type-specific characteristics of light-regulated PTS events, we quantified the enrichment levels of genes involved in light-associated PTS introns via snRNA-seq data (51). Intriguingly, we found that genes specifically in clusters 1, 5, and 6 (Fig. 1 B and C) were highly and significantly expressed within mesophyll cells under DL1 and DL6 conditions (Fig. 3A and SI Appendix, Fig. S9). Mesophyll cells represent a substantial portion of internal leaf tissue, contain abundant chloroplasts in light-grown seedlings, and play a crucial role in photosynthesis. To investigate the influence of de-etiolation on mesophyll cells, we analyzed 3,826 nuclei classified as mesophyll cells from the above snRNA-seq data with marker genes and performed a re-clustering analysis, clustering these nuclei into nine distinct groups (Dataset S3). We classified these nine groups of nuclei into three larger categories: Nuclei in groups A, E, and F specifically originated from Dark samples, while nuclei from groups C and I were of DL1 origin, with nuclei from groups B, D, G and H being of DL6 origin (Fig. 3B). This delineation highlights a developmental progression from Dark to DL6 conditions across these three larger categories. To investigate the developmental trajectory of mesophyll in response to light, we first identified differentially expressed genes for each group based on their distinct expression patterns (Dataset S4). For example, in dark category, we observed PHYTOCHROME A (phyA), a red and far-red light receptor (52–55), specifically expressed in A, E, and F groups (SI Appendix, Fig. S10). In addition, PHYTOCHROME INTERACTING FACTOR 3-LIKE 1 (PIL1) exhibits an expression pattern inversely linked to de-etiolation (56), confined to the E and F groups (SI Appendix, Fig. S10). In addition, a chlorophyll biosynthetic related gene, PROTOCHL OROPHYLLIDE OXIDOREDUCTASE A (PORA) (57, 58), only expressed in A group. Since exposed to light, several well-studied circadian oscillator–associated genes, such as CIRCADIAN CLOCK ASSOCIATED 1 (CCA1) (59–63), LATE ELONGATED HYPOCOTYL (LHY) (62, 63), and REVEILLE 8 (RVE8) (64), were induced in DL1 while dismissed in DL6 condition, suggesting their potential role in the early light signal response (SI Appendix, Fig. S10). Under DL6 condition, we observed genes involved in light signaling, such as B-box type zinc finger proteins, namely BBX14, BBX15 (65, 66), and SMALL AUXIN UPREGULATED RNA 14 (SAUR14) (67) were highly expressed (SI Appendix, Fig. S10). Additionally, amounts of genes associated with photosynthesis were also expressed in DL6, such as components of the photosynthesis NDH subcomplex B (PnsB3, PnsB4), and photosynthesis NDH luminal subcomplex 1 (PnsL1) (68) (SI Appendix, Fig. S10). More interestingly, we found chloroplast phosphate transporters, AtPHT4;1 (69) and AtPHT2;1 (70), along with downstream purple acid phosphatase PURPLE ACID PHOSPHATASE 3 (PAP3) were induced in DL6 (SI Appendix, Fig. S10). Phosphorus is a crucial element of ATP and NADPH, which plays a crucial role in photosynthesis (71). Phosphorus deficiency directly resulting in reduced photosynthesis efficiency (72). The upregulation of photosynthesis and phosphate homeostasis–related genes under DL6 conditions facilitates the transition from etiolation to growth in plants. These genes play a crucial role in sustaining photosynthetic efficiency, thereby generating energy for subsequent plant growth and development.

To investigate the developmental trajectory of mesophyll in response to light, we used the software Monocle to assign a pseudotime value (the nuclei’ progress in dark to light response) to each nucleus (Dataset S3) (73, 74). By mapping each nucleus onto this pseudotime trajectory, we reconstructed the developmental trajectory using all nuclei marked as originating from a mesophyll cell (Fig. 3 C and D). This pseudotime analysis revealed a predominant developmental trajectory for the mesophyll nuclei (87.7%, comprising 3,354 out of 3,826 cells). Along this trajectory, thousands of genes exhibited distinct expression patterns, of which 937 also exhibited PTS events (Dataset S5). Along the developmental trajectory, we observed distinct expression patterns for these 937 PTS genes. For example, genes involved in ABA signal transduction, such as ABA-RESPONSIVE ELEMENTS-BINDING FACTOR 3 (ABF3) (75) and HYPERSENSITIVE TO ABA1 (HAB1) (76, 77), remained constantly expressed. By contrast, development-related genes like URIDINE/CYTIDINE KINASE 3 (UCK3) (78) and genes involved in nutrient metabolism, such as ASPARAGINE SYNTHASE 1 (ASN1) (regulation of nitrogen reutilization) (79), and XYLOGLUCAN ENDOTRANSGLUCOSYLASE/HYDROLASE 15 (XTH15, involved in cell wall modification to promote shoot elongation) (80), initially showed an increase in their expression in the early stages (in the nuclei of mesophyll cells from etiolated seedlings), followed by a gradual decline after exposure to light. Additionally, genes such as EMBRYO DEFECTIVE 1211 (EMB1211), involved in embryo development and chloroplast biogenesis (81, 82) and RUBISCO ACTIVASE (RCA), encoding a chloroplast protein crucial for the light activation of Rubisco (83), were gradually up-regulated and their expression remained high following light exposure (Fig. 3E).

To investigate the regulatory mechanism underlying PTS regulation in light response, we focused on splicing-related factors, particularly AtPRMT5, known as a key factor in PTS (20). We found that AtPRMT5 affects 536 light-responsive PTS events (23) (SI Appendix, Fig. S11); in addition, atprmt5 mutants exhibited striking defects in cotyledon opening and delayed greening during de-etiolation (Fig. 4 A and C and SI Appendix, Fig. S12). Moreover, the cotyledon area of atprmt5 mutant seedlings was smaller than that of wild-type (Col-0) seedlings grown in continuous white light (Fig. 4 B and D). In particular, atprmt5 mutant seedlings had a longer hypocotyl when grown in the light but a shorter hypocotyl in continuous darkness compared to that of Col-0 seedlings (Fig. 4 E and F). These results indicate that AtPRMT5 is an essential factor that enables seedlings to properly respond to darkness and light conditions by orchestrating light-regulated PTS events, thereby facilitating prompt transition from skotomorphogenesis to photomorphogenesis.

Components of the light signaling pathway, such as the red-light photoreceptor phytochrome B (phyB) and COP1, are involved in light-controlled alternative splicing (AS) with various splicing factors (6–10, 12). Therefore, we integrated and analyzed intron retention events from various public RNA-seq datasets gathered from light signaling mutants, including mutants defective in phytochromes (84), cryptochromes (85), COP1 (86), DE-ETIOLATED 1 (DET1) (87), ELONGATED HYPOCOTYL 5 (HY5) (88), MYB DOMAIN PROTEIN 30 (MYB30) (89), PIFs (90), and B-BOX PROTEIN 11 (BBX 11) (91), to uncover potential regulators and regulatory networks. This analysis revealed that COP1 regulates a substantial number of light-associated PTS events (Fig. 5A). Unlike splicing factors, COP1 exerted both positive and negative roles in regulating IR on different transcripts. For example, several IR events observed in the cop1-4 mutant grown in continuous darkness coincided with those induced under DL6 conditions (121 genes). This correlation suggests a potential role for COP1 in facilitating the splicing of PTS introns in dark conditions (Fig. 5B and SI Appendix, Fig. S13). Of 273 PTS events regulated by COP1, 141 appeared to be co-regulated by AtPRMT5 in response to light (SI Appendix, Fig. S14). Furthermore, although the expression of COP1 and AtPRMT5 are not cell-type specific (SI Appendix, Fig. S15), PTS-associated transcripts regulated by both COP1 and AtPRMT5 were enriched within mesophyll cells (SI Appendix, Fig. S16), such as At2g44920 (a pentapeptide-repeat protein) (92) and protoporphyrinogen oxidase PROTOPORPHYRINOGEN OXIDASE 2 (PPO2) (93) (Fig. 5 C and D and SI Appendix, Fig. S17). We noted significant intron retention within the nuclei of cop1-4 and atprmt5-2 mutants, particularly in the darkness, for both At2g44920 and PPO2 genes, compared to the Col. This observation confirms that COP1 and AtPRMT5 regulate the post-transcriptional splicing of these genes within the nucleus under dark conditions (SI Appendix, Fig. S17). These results highlight the pivotal role of light-responsive PTS regulated by AtPRMT5 and COP1 in governing greening, photosynthesis, and morphogenesis during the transition from darkness to light.

Arabidopsis thaliana Columbia-0 (Col-0) ecotype was utilized as the wild type (WT). Mutant lines (atprmt5-1 and atprmt5-2) along with the complementary line 1 (COM-1) that have been previously characterized (21, 102, 103) were employed. These seeds were germinated on Murashige and Skoog (MS) medium adjusted to a pH of 5.8, supplemented with 1% sucrose and 0.8% agar. The growth protocol for Arabidopsis seedlings involved surface sterilization of the seeds, followed by a 3-d period of cold stratification in darkness at 4 °C. Subsequently, the seeds were cultivated on MS plates. To induce germination, a 14-h exposure to continuous white light was given. Then, the seedlings were allowed to grow for 5 d at 22 °C under either complete darkness or white light conditions. Notably, continuous darkness grown seedlings were exposed to white light spanning the 380 to 780 nm wavelength range for durations of either 1 or 6 h before sampling. For observation of hypocotyl length and cotyledon phenotypes, a light fluence rate of approximately 10 μmol·m⁻²·s⁻¹ was applied, while for other experiments, it was 45 μmol·m⁻²·s⁻¹.

To obtain complementary line 2 (COM-2), the genomic fragment of AtPRMT5 was subcloned into a binary vector that contained the sYFP-FALG tag. In detail, a 3′ UTR DNA fragment comprising 184 bp following the stop codon of AtPRMT5 and its stop codon was amplified using primers (SI Appendix, Table S3). This fragment was then inserted into the pCAMBIA1300-FLAG vector to generate pCAMBIA1300-FLAG-3′UTR. Subsequently, a 4,971 bp genomic fragment spanning 590 bp upstream of the start codon of AtPRMT5 and a 4,381 bp genomic region without the stop codon was amplified and cloned into pCAMBIA1300-FLAG-3′UTR, generating pCAMBIA1300-AtPRMT5p-AtPRMT5-FLAG-3′UTR. The sYFP segment, derived from pPLV18, was amplified and inserted into pCAMBIA1300-AtPRMT5p-AtPRMT5-FLAG-3′UTR to generate pCAMBIA1300-AtPRMT5p-AtPRMT5-sYFP-FLAG-3′UTR. Validation of all constructs was performed through DNA sequencing, and the primer details are provided in SI Appendix, Table S3. Subsequently, the plasmid was introduced into the EHA105 strain of Agrobacterium tumefaciens and transformed into atprmt5-2 plants to generate complemented line AtPRMT5-sYFP-FALG.

The quantification of chlorophyll content followed previously established methods (104). Briefly, Arabidopsis seedlings, having grown for 5 d in continuous darkness followed by exposure to light conditions for 24 h, were weighed (approximately 0.05 to 0.07 g) and ground in liquid nitrogen within a darkened room under dim green light. After that, chlorophyll was extracted from the samples using 80% acetone in water and incubated at 4 °C for 3 h to ensure complete dissolution. The concentration of chlorophyll was calculated by measuring the absorption at 663 and 645 nm.

The extraction of RNA from different cellular fractions following the previously published methods (26) with minor modifications. The frozen materials were ground with liquid nitrogen to fine powders using a mortar and pestle. For total RNA extraction, appropriate powders were utilized followed by RNA extraction. For nuclear fraction isolation, the rest powders were incubated with 10 mL Honda buffer [0.44 M sucrose (Sigma), 1.25% (w/v) Ficoll (Sigma), 2.5% (w/v) dextran T40 (Macklin), 20 mM HEPES-KOH pH 7.4, 10 mM MgCl₂, 0.5% (w/v) Triton X-100 (Sigma), 1 mM dithiothreitol (DTT, Thermo Fisher), 1× protease inhibitor (Sigma), and 100 ng μl⁻¹ tRNA (Mei5bio)], followed by homogenization through vortex. To remove tissue debris, the homogenate was filtered through a layer of Miracloth (Millpore). Then, the nuclei fraction was pelleted by centrifugation at 3,500×g and 4 °C for 5 min, and the pellet was resuspended in 15 mL Honda buffer for two rounds of wash. Next, the pellet was resuspended in 1 mL Honda buffer supplemented centrifugated at 4 °C and 8,000×g for 1 min. The supernatant was discarded, and the pellet was weighted, which was followed by nuclear RNA extraction. Total RNA and nuclear RNA were extracted using RNAprep Pure Plant Kit (Polysaccharides & Polyphenolics-rich) (Tiangen).

ANOVA analyses were performed using SPSS statistical software, and Student’s t tests were performed in Microsoft Excel. Distinct letters indicate statistical significance determined by ANOVA (P-value < 0.05) for multiple comparisons. The same letter denotes levels that do not show significant differences.

The mRNA libraries preparation and sequencing were performed at the Beijing Genomics Institute (BGI). In brief, polyadenylated RNAs from nuclear and total fractions were enriched using Oligo-dT beads and subsequently reverse-transcribed with N6 random primers. This was followed by adaptor ligation and synthesis of double-stranded DNA (dsDNA). The dsDNA was then denatured, leading to the production of single-stranded DNA which underwent circularization. DNA nanoballs were subsequently generated via rolling circle amplification. Then, the synthesized DNA nanoballs were subjected to strand-specific mRNA sequencing (PE150) on a DNBSEQ platform. Processing of NGS data were performed as previously described (20) with minor modifications. Initially, paired-end sequencing reads were aligned to Arabidopsis TAIR10 reference genome (105) using Hisat2 (v2.1.0) (106) with the parameters “--min-intronlen 20 --max_intronlen 12000.” The reads aligned to rDNA, mitochondria and chloroplast genomes were excluded by a custom Python script “filter_rRNA_bam.py” (20). Then, the PCR duplicates were removed using Picard with default paraments, and the FPKM was calculated by Stringtie with the paraments “-e --rf -B” using the Araport11 annotation file (107). Next, to calculate the percentage of intron retention (PIR) value for each intron, a customized Python script “ASCaller.py” was employed (20). Then, DEseq2 (v.1.19.31) (108) was used to conduct differential IR test between different samples, and the resulting P-values were adjusted by the Benjamini–Hochberg method to control the False Discovery Rate (FDR). Introns with IR ratio > 0.1 and P-value < 0.05 were identified as introns with differential retention ratio. Finally, the light-regulated pts-introns were identified by intersecting a pts-introns list identified in a previous study (20) with the differential retained intron list and further grouped into six clusters using K-means in R.

The library construction and sequencing procedures were conducted as previously reported with minor modifications (26). In brief, 500 to 1,000 ng nuclear RNA was used for library construction followed by rRNA removal. Then, a 3′ linker (NEB, 5′ rAppCTGTAGGCACCATCAAT–NH2 3′) was ligated to the 3′ terminus of purified rRNA-depleted RNA. The subsequent library construction was carried out following the PCR-cDNA Barcoding Kit (SQK-PCB109) and sequenced on a MinION device using an R9.4 flow cell (Oxford Nanopore Technologies). The nanopore data processing and calculation of the PIR value of each intron was performed following the FLEP-seq pipeline (https://github.com/ZhaiLab-SUSTech/FLEPSeq) (20).

To convert Nanopore raw signal to fasta sequence, Guppy basecaller (v4.0.11) was used with parameters “--c dna_ r9.4.1_450bps_hac.cfg, --qscore_filtering, --device ‘cuda:all:100%’ --barcode_kits ‘SQK-PCB109’.” The basecalled reads were mapped to TAIR10 reference genome using Minimap2 (v2.10-r761) (109) with the parameters “-ax splice, –secondary = no, -G 12000.” Then, the 3′ linker was identified by a customized Python script “adapterFinder.py” with parameter “--mode 1” as previously described (26). Only the full-length polyadenylated reads were used for PIR value calculation by using a customized R script “cal_polya_transcript_ir.R” (20).

The nuclei isolation procedures were performed following a previously reported method with some modifications (110). Initially, nuclei isolation buffer (1×NIB, MilliporeSigma) was freshly prepared on ice, supplemented with 1 mM DTT (Thermo Fisher), 1× protease inhibitor (Sigma), and 0.4 U μL⁻¹ murine RNase inhibitor (Vazyme). All subsequent procedures were conducted on ice or at 4 °C. The frozen seedlings were chopped in 500 to 1,000 μL NIB with a razor blade for 5 min to release nuclei and incubated with 10 mL NIB. The nuclei extract was subjected to a brief centrifugation at 100×g for 1 min to pellet tissue debris. Then, the supernatant was filtered through a 40-μm strainer (Sigma) and centrifuged at 1,000×g for 5 min. Next, the pelleted nuclei were gently resuspended in an appropriate volume of NIB by pipetting and filtered through a 20-μm strainer (pluriStrainer Mini) to eliminate remaining debris. For the isolation of nuclei, frozen materials from 5-d-old Arabidopsis seedlings were utilized. For cell sorting, the nuclei suspension was stained with DAPI and loaded to a flow cytometer (SONY MA900) equipped with a 100-μm nozzle. Approximately 100,000 nuclei were sorted for each sample based on DAPI signal intensity and nuclear size. Then, the sorted nuclei were collected into a 2-mL collection buffer composed of 1× PBS (CORNING), 1% BSA (Sigma), and 0.4 U μL⁻¹ murine RNase inhibitor (Vazyme). The sorted nuclei were centrifuged at 1,000×g for 5 min and resuspended in 50 μL collection buffer. After checking the quality and number of nuclei under microscopy manually, ~10,000 nuclei were loaded onto 10× Genomics Chip G and processed by 10× Genomics Controller for nuclei barcoding. The library construction was carried out following the manufacturer's instructions of 10× Chromium Single Cell 3′ Solution v3.1 kit. Each library was sequenced with Illumina NovaSeq 6000.

Initially, the raw reads were processed using Cell Ranger (v6.1.2) (111), and aligned to Arabidopsis TAIR10 genome under default settings, generating h5ad files for each sample. The outputs were analyzed using the Python package SCANPY (v1.9.0) (112). Next, we used ScDblFinder (v1.10.0) (113) to remove the doublets, followed by detailed quality control procedures for each sample. High-quality nuclei from three libraries were kept with a cutoff that has genes 400 to 4,000, UMI 400 to 7,500, and containing transcripts from mitochondria or chloroplasts less than 10%. After the quality control step, the count matrices from the three samples were integrated scVI (v0.16.0) (114), in accordance with the official tutorial, except that “n_top_genes” was set to 5,000. Data normalization was achieved using the “scanpy.pp.normalize_total” function. Then, the integrated data were subjected to “scanpy.pp.neighbors” function to construct the nearest-neighbor graph with a parameter “n_pcs=30” (the number of principal components, estimated by JackStraw function, SI Appendix, Fig. S6F). Next, cell clustering utilized the Leiden algorithm (115) via the “scanpy.tl.leiden” function (with parameters “resolution= 0.5”). Data visualization was facilitated through the Uniform Manifold Approximation and Projection (UMAP), employed through “scanpy.tl.umap” function (with parameters “min_dist = 0.2”). Cluster-specific genes were identified using the CELLEX (v1.2.2) algorithm (116), applying the specificity score criterion of greater than 0.9. For functional analysis of cluster-specific genes, an R package clusterProfiler (v4.6.0) (117) was used. The DEGs of each group were identified using DESeq2 method with default parameters from Delegate, the number of pseudo-replicates was set to 3, and filtered with FDR < 0.05. Finally, the enrichment level of a gene set was evaluated by area under the curve (AUC) score calculated by AUCell package in R.

The gene expression matrix of mesophyll cells was extracted from the object created by Scanpy and underwent further analysis using the Seurat (v4.3.0.1) (118) and Monocle (v2.28.0) (73) pipelines. To begin, this expression matrix was initially loaded to create a Seurat object using the CreateSeuratObject function with default parameters. The JackStraw function was then employed to estimate the optimal number of dimensions for subsequent UMAP algorithm implementation (RunUMAP function with patameter dims = 1:20). Given that we have successfully clustered the mesophyll cells into eight distinct groups, we proceeded with the Monocle pipeline, which loaded metadata and marker genes for each group to sort the cells based on their pseudotime values, facilitated by the setOrderingFilter function. Finally, we visualized the developmental trajectories of mesophyll cells using the plot_cell_trajectory function.

Y.Y., H.L., Y.Q., J.Z., and X.C. designed research; Y.Y., Y.Q., T.Y., Y.H., and Z.L. performed research; J.J. and Y.L. contributed new reagents/analytic tools; Y.Y., H.L., and Y.Q. analyzed data; and Y.Y., H.L., Y.Q., J.Z., Y.L., X.D., and X.C. wrote the paper.

J.Z. has co-published a paper with Dr. Robert Schmitz on December of 2019 (doi: https://doi.org/10.1038/s41477-019-0547-0) with no direct collaboration.