Weeping candidate genes screened using comparative transcriptomic analysis of weeping and upright progeny in an F1 population of Prunus mume.

Weeping is a specific plant architecture with high ornamental value. Despite the considerable importance of the weeping habit to landscaping applications and knowledge of plant architecture biology, little is known regarding the underlying molecular mechanisms. In this study, growth and phytohormone content were analyzed among the progeny of different branch types in an F1 mapping population of Prunus mume with varying plant architecture. Bulked segregant RNA sequencing was conducted to compare differences among progeny at a transcriptional level. The weeping habit appears to be a complex process regulated by a series of metabolic pathways, with photosynthesis and flavonoid biosynthesis highly enriched in differentially expressed genes between weeping and upright progeny. Based on functional annotation and homologous analyses, we identified 30 candidate genes related to weeping that merit further analysis, including 10 genes related to IAA and GA3 biosynthesis, together with 6 genes related to secondary branch growth. The results of this study will facilitate further studies of the molecular mechanisms underlying the weeping habit in P. mume.

4‐coumarate:CoA ligase

abscisic acid

bulked segregant analysis

bulked segregant RNA‐Seq

clusters of orthologous groups of proteins

ent‐copalyl diphosphate synthase

cellulose synthase like D2

CYTOCHROME P450 FAMILY 79 SUBFAMILY B POLYPEPTIDE 2

cis‐zeatin

differentially expressed genes

dihydrozeatin

gibberellin 20‐oxidase

gibberellin

GIBBERELLIN INSENSITIVE DWARF

gene ontology

indole‐3‐acetic acid

3‐indolebutyric acid

indole‐3‐carboxaldehyde

N6‐isopentenyladenine

ent‐kaurene synthase

longly guy

methyl indole‐3‐acetate

phenylalanine ammonia lyase

PIN‐FORMED

peroxidase

TRYPTOPHAN AMINOTRANSFERASE OF ARABIDOPSIS

TILLER ANGLE CONTROL

transcription factor

trans‐zeatin

YUCCA6

Introduction

Prunus mume (mei) is a traditional ornamental tree that has been cultivated for more than 3000 years in China. Mei has many valuable ornamental traits, including pleasant floral scent and weeping habit. The weeping habit in woody species represents a unique form of tree architecture that has received increasing attention during recent years (Zhang et al. 2015); however, few systematic studies have examined its molecular and physiological aspects. The weeping habit is generally attributed to a lack of branch structural integrity and insufficient tensile strength to support upright branch orientation. The development of weeping branches in mei and peach (P. persica) may be related to differential secondary branch growth characteristics, including secondary cell wall thickening and branch lignification. Secondary growth involves several biological processes including cell division and expansion, secondary cell wall biosynthesis, lignification and programmed cell death (Oh et al. 2003). The study of weeping mei architecture has focused on processes related to secondary cell wall biosynthesis and lignification. In mei, secondary cell walls are thinner and develop more slowly in the elongation zone of weeping branches than in the elongation zone of upright branches; furthermore, lignification is slower in mei weeping branches (Zhang 2016, Hou et al. 2017). The process of secondary wall thickening is accompanied by cellulose and hemicellulose biosynthesis, whereas lignification is associated with lignin biosynthesis and accumulation. Previous studies by our research group revealed differences in cellulose, hemicellulose, and lignin content between weeping and upright mei branches (Zhang 2016, Hou et al. 2017). We also found lower lignin content on the adaxial side of upright branches than on the abaxial side, whereas the opposite trend was observed in weeping branches; cellulose content distribution did not differ between mei branch types. Thus, the biosynthesis of cellulose, hemicellulose and lignin were related to weeping branch architecture in mei.

Differences in phytohormone content have also been observed between weeping and upright mei branches. Differential distribution of auxin (indole‐3‐acetic acid, IAA) and gibberellin (GA3) within branches appear to affect the development of weeping branches in mei (Wang 2014). The IAA and GA content of base part of upright branch in the adaxial side was lower than those on the abaxial side, whereas in the base part of weeping branch only the content of GA was higher in the adaxial side than the abaxial side (Wang 2014). GAs are among the key phytohormones involved in formation of the weeping habit in peach; non‐uniform distribution of IAA due to polar transport promotes upright branch growth in rice (Oryza sativa; Wang et al. 2018). Hormone‐related genes were found to be differentially expressed between two willow (Salix matsudana) varieties with different branch types. In such species, genes associated with IAA and GA are highly likely to be responsible for the weeping habit (Liu et al. 2017). Both IAA and GA3 appear to play important roles in branch elongation among many angiosperms (Tanimoto 2005), which is consistent with the longer internode length of the weeping peach shoot, compared with the upright shoot. Further studies are required to comprehensively evaluate the roles of phytohormones in formation of the weeping habit. In addition to IAA and GA3, two other important phytohormones – cytokinin and abscisic acid (ABA) – should be further studied.

The exact molecular mechanism and regulatory pathway responsible for the weeping habit in plants are presently unclear. Genes controlling the weeping habit have been identified in a few species; most reports thus far have focused on genes related to phytohormones and secondary development based on physiological studies (Sugano et al. 2004, Joshi et al. 2011, Wang 2014, Zhang et al. 2018). Several key genes in the GA and IAA biosynthesis pathways have been cloned; relationships of these genes with the weeping habit have been explored by comparing differential gene expression levels between weeping and upright branches. The expression of Ps3ox, a key gene in the GA biosynthesis pathway, was upregulated in weeping shoots, compared with upright shoots, in weeping Japanese flowering cherry (P. spachiana)at different development stages. Two GA biosynthesis genes [P. mume gibberellin 20‐oxidase (PmGA20ox) and PmGA3ox] and three IAA transport genes [PmAUX1, P. mume PIN‐FORMED 1 (PmPIN1), and PmPIN3] have also been cloned in mei. The expression patterns of these genes have been found to differ significantly between weeping and upright shoots (Wang 2014). LAZY, a member of IGT family which is conserved in GφL(A/T)IGT motif, was reported to promote asymmetric IAA localization in rice (Dardick et al. 2013, Zhang et al. 2018); TILLER ANGLE CONTROL (TAC), another IGT family member, was suggested to act upstream of WEEP together with LAZY. Overexpression of a secondary wall‐associated cellulose synthase gene (PtdCesA8) leads to the weeping habit in aspen (Populus tremuloides; Joshi et al. 2011). Recently, WEEP was shown to cause the weeping peach architecture; this was the first report of a gene verified to partially control the weeping habit (Hollender et al. 2018). The inheritance mode of the weeping habit has been studied by cross‐breeding in various plants, especially woody plants (Zhang 2016). In peach and mulberry (Morus alba), the weeping habit is likely controlled by a single recessive gene, according to early cross‐breeding results (Bassi and Rizzo 2000). A mutant allele frequency analysis and density mapping showed that four genomic regions were significantly associated with weeping in Malus domestica (Dougherty et al. 2018).

Thus, differences of secondary development and plant hormone content among branches are known to be related to the development of the weeping habit. However, the genes responsible for these physiological differences and ultimate development of the weeping habit remain unidentified. The weeping habit is a complex phenotype that involves complex biological processes. Our research group has conducted mei crossbreeding for many years; in a previous study, an F1 mapping population with variation in plant architecture was generated using female ‘LiuBan' (upright type) and male ‘FenTaiChuiZhi’ (weeping type) mei (Zhang et al. 2015). Based on χ2 test results for hybrid progeny segregation ratios in different hybrid combinations, the weeping habit was found to be controlled by a major recessive gene, together with some minor loci (Zhang et al. 2015). The major locus for the weeping habit was located on a 1.14‐Mb interval, with 159 predicted protein encoding genes on pseudochromosome 7 (Zhang et al. 2015). Thus far, there have been few physiological and molecular studies of the weeping habit; therefore, it is difficult to identify the key candidate gene among the 159 identified genes. Several studies have combined the power of bulk segregant analysis (BSA) with the ease of RNA sequencing (RNA‐Seq) into a new strategy, BSR‐Seq; this strategy is both inexpensive and efficient for large genomes due to its use of high throughput for deep coverage of the transcriptome and its strong ability to detect the genetic differences underlying traits (Wang et al. 2013a). BSR‐Seq has been successfully applied to gene cloning and molecular breeding of important phenotypes of economic crops, such as maize (Zea mays) and wheat (Triticum aestivum; Liu et al. 2012, Ramirez‐Gonzalez et al. 2015). Thus, BSR‐Seq can be used to further screen candidate genes for the weeping habit. In this study, we detected differences in branch growth characteristics and phytohormones between different sides of two branch types in a mei F1 population. We then applied BSR‐Seq to screen candidate genes related to the weeping habit in mei. The transcriptional‐level expression information regarding candidate genes determined in this study will provide a reference for the establishment of a gene regulation network for the weeping habit in mei.

Materials and methods

Sampling, RNA isolation and Illumina sequencing

The F1 population was generated by the female ‘LiuBan' (upright branch type) from QingDao, China (36.20169°N, 120.4162°E) and the male ‘FenTaiChuiZhi’ (weeping branch type) from WuHan, China (30.54526°N, 114.39511°E) in 2011. In spring 2013, a clonally replicated plantation of this entire pedigree was established in a randomized complete block design with three replicates and two‐tree plots at Hangzhou, China (30.566389°N, 119.879582°E; Zhang et al. 2015). Fifty weeping progenies and 50 upright progenies with similar growth condition were selected to collect annual shoot before lignification in 2018 [the fourth stage of branch development as described by Zhang (2016), i.e. the leaves on the new shoot are almost fully spread with the length of the new shoot being 10–14 cm], and then quickly divided into adaxial and abaxial sides along the central axis using a knife (Fig. 1). 0.3 cm width ×9 cm length of branch tissue of the adaxial and abaxial sides of shoots of 50 weeping progenies and 50 upright progenies were separately collected and mixed, and transferred to liquid nitrogen immediately for preservation (Supporting Information, Fig. S1). Equal amount of branch tissue of the adaxial and abaxial sides of shoots were mixed separately to form four sample pools with three replicates. Namely, the adaxial side pool of 50 weeping progenies (WAd), the abaxial side pool of 50 weeping progenies (WAb), the adaxial side pool of upright progenies (UAd) and the abaxial side pool of upright progenies (UAb). Then, total RNA was isolated using an RNA extraction kit (Aidlab) according to the manufacturer's instructions. The RNA integrity number (RIN) was greater than 9.0 for all samples. One microgram RNA of each sample was used for sequencing libraries construction, with the RNA sequencing was performed on Illumina HiSeq 2000 platform (Majorbio, Shanghai, China).

Fig. 1

The phenotypes of Prunus mume (mei) upright (A) and weeping (B) progenies and the growth characteristics of annual branches.

Measurement of annual branch length, plant height and hormone content in weeping and upright progenies of the ‘LiuBan'×’FenTaiChuiZhi’ F1 population

Thirty annual branches length of 25 weeping progenies and 25 upright progenies with similar growth condition were selected and measured after they stopped growth in November 2018. The difference of the annual branch‐length between the weeping and upright branches was calculated by IBM SPSS statistics 21 program. Graphs were generated by Excel 2016.

Relative quantities of 10 plant hormone [IAA, methyl indole‐3‐acetate (ME‐IAA), 3‐indolebutyric acid (IBA), indole‐3‐carboxaldehyde (ICA), N6‐isopentenyladenine (IP), cis‐zeatin (cZ), trans‐zeatin (tZ), and dihydrozeatin (DZ), GA3 and ABA] in the WAd, WAb, UAd and UAb pools which were consistent with the pools of RNA‐Seq were analyzed using a liquid chromatography–tandem mass spectrometry (LC–MS/MS) system comprising the Shim‐pack UFLC SHIMADZU CBM30A system (Shimadzu) and the 6500 Quadrupole TRAP system (Applied Biosystems). The samples were first freeze‐dried and then grounded to powder using a grinder (MM 400; Retsch) at 30 Hz for 1.5 min. Then, 120 mg of the powder was placed in 1.2 ml of 80% aqueous methanol at 4°C overnight, and vortexed for six times (every half an hour) to increase extraction efficiency. Next, the extract was centrifuged at 12 000 g for 15 min, and the supernatant was obtained which was then dried at 35°C using liquid nitrogen. Samples were placed in 100 μl 30% aqueous methanol and swirled to fully dissolved again, then centrifuged at 12 000 g for 15 min to obtain the supernatant. All the samples were stored in an injection bottle before LC–MS/MS analysis. The analytical conditions were as follows: LC: column, Waters ACQUITY UPLC HSS T3 C18 (1.8 μm, 2.1  × 100 mm); solvent system, water (0.04% acetic acid), acetonitrile (0.04% acetic acid); gradient program, 95:5 (v/v) at 0 min, 5:95 (v/v) at 11.0 min, 5:95 (v/v) at 12.0 min, 95:5 (v/v) at 12.1 min, 95:5 (v/v) at 15.0 min; flow rate, 0.35 ml min−1; temperature, 40°C; injection volume, 5 μl.

Data filtering, de novo assembly, SNP identification and functional annotation analysis

The raw sequencing reads were trimmed to remove the adaptor sequences, duplicated sequences, ambiguous reads (‘N'), and low‐quality reads using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle). The reference genome of Mei and corresponding gene annotation files were downloaded from http://prunusmumegenome.bjfu.edu.cn (Zhang et al. 2012). Then, clean reads were used for sequence alignment against reference genome of Mei (Zhang et al. 2012) using the TopHat2 (http://tophat.cbcb.umd.edu/; Kim et al. 2013) and de novo transcriptome assembly was carried out using Cufflinks (http://cole-trapnelllab.github.io/cufflinks/; Trapnell et al. 2010) to obtain unigenes. Meanwhile, the gene mapping rates (the ratio of reads which were mapped to a corresponding gene in the gene annotation file) were calculated. Samtools (Version 1.9; https://github.com/samtools/samtools.git) and GATK (Version 3.8; https://software.broadinstitute. org/gatk/download/; McKenna et al. 2010) were used to perform SNP‐calling and InDel‐calling, respectively. The difference in SNP allele frequencies between the weeping and upright pools was evaluated using SNP‐index methods (Takagi et al. 2013). The difference in proportions of one genotype in weeping pool was compared with upright pool and was defined as Δ(SNP‐index). And Δ(SNP‐index) was calculated for all the SNP positions. The expression level of each unigene was measured as FPKM (fragments per kilobase per million reads) by RSEM (RNA‐Seq by Expectation Maximization) software (Li and Dewey 2011). Annotation of all unigenes were carried out using BLASTX (http://blast.ncbi.nim.nih.gov/Blast.cgi/; E‐value < 1e−5) against NR (non‐redundant protein sequence in NCBI, http://www.ncbi.nlm.nih.gov/), Swiss‐Prot (http://www.expasy.ch/sprot/), Pfam (http://pfam.xfam.org/), COG (Clusters of Orthologous Groups of proteins, http://www.ncbi.nim.nih.gov/COG/), GO (Gene Ontology, http://geneontology.org/) and KEGG (Kyoto Encyclopedia of Genes and Genomes, https://www.genome.jp/kegg/) databases.

Detection and functional enrichment analyses of differentially expressed genes (DEGs)

Unigenes with the false discovery rate (FDR)/p adjusted value (padj) ≤ 0.05 and the FC ≥1.5 between any two average libraries (WAd vs UAd, WAb vs UAb, WAd vs WAb or UAd vs UAb) were considered as differentially expressed gene (DEG). And the DESeq2 implemented in R package was used to identify all the DEGs. GO enrichment analysis of DEGs was further implemented by the Goatools (Faria‐Blanc et al. 2018). The KEGG pathway enrichment analysis of DEGs was performed using KOBAS (KEGG Orthology Based Annotation System, Peking University, Beijing, China) software (Xie et al. 2011).

Phylogenetic analyses of DEGs and the well‐characterized homologous genes

To further clarify the specific function of DEGs, genes which had been verified to be plant hormone‐related or secondary‐growth‐related in Arabidopsis thaliana (Arabidopsis) were obtained from the TAIR website (https://www.arabidopsis.org/) for phylogenetic analysis. Then, multiple sequence alignment was performed using MUSCLE implemented in MEGA 6.0 and the phylogenetic tree was inferred using Neighbor‐Joining method with 1000 bootstrap in MEGA 6.0 (Tamura et al. 2013).

Gene expression analysis by qRT‐PCR

In order to verify the reliability of the gene expression results of the RNA‐seq, qRT‐PCR experiment was conducted for nine candidate genes. Total RNA was isolated using an RNA extraction kit (Aidlab) according to the manufacturer's instructions. Two thousandng total RNA was used for reverse transcription using 5X All‐In‐One RT Mastermix (Applied Biological Materials). Synthetic first‐strand cDNAs were then diluted 10‐fold for gene expression analyses. Primers for quantitative real‐time PCR (qRT‐PCR) were designed using Primer Premier 5 software (Lalitha 2000; Table S9). qRT‐PCR was conducted using the LightCycler/LightCycler96 System Real Time PCR (Roche) with SYBR Premix Ex Taq II (TaKaRa). Three biological replicates and three technical repeats were contained for each gene and sample. Relative gene expression was normalized by comparison with the expression of the ACTIN gene and ER1‐α was analyzed using the 2–ΔΔCT method (Livak and Schmittgen 2001).

Results

Variations in physiological characteristics of annual branches in weeping and upright progeny of the ‘LiuBan' × ‘FenTaiChuiZhi’ F1 population

As shown in Fig. 1A,B, annual branches of upright progeny of the ‘LiuBan' × ‘FenTaiChuiZhi’ F1 population grew upward, whereas those of the weeping progeny were naturally pendulous. To further characterize physiological differences between weeping and upright progeny of the ‘LiuBan' × ‘FenTaiChuiZhi’ F1 population, the lengths of 30 annual branches from 25 weeping and 25 upright progenies under similar growing conditions were selected and measured. The longest average annual branch lengths of weeping and upright progeny were 151.58 and 57.65 cm (Table S1), respectively. Among the upright progeny, the longest branch was 64.94 cm, while the shortest was 44.83 cm. Annual branch length was significantly greater in weeping progeny than in upright progeny (P < 0.01; Table S1, Fig. 2A). Plant heights were then measured in 25 weeping and 25 upright progeny under similar growing conditions to detect differences between progeny with different plant architectures. Plant height was highly significantly greater in upright progeny than in weeping progeny (P < 0.01; Fig. 2B).

Fig. 2

Boxplot of the length of annual branches (A) and plant height (B) of 25 weeping progenies and 25 upright progenies. (A) The branche length are presented as the means ± sd of the lengths of 30 annual branches from 25 weeping and 25 upright progenies. **P < 0.01, t‐test. (B) The plant height are presented as the means ± sd of 25 weeping and 25 upright progenies. **P < 0.01, t‐test.

Phytohormone content was detected in pooled tissues. In total, four auxins (IAA, ME‐IAA, IBA and ICA), four cytokinins (IP, cZ, tZ and DZ), one gibberellin (GA3), and ABA were detected. Among the four auxins, only IAA and ICA contents significantly differed between weeping and upright progeny (both P < 0.01; Fig. 3). IAA contents in the adaxial (WAd pool) and abaxial (WAb pool) sides of weeping progeny were significantly greater than in upright progeny. Among plants with the same branch architecture type, IAA content on the abaxial side of upright branches (UAb pool) was significantly greater than that on the adaxial side (UAd pool); we could not detect a statistically significant difference between WAd and WAb pools (Fig. 3). ICA content was also significantly greater in the WAb pool than in the UAb pool. Significant differences in cytokinin (IP, tZ, cZ, and DZ) contents were observed between different branch types on the same side of the branch (WAd vs UAd, WAb vs UAb), except for the DZ content between the WAb and UAb pools (all P < 0.01). In detail, IP, tZ, and cZ contents were significantly greater in the UAb pool than in the WAb pool; IP, tZ, cZ and DZ contents significantly (P < 0.05) differed between different sides of weeping branches. GA3 content differed significantly between WAd and UAd, WAb and UAb, WAd and WAb, and UAd and UAb (all P < 0.01). Specifically, GA3 content was significantly greater in weeping progeny pools than in upright progeny pools; it was also significantly greater on the adaxial side than on the abaxial side in both branch types (Fig. 3). ABA distributions in both branch types were similar, with significantly (P < 0.01) greater ABA content on the abaxial side than on the adaxial side.

Fig. 3

Phytohormone content diagram of the adaxial and abaxial sides of mei annual weeping and upright branch pools. *P < 0.05, **P < 0.01, t‐test, n = 3.

De novo transcriptome assembly, single‐nucleotide polymorphism (SNP) identification, and functional annotation of pools of weeping and upright progeny of the F1 population

To further elucidate differences between weeping and upright mei progeny at the transcriptional level, 12 high‐quality libraries (Q20 and Q30 values >90%) were constructed for the WAd, WAb, UAd and UAb pools with three biological replicates. In total, > 670 million raw reads were acquired from the 12 libraries with 49–64 million reads per library (Table S2). After removal of low‐quality reads, adaptors, ambiguous reads, and duplicates, an average of 99% clean reads were obtained for each library (Table S2). The clean reads were aligned to the mei reference genome (Zhang et al. 2012). As shown in Table S2, >90% of the clean reads were mapped onto the mei reference genome; >88% of the clean reads were uniquely mapped, demonstrating the high quality and confidence level of the sequencing results. In total, 31 390 unigenes were obtained, including 32 787 SNPs and 477 indels between weeping and upright progeny. Pseudo‐chromosome 2 (Pm2) contained the largest number of SNPs (6473) and indels (95); 2311 SNPs and 47 indels were detected on Pm8 (Table 1). The ΔSNP index was then calculated between the weeping and upright progeny pools; Pm7 was found to contain the most SNPs, with high ΔSNP index values (≥0.8). For comprehensive unigene annotation, all unigenes were aligned to six public databases: non‐redundant (NR), Swiss‐Prot, Pfam, Clusters of Orthologous Groups (COG), Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). Among all unigenes, 17 350 (55.27%) were annotated in the GO database, 11 214 (35.72%) were annotated in KEGG, 25 048 (79.80%) were annotated in COG, 29 191 (92.99%) were annotated in NR, 20 887 (66.54%) were annotated in Swiss‐Prot, and 22 234 (70.83%) were annotated in Pfam. In total, 29 323 (93.42%) unigenes were annotated from at least one of these databases. Based on gene annotation, variations with high ΔSNP index values on Pm7 were distributed in 42 genes.

Table 1

Statistics of SNP and InDel between weeping pools and upright pools on each pseudo‐chromosomes of mei.

Chromosome	Total		Δ(SNP‐index) ≥ 0.8
	SNP	InDel
Pm1	3965	64	2
Pm2	6473	95	28
Pm3	2749	31	3
Pm4	3589	43	2
Pm5	3568	55	1
Pm6	3408	54	5
Pm7	3574	43	65
Pm8	2311	47	2
Scaffolds	3150	45	12

Functional analysis of differentially expressed genes (DEGs) among mei F1 progeny pools

In total, 2208 DEGs were obtained from comparisons between WAd and UAd, WAb and UAb, WAd and WAb, and UAd and UAb, which yielded 1782, 845, 173, and 3 DEGs, respectively (Fig. 4E). Among these, WAd vs UAd and WAb vs UAb produced the largest number of DEGs. In total, 1038 and 744 genes were upregulated and downregulated, respectively, in the WAd vs UAd pool (Fig. 4A); 373 and 472 genes were upregulated and downregulated in the WAb vs UAb pool, respectively (Fig. 4B). Few DEGs were detected between WAd and WAb, with two upregulated genes and one downregulated gene, respectively (Fig. 4C). Furthermore, 38 and 135 genes were upregulated and downregulated, respectively, in the UAd vs UAb pool (Fig. 4D).

Fig. 4

Overview of DEGs numbers in four comparisons of mei. Volcano plot of DEGs in WAd and UAd (A), WAb and UAb (B), WAd and WAb (C), UAd and UAb (D) comparisons, VEEN analysis of DEGs in four comparisons (E). Red diamond represents up‐regulated DEGs, while blue square represents downregulated DEGs.

GO and KEGG pathway function enrichment analyses were then performed on the four compared groups. As shown in Fig. 5A, the top five significantly enriched GO terms (Padj < 0.05) in DEGs between WAd and UAd were all related to photosynthesis: photosynthesis, light harvesting, photosystem II, chlorophyll binding, and photosystem I; further details are listed in Table S3. Many DEGs were enriched in the oxidation–reduction reaction and stress response processes; some biosynthesis and metabolism processes were also enriched including the trehalose metabolic process, flavonoid biosynthetic process, and disaccharide biosynthetic process. Interestingly, the most significantly enriched terms for DEGs between WAb and UAb were also related to photosynthesis: photosynthesis, light harvesting, photosystem I and photosystem II (Fig. 5B, Table S3). Among DEGs between UAd and UAb, the most significantly enriched term was DNA‐templated transcription, followed by nucleic acid‐templated transcription (Fig. 5C). The remaining enriched GO terms were involved in the biosynthesis of various substances such as cellular nitrogen compounds, vitamin B6 and glutamine (Fig. 5C). Gene function annotation was performed for the three DEGs between WAd and WAb; based on homologous alignment, Pm026064 was predicted to be CALCIUM‐DEPENDENT PROTEIN KINASE 1, Pm019379 was predicted to be a dihydrodipicolinate reductase‐like protein related to photosynthesis, and Pm000744 was predicted to be DNA‐DIRECTED RNA POLYMERASE III SUBUNIT 1.

Fig. 5

GO and KEGG pathway enrichment analysis of DEGs in WAd and UAd (A/D), WAb and UAb (B/E), UAd and UAb (C/F) comparisons. Y‐axis represent GO/KEGG terms while X‐axis represent the rich factor of each term. The size of each point represents the number of DEGs in the term/pathway, and the color of the point represents the FDR (corrected P‐value).

KEGG enrichment analyses of DEGs between groups showed that photosynthesis‐antenna protein was the most significantly enriched DEG pathway in the WAd vs UAd comparison, followed by flavonoid biosynthesis, plant–pathogen interaction, and brassinosteroid biosynthesis (Fig. 5D). Many DEGs were enriched in the plant hormone signal transduction pathway. In the WAb vs UAb comparison, photosynthesis‐antenna protein was the most significantly enriched pathway, followed by photosynthesis, and protein processing in the endoplasmic reticulum (Fig. 5E). The most enriched pathways in UAd vs UAb were involved in metabolism‐related processes (Fig. 5F) including glyoxylate and dicarboxylate metabolism; carbon fixation in photosynthetic organisms; and alanine, aspartate, and glutamate metabolism. The three DEGs between WAd and WAb were assigned to four pathways: purine metabolism, RNA polymerase, pyrimidine metabolism and plant–pathogen interaction. Detailed KEGG pathway enrichment results are provided in Table S4.

Functional analysis of unique DEGs between different branch architectures

Many DEGs were differentially expressed only in a single comparison (termed unique DEGs) based on VEEN analysis of DEGs among the four compared groups (Fig. 4E). In contrast, WAd vs UAd contained the maximum number (1272) of unique DEGs; WAb vs UAb contained 289 unique DEGs (Fig. 4E).

To further elucidate the functions of DEGs on the same side of the branch, GO enrichment analysis was performed for unique DEGs in WAd vs UAd and WAb vs UAb. In WAd vs UAd, the top enriched GO terms were related to flavonoid and hormones biosynthesis and metabolism. The flavonoid biosynthetic process, flavonoid metabolic process and hormone metabolic process were enriched (Fig. 6A, Table S5). The most enriched GO terms in WAb vs UAb were related to photosynthesis, including photosystem, thylakoid part, and photosynthetic electron transport chain. KEGG pathway enrichment analyses revealed that flavonoid biosynthesis was the most significantly (Padj < 0.05) enriched pathway of unique DEGs between WAd and UAd, followed by brassinosteroid biosynthesis and cutin, suberin, and wax biosynthesis (Fig. 6C). Many unique DEGs between WAd and UAd were assigned to plant hormone signal transduction and plant–pathogen interaction. The remaining most enriched pathways were mainly related to metabolite synthesis, including monoterpenoid biosynthesis, galactose metabolism and anthocyanin biosynthesis. Among the top enriched pathways of unique DEGs between WAb and UAb, photosynthesis was the most significantly (Padj < 0.05) enriched pathway, followed by pathways related to metabolite synthesis, including diterpenoid biosynthesis, fatty acid biosynthesis, and porphyrin and phenylalanine metabolism (Fig. 6D, Table S6). Unique DEGs between WAd and UAd were functionally enriched in plant hormones and flavonoid metabolism‐related terms, whereas those between WAb and UAb were mainly related to photosynthesis.

Fig. 6

GO and KEGG pathway enrichment analyses of unique DEGs in WAd and UAd comparison (A/C), WAb and UAb comparison (B/D). Y‐axis represent GO/KEGG terms, while X‐axis represent the rich factor of each term. The size of each point represents the number of DEGs in the term/pathway, and the color of the point represents the FDR (corrected P‐value).

Screening for DEGs related to IAA and GA

Previous studies demonstrated that IAA and GA play important roles in branch development and formation in rice and Japanese cherry (Sugano et al. 2004, Wang et al. 2018). Plant hormone analyses of the four mei F1 population pools also indicated that IAA and GA highly significantly differed between the adaxial and abaxial sides of branches in the weeping and upright progeny pools (P < 0.01). To elucidate hormone functions in formation of the weeping habit in mei, DEGs related to IAA and GA were selected from the WAd vs UAd and WAb vs UAb comparisons for further analyses (Table S7). In WAd vs UAd, 22 IAA‐related DEGs and nine GA‐related genes were detected, whereas only six IAA‐related DEGs and four GA‐related DEGs were found in WAb vs UAb (Table S7). Notably, 61% of the hormone‐related DEGs were upregulated in WAd vs UAd (Table S7), which was consistent with the accumulation pattern of hormone content on the adaxial side of branches in the weeping and upright progeny pools (Fig. 3).

To further investigate the functions of DEGs related to IAA and GA in WAd vs UAd and WAb vs UAb, phylogenetic analyses of these DEGs were performed, together with functionally verified key genes that regulate the biosynthesis and transport of IAA (Jiang et al. 2017) and GA (Gao et al. 2017) in Arabidopsis. As shown in Fig. S3A, Pm005112, Pm013243 and Pm024306 appear to be involved in the regulation of IAA biosynthesis; they clustered with YUCCA6 (YUC6, At5g25620), CYTOCHROME P450 FAMILY 79 SUBFAMILY B POLYPEPTIDE 2/3 (CYP79B2/3, At4g39950/At2g22330), and TRYPTOPHAN AMINOTRANSFERASE OF ARABIDOPSIS (TAA, At1g23320), which has been identified as a key gene in the regulation of IAA biosynthesis. Pm007046 and Pm020838 may be related to polar transport of IAA because they clustered with LAZY and AUX1 (At2g38120), which exhibits polar transport of IAA in Arabidopsis (Jiang et al. 2017). Pm011672 and Pm004966 may participate in GA biosynthesis because they clustered with ENT‐COPALYL DIPHOSPHATE SYNTHASE/ENT‐KAURENE SYNTHASE (CPS/KS, At4g02780/ At1g79460) and GA3OX (At1g15550), which have been reported to regulate GA biosynthesis in Arabidopsis (Gao et al. 2017; Fig. S3B). Pm019920, which clustered with GIBBERELLIN INSENSITIVE DWARF (GID, At3g63010), may act as a GA receptor protein (Wang et al. 2017).

Screening for DEGs involved in secondary branch growth

Secondary branch architecture in mei shows substantially different characteristics from primary architecture (Zhang 2016, Hou et al. 2017), as well as secondary wall thickness and degree of lignification (Zhang 2016, Hou et al. 2017). Based on our functional annotation results, 86 DEGs involved in secondary branch growth were detected between WAd and UAd; among these, 64 and 22 DEGs were upregulated and downregulated, respectively (Table S7). Among the 34 DEGs between WAb and UAb that were related to secondary branch growth, 22 and 12 were upregulated and downregulated, respectively (Table S7).

Phylogenetic analyses were performed among these DEGs related to secondary branch growth in the WAd vs UAd and WAb vs UAb comparisons, together with verified functional genes reported to regulate secondary branch growth in Arabidopsis. As shown in Fig. S3C, Pm030127, Pm002521 and Pm019600 may be involved in lignin biosynthesis because they clustered with PHENYLALANINE AMMONIA LYASE1/2 (PAL1/2, At2g37040/At3g53260), Arabidopsis PEROXIDASE52 (AtPrx52, At5g05340) and 4‐COUMARATE:CoA LIGASE1/2 (4CL1/2, At1g51680/AT3G21240), which are reportedly involved in lignin biosynthesis in Arabidopsis (Zhao and Dixon 2011). Pm023145 may be related to the biosynthesis of mannose, because it clustered with CELLULOSE SYNTHASE LIKE D2/3/5 (CSLD2/3/5, At5g16910/At3g03050/At1g02730), which is involved in mannose biosynthesis in Arabidopsis (Verhertbruggen et al. 2011).

Analysis of differential expression transcription factors (TFs)

To identify various TFs in the mei transcriptome, genes were blasted to the PlantTFDB database (http://planttfdb.cbi.pku.edu.cn/index.php). In total, 196 differentially expressed TFs that belonged to 23 gene families were detected in WAd vs UAd (Table S8). Among the various families of TFs, MYB, WRKY, bHLH and NAC contained the largest numbers of differentially expressed TFs (Fig. S2A). Notably, 98 differentially expressed TFs were found, which belonged to 24 gene families in WAb vs UAb (Table S8). Among the various families of TFs, MYB, bHLH and C2C2 contained the largest number of differentially expressed TFs (Fig. S2B).

TF factors are critical molecular switches that regulate mRNA expression in a vast majority of associated genes (Joshi et al. 2016). Hence, TFs potentially related to transcriptional regulation of plant hormones and secondary growth in mei were identified by phylogenetic analyses of functionally verified genes in Arabidopsis. Pm030202 may be involved in IAA biosynthesis because they clustered with TCP13 (At3g02150; Fig. S3F), which regulate IAA biosynthesis in Arabidopsis (Shin et al. 2007, Huh et al. 2012, Liu et al. 2012, Zhang et al. 2018). Pm010085 may participate in GA biosynthesis because it clustered with ERF11 (At1g28370; Fig. S3D), which regulates GA biosynthesis in Arabidopsis (Zhou et al. 2016). Pm013020 and Pm009364, which clustered with MYB85 (At4g22680) and MYB58/63 (At1g16490/At1g79180; Fig. S3G), may be involved in lignin biosynthesis (Zhong et al. 2010).

Validation of candidate genes by quantitative reverse‐transcription polymerase chain reaction (qRT‐PCR)

To verify the expression levels of candidate genes, qRT‐PCR was performed using samples that were consistent with our BSR‐Seq analysis. For the nine randomly selected candidate genes, the same expression patterns were detected in the qRT‐PCR and BSR‐Seq results. The Pearson correlation coefficients of Pm013243, Pm030202, Pm010085, Pm007046, Pm011672, Pm030127, Pm020838, Pm019920 and Pm019600 were 0.78, 0.96, 0.88, 0.88, 0.88, 0.87, 0.84, 0.77 and 0.89, respectively; the overall Pearson correlation coefficient between the qRT‐PCR and BSR‐Seq data was 0.7695, indicating that the transcriptome data were reliable (Fig. S4).

In detail, Pm013243 and Pm030202 were significantly upregulated on both branch sides in weeping progeny, compared with upright progeny (P < 0.05). In contrast, Pm030127 and Pm019600 were both significantly (P < 0.05) downregulated on both branch sides in weeping progeny, compared with upright progeny. There was no observable pattern in the differences in gene exp. ression among Pm010085, Pm020838, and Pm019920 between weeping and upright progeny (Fig. 7). Within the same branch architecture, Pm011672 was significantly (P < 0.01) upregulated on the adaxial side compared with the abaxial side of branches in weeping progeny, and the opposite trend was observed among upright progeny; while the expression trend of Pm020838 was just contrary to Pm011672. Pm010085, Pm019920 and Pm019600 were significantly (P < 0.01) downregulated on the adaxial side, compared with the abaxial side, in upright progeny; an opposite trend was observed in Pm007046. There were no significant differences in expression between Pm013243, Pm030202 and Pm030127 on the adaxial and abaxial sides in weeping or upright progeny.

Fig. 7

Validation of candidate genes by qPT‐PCR experiment. Stars denote significant difference (*P < 0.05, **P < 0.01, t‐test, n = 3).

Discussion

DEGs were enriched in photosynthesis and flavonoid biosynthesis pathways between annual branches of weeping and upright progeny of the mei F1 population at the transcriptome level

Studies of the weeping habit in mei are limited to biological differences between weeping and upright branches; it has been speculated that the weeping branch habit is related to differences in lignin formation and hormone distribution among branches in mei and peach (Li 2006, Wang 2014). Physiological and molecular studies of the weeping habit mainly focus on these aspects of its generation; however, the key genes and regulatory networks responsible for weeping have not been clarified in mei (Lü and Chen 2003, Wang 2014, Hou et al. 2017). Using mei F1 population materials and a relatively consistent genetic background, we applied BSR‐Seq to further illuminate the biological processes that control traits related to weeping branches at the transcriptome level. As shown in Fig. 4, the numbers of DEGs in WAd vs UAd and WAb vs UAb were significantly greater than those in WAd vs WAb and UAd vs UAb. Analyses of unique DEGs also indicated greater numbers of DEGs in WAd vs UAd and WAb vs UAb, compared with groups that exhibited the same branch architecture. Among progeny with the same branch architecture, greater differences were detected among the same sides of the branch than among different sides of the branch. There may be different biological or molecular processes involved in branch architecture on different sides of branches between weeping and upright progeny.

To further study functions of DEGs between weeping and upright progeny, GO and KEGG enrichment analyses were conducted in this study. As shown in Fig. 5, photosynthesis/light harvesting (GO:0009765) was the most common GO enrichment term in DEGs of WAd vs UAd and WAb vs UAb. Photosynthesis/antenna proteins (map00196) was the most common pathway identified by KEGG pathway enrichment analysis of DEGs of WAd vs UAd and WAb vs UAb. GO enrichment results were similar for unique DEGs of WAb vs UAb, of which photosystem (GO:0009521) was the most common term; many of the top 15 enriched terms were also involved in the photosynthesis process (Fig. 6). Photosynthesis (map00195) was also the most commonly enriched pathway among unique DEGs of WAb vs UAb. However, the top 15 enriched GO terms and KEGG pathways for unique DEGs of WAd vs UAd were quite different, but exhibited highly enriched flavonoid‐ and phytohormone‐related GO terms and KEGG pathways. The plant heights of weeping progeny were significantly smaller (P < 0.01) than those of upright progeny in the mei F1 population (Table S1). Notably, DEGs assigned to photosynthesis‐related terms were all downregulated in the WAd pool, compared with the UAd pool; they were also downregulated in the WAb pool, compared with the UAb pool. The significant difference in plant height may explain why expression levels of photosynthesis‐related genes on the adaxial and abaxial sides were lower on weeping progeny branches than on upright progeny branches. In WAd vs UAd, 66.1% (80/121) of DEGs related to a stimulus response term (GO:0050896) were upregulated in the WAd pool, compared with the UAd pool, suggesting that the weeping branch habit may be an adverse variation of the upright branch habit caused by dwarf plant architecture.

DEGs between WAd and UAd were significantly enriched in flavonoid biosynthesis (GO:0009813; Padj < 0.05); the top two highly enriched terms of unique DEGs in WAd vs UAd were also related to flavonoid biosynthesis and metabolic processes. Flavonoid biosynthesis (map00941) was also highly enriched in both WAd vs UAd and WAb vs UAb. Most DEGs related to flavonoid biosynthesis were upregulated in the UAd pool, compared with the WAd pool; higher expression levels of DEGs related to flavonoid biosynthesis were observed in the UAb pool, compared with the WAb pool. Flavonoid participates in many aspects of plant physiology, including IAA transport (Peer and Murphy 2007) and stress response (Winkel‐Shirley 2002). Studies have suggested that flavonoids may indirectly regulate plant architecture by modulating the levels of IAA in many species (Peer and Murphy 2007, Buer and Djordjevic 2009); however, the specific mechanism remains unclear. The relationship between weeping branch architecture and flavonoid biosynthesis in mei requires further investigation.

Weeping habit‐related phytohormone and possible candidate analyses in mei

Phytohormones play important roles in the development of plant branch architecture in rice (Wang et al. 2018); branch morphogenesis has been shown to involve coordination of various plant phytohormones in many species (Mcsteen and Leyser 2013). In this study, LC–MS/MS analyses were conducted to further detect differences in phytohormones between weeping and upright progeny of the mei F1 population. We determined the contents of 10 phytohormones; IAA, ICA, IP, tZ, cZ, DZ, and GA3 contents significantly (P < 0.05) differed between weeping and upright annual branches (Fig. 3). The non‐uniform distributions of IAA and GA may responsible for the weeping habit in willow and peach (Li 2006, Liu et al. 2017). In this study, ABA, IAA, and GA3 showed the highest concentrations of hormones, whereas ABA content did not differ significantly between weeping and upright branches. Branches were highly significantly longer in weeping progeny of the F1 population than in upright progeny (P < 0.01; Table S1). Because GA3 and IAA have been reported to act as strong accelerators of branch growth in many species (Tanimoto 2005), longer branches in weeping progeny may be the result of high IAA and GA contents in weeping branches. The rapid elongation of weeping branches in P. spachiana is reportedly a vital factor for the development of weeping branches, due to the unmatched timing between branch elongation and the formation of lignified secondary xylem (Sugano et al. 2004); however, the detailed effects of IAA and GA3 in weeping mei branches require further investigation. Cytokinin content was lower than that of IAA and GA3; however, the roles of cytokinins in the regulation of branch architecture remain unclear.

IAA and GA3 contents were significantly greater on both sides of weeping branches than on the corresponding sides of upright branches (P < 0.01). This result is consistent with previous findings that IAA contents were greater in the base and middle parts of the branch in weeping mei branches, compared with upright mei branches; moreover, GA3 contents were greater in weeping branches of Japanese cherry and peach (Sugano et al. 2004, Li 2006, Wang 2014). This study is the first to determine GA3 and IAA contents on the adaxial and abaxial sides of weeping and upright mei progeny. GA3 content was significantly greater on the adaxial side of the branch than on the abaxial side in both weeping and upright branches (P < 0.01); IAA content was significantly greater only on the abaxial side in upright branches (P < 0.01; Fig. 3). Non‐uniform distribution of GA may affect formation of the pendulous branch phenotype in mei (Wang 2014). However, external application of GA3 did not relieve the weeping habit in mei branches (Lü and Chen 2003). Although non‐uniform distribution of GA3 and IAA in mei branches has been reported, further experiments are needed to confirm this finding.

DEGs related to IAA and GA3 were screened firstly in each group comparison. Based on homologous and phylogenetic analyses, six genes (Pm005112, Pm013243, Pm007046, Pm024306, Pm020838 and Pm030202) and four genes (Pm011672, Pm019920, Pm004966 and Pm010085) that may play important roles in IAA and GA3 biosynthesis were identified separately (Fig. S3). Because the expression levels of Pm005112, Pm024306 and Pm004966 were extremely low (RPKM <20), we performed qPCR for the remaining four IAA‐related DEGs and three GA‐related DEGs. As shown in Fig. 7, the expression of Pm013243, which may be a key gene in IAA biosynthesis, and Pm030202, a TCP family member that promotes IAA biosynthesis, were significantly (P < 0.05) upregulated on both sides of annual shoots in weeping progeny, compared with upright progeny. The higher expression levels of Pm013243 and Pm030202 in weeping progeny may contribute to greater IAA content in weeping progeny. Pm007046 and Pm020838, two putative IAA polar transport genes, were both significantly (P < 0.05) upregulated in the UAd pool, compared with the UAb pool; however, they were downregulated in the WAd pool, compared with the WAb pool. Higher expression levels of Pm007046 and Pm020838 in the UAd pool may be associated with greater IAA content on the abaxial side of branches in upright progeny. As shown in Fig. 7, Pm011672, a potential key gene in GA biosynthesis, and Pm010085, an AP2/ERF family member that promotes GA3 biosynthesis, were highly significantly (P < 0.01) upregulated on adaxial sides of annual branches in weeping progeny, compared with upright progeny. Pm019920, a putative GA receptor, was significantly (P < 0.05) upregulated in the WAd pool, compared with the UAd pool. The higher expression levels of Pm011672, Pm010085 and Pm019920 in the adaxial side of weeping progeny may be related to greater GA3 content in weeping progeny, compared with upright progeny. Only one key gene in the cytokinin synthesis pathway, LONELY GUY (LOG, Pm005529), was differentially expressed; it exhibited a higher expression level in weeping progeny.

Weeping habit may be related to differential secondary growth during branch elongation

Field observation of the ‘LiuBan' × ‘FenTaiChuiZhi’ mei F1 population showed no clear weeping habit in the early stage of annual branch growth. However, the weeping habit becomes obvious during the lignification process in annual branches of mei (Zhang 2016), as well as in weeping cherry (Sugano et al. 2004). In peach, elongation growth is faster than secondary growth in weeping branches; this results in the weeping habit, because the branch cannot support its own weight (Li 2006). Previous studies have suggested that the secondary cell wall of the extension zone is thinner in weeping mei branches than in upright branches; moreover, lignification is delayed in weeping branches (Zhang 2016, Hou et al. 2017). Lignification is greater on the abaxial side than on the adaxial side in upright mei branches, whereas the two sides do not differ significantly in weeping branches (Zhang 2016, Hou et al. 2017). Based on these findings, we speculated that differential secondary branch growth may occur between weeping and upright progeny in the mei F1 population. In many woody plants, secondary cell wall thickening and branch lignification during secondary branch growth provide strong mechanical support in the branch (Zhong et al. 2010, Wang et al. 2013b).

Homologous and phylogenetic analyses identified six genes (Pm030127, Pm002521, Pm019600, Pm023145, Pm013020 and Pm009364) that may play important roles in secondary branch growth in mei. Pm030127 and Pm019600 were verified by qPCR analysis; Pm002521, Pm023145, Pm013020 and Pm009364 showed low expression (RPKM <20). As shown in Fig. 7, Pm030127 and Pm019600, which are both key enzymes in lignin biosynthesis, were significantly (P < 0.05) upregulated on both sides of annual branches in upright progeny, compared with weeping progeny. Downregulation of Pm030127 and Pm019600 in weeping branches, compared with upright branches, may be related to delayed secondary growth of weeping branches compared with upright branches; further investigation is needed to test this hypothesis.

In our previous study, a locus conferring the weeping habit was detected on a 1.14‐Mb interval of pseudochromosome 7 in mei (Zhang et al. 2015). SNP analysis showed that 55% (68/125) of SNPs with ΔSNP index values ≥0.8 were also distributed on pseudochromosome 7. Annotation results suggested the 1.14‐Mb interval contains 159 predicted protein encoding genes (Zhang et al. 2015). Here, we determined the expression levels of these weeping candidate genes. Nine genes (Pm024143, Pm024158, Pm024169, Pm024201, Pm024226, Pm024234, Pm024247, Pm024271 and Pm024298) in the 1.14‐Mb interval were differentially expressed between WAd and UAd; among these, Pm024169, Pm024247, Pm024271 and Pm024298 were upregulated in WAd, compared with UAd. Five genes (Pm024123, Pm024148, Pm024234, Pm024247 and Pm024293) were differentially expressed between WAb and UAb; Pm024148 and Pm024247 were upregulated in WAb, compared with UAb. Only two genes (Pm024169 and Pm024240) were differentially expressed between UAd and UAb; both were upregulated in UAb, compared with UAd. No DEGs were detected between WAd and WAb. Notably, Pm024158, Pm024201 and Pm024293 were all related to GO:0016021, which was related to integral membrane components. Recently, increasing effort has been expended in identification of weeping habit candidate genes. WEEP, an extremely conserved protein among many plants—for which the function is unknown—was determined to cause the weeping habit in weeping peach shoots (Hollender et al. 2018). Homologous analysis suggested that Pm015033 in mei was highly conserved with WEEP. However, its differential expression was not detected in any of the compared groups. TILLER ANGLE CONTROL 1 (TAC1) has been reported to promote outward lateral shoot growth and orientation in poplar, peach and Arabidopsis (Xu et al. 2017, Hollender et al. 2018). Weeping mei trees have a larger shoot angle than upright mei trees (Li et al. 2019), and TAC1 may act upstream of WEEP (Hill and Hollender 2019). Pm018319, which is conserved with TAC1, was differentially expressed between WAd and UAd, as well as WAb and UAb; FC ratios were 1.738 and 1.502, respectively.

In summary, the weeping habit is a complex process regulated by a series of metabolic pathways (Figs 4 and 5). Functional analyses showed that photosynthesis and flavonoid biosynthesis were highly enriched among DEGs of weeping and upright progeny. Differences in growth characteristics and phytohormone content were detected among progeny with different branch types in the F1 population. Functional annotation and homologous analysis showed that 10 genes related to IAA and GA3 biosynthesis, together with 6 genes related to secondary branch growth, merit further research; 13 DEGs found in the weeping candidate region of chromosome 7, as well as TAC1, also warrant detailed functional analysis.

Author contributions

The order of authorship is a joint decision of the co‐authors. J.Z. conceived and designed the experiments. T.‐Y.M., Y.‐Y.L., H.‐H.Z. analyzed the data. Q.‐X.Z., M.‐Z.B. and R.‐Q.J. contributed reagents/materials/analysis tools. T.‐Y.M and J.Z wrote the paper. J.‐W.Z. revised the manuscript. All authors read and approved the final manuscript.

Fig. S1. Schematic diagram of sampling of upright (A) and weeping (B). The arrows indicate the direction of gravity.

Fig. S2. Statistics of gene families of differential expressed transcription factors in WAd and UAd comparison (A) or WAb and UAb comparison (B).

Fig. S3. Phylogenetic analysis of candidate genes and its functional verified homologs in Arabidopsis.

Fig. S4. Pearson correlation analysis of FPKM obtained from RNA‐Seq and qRT‐PCR results of nine selected genes.

Click here for additional data file.

Table S1. The length of annual branches and plant height of selected weeping progenies and upright progenies.

Click here for additional data file.

Table S2. Summary of BSR‐Seq and mapping data using the Prunus mume genome as a reference.

Click here for additional data file.

Table S3. GO function enrichment of DEGs in WAd vs UAd, WAb vs UAb, WAd vs WAb and UAd vs UAb comparisons.

Click here for additional data file.

Table S4. KEGG pathway enrichment analyses of DEGs in different comparisons.

Click here for additional data file.

Table S5. GO function enrichment of unique DEGs in WAd vs UAd and WAb vs UAb comparisons.

Click here for additional data file.

Table S6. KEGG pathway enrichment analyses of unique DEGs in WAd vs UAd and WAb vs UAb comparisons.

Click here for additional data file.

Table S7. Detail information of DEGs involved in the biosynthesis of auxin and GA3 and secondary growth of branch in WAd vs UAd and WAb vs UAb comparisons.

Click here for additional data file.

Table S8. Detailed information of the differentially expressed TFs in WAd vs UAd, WAb vs UAb and UAd vs UAb.

Click here for additional data file.

Table S9. Primers used in the qRT‐PCR validation experiment.

Click here for additional data file.