Genomic Resequencing Unravels the Genetic Basis of Domestication, Expansion, and Trait Improvement in Morus Atropurpurea

Abstract

Mulberry is an economically important plant in the sericulture industry and traditional medicine. However, the genetic and evolutionary history of mulberry remains largely unknown. Here, this work presents the chromosome-level genome assembly of Morus atropurpurea (M. atropurpurea), originating from south China. Population genomic analysis using 425 mulberry accessions reveal that cultivated mulberry is classified into two species, M. atropurpurea and M. alba, which may have originated from two different mulberry progenitors and have independent and parallel domestication in north and south China, respectively. Extensive gene flow is revealed between different mulberry populations, contributing to genetic diversity in modern hybrid cultivars. This work also identifies the genetic architecture of the flowering time and leaf size. In addition, the genomic structure and evolution of sex-determining regions are identified. This study significantly advances the understanding of the genetic basis and domestication history of mulberry in the north and south, and provides valuable molecular markers of desirable traits for mulberry breeding.

Keywords: domestication; evolutionary history; flowering time; genome; mulberry.

Figure 1

De novo assembly of the M. atropurpurea ‘Huiqiu1’genome and gene family expansions and contractions throughout evolution. a) The plant of M. atropurpurea was used for genome assembly in this study. i) tree; ii) dormant bud; iii) axillary bud; iv) male flower; v) female flower. b) Circos plot showing the landscape of the M. atropurpurea genome. Outermost to innermost tracks indicate the a) pseudochromosomes, b) GC content, c) gene density, d) repeat density, e) gene expression, f) pseudogene density , and g) single nucleotide polymorphism (SNP) density. c) Species tree with molecular dating and gene family expansion and contraction indicated. The phylogenetic tree was constructed, and the divergence time (in MYA) was estimated based on single‐copy genes from 12 gene families. Gene family expansions and contractions on different lineages and species tips are indicated with green and red colors, respectively. The red dot represents the correction of time point using the speciation time of A. thaliana and P. trichocarpa (102.0–113.8 MYA), M. truncatula and C. lanatus (89.2–104.5 MYA), M. domestica and C. sativa (73.6–90.2 MYA) based on TIMATREE5 (http://timetree.org/).

Figure 2

Phylogenetic relationships and population structure of 425 resequenced mulberry accessions. a) Phylogenetic tree of mulberry accessions. Cannabis sativa was used as an outgroup. Color codes of accessions are consistent throughout Figure a–d (orange, M. alba (MA); purple, M. multicaulis (MM); brilliant blue, interspecific hybrid from Japan (JIH); pink, interspecific hybrid from China (CIH); red, modern elite cultivars M. atropurpurea (MECMA); green, Landrace2; dusty blue, Landrace1; beige color, wild). b) LD decay‐distance analysis. c) Principal component analysis (PCA) plots of the mulberry populations. d) Population structure based on different numbers of clusters (K = 2–5). The x axis indicates the mulberry groups with all accessions arranged in the same order as in a. The left y axis quantifies genetic diversity in each accession, which is represented by a vertical color‐coded column.

Figure 3

Demographic history of mulberry. a) Geographical distribution of mulberry accessions used in this study. Arrows indicate potential dispersal routes of mulberry. The red star indicates southwest China as a potential origin region of domesticated M. atropurpurea. b) Nucleotide diversity (Pi) among six different mulberry groups. The numbers within the circles represent the mean value of nucleotide diversity in each group. ANOVA analysis using LSD and Dunnett's T3 tests showed that significant difference (p < 0.05) was detected in each pairwise comparison except Landrace1 versus Landrace2 (p = 0.620) in Pi. c) Phylogenetic tree and genetic drift between different species using the f3‐statistic method. The number on the branch indicates branch length and the percentage indicates the proportion of admixture. d) Historical effective population size of mulberry. The shaded pink column represents the time interval of the LGM and PGP. LGM, Last Glacial Maximum (26.5–19 KYA). PGP, Penultimate Glacial Period (130–115 KYA). e–g) Admix proportion of accessions from interspecific hybrids in e) JIH (interspecific hybrid from Japan) group and f) CIH (interspecific hybrid from China). Different lowercase letters indicate a significant difference between each group (LSD and Dunnett's T3, p < 0.05). g) Average admix proportion of JIH and CIH from different mulberry subpopulations. MA (Morus alba), MM (M. multicaulis), and M. atropurpurea (including Landrace1, Landrace 2, and MECMA (elite cultivars) subgroups).

Figure 4

Genomic loci associated with leaf size during M. atropurpurea domestication. a,b) Box plots showing significant differences in leaf size (a) and leaf weight (b) between the MECMA and landrace2 and landrace1 groups. The MECMA group with larger leaves is widely cultivated in Guangdong. ** indicates p < 0.01 with Student's t‐test. c) Manhattan plot of GWAS and f) selective sweeps for leaf size and leaf weight. Significant overlapping loci were identified on chromosome 7 and highlighted by the red dashed column. The blue line indicates significance threshold −log₁₀(p) ≥ 8 in GWAS. d) Comparison of leaf development in four different stages between large leaf variety (“Tang10”) from the MECMA subgroup and small leaf variety (“Luozhi4”) from the Landrace2 subgroup. e) Bar plot showing significant differences between large‐ and small‐leaf varieties in the S4 development stage (12 d after bud sprouting). ** indicates p < 0.01 with Student's t‐test. g) Expression pattern of candidate genes across four different leaf development stages. L/S, Large leaf/Small leaf. h) Leaf size‐based association mapping located candidate loci between 9 747 000 and 9 793 000 on chromosome 7 and pairwise linkage disequilibrium (LD) analysis between the significant SNPs. Five genes encoding endo‐1,4‐beta‐xylanase (MaBXY) were tandemly duplicated and indicated by blue arrows (all gene annotations are shown in Table S10, Supporting Information). Pairwise LD block analyses (r² > 0.9) of SNPs with −log₁₀(p) ≥ 8. i) Leaf size feature of MaBXY transgenic Arabidopsis seedlings.

Figure 5

Identification of the MaERF110 gene associated with flowering time. a) Flowering time in different groups. Statistical significance was determined using a two‐sided t‐test. b) Manhattan plot for GWAS on flowering time. The blue line indicates the significance threshold with −log₁₀(p) ≥ 8. c) Selection signatures by calculating genetic differentiation (F _ST) between early and late flowering accessions across the genome in 100‐kb sliding windows with a step size of 20 kb. The orange dashed column indicates significant SNPs overlapping with GWAS hits, which located the MaERF110 gene. d) MaERF110‐based association mapping and pairwise LD analysis. Triangles show SNPs within the MaERF110 gene. Haplotype SNPs in the coding region have large effects and were highlighted in red. Strong LD with the lead SNP is connected to the pairwise LD with grey solid lines. e) Haplotypes (Hap) of MaERF110 among mulberry varieties. Statistical significance was determined using a two‐sided t‐test. The flowering time distribution of each haplotype group is displayed as a bar plot. f) Proportion of haplotypes of MaERF110 in different mulberry groups. g) Expression levels of MaERF110 in two early flowering time mulberry trees and two late flowering time mulberry trees. h) Expression levels of MaAP1 in flowering tissues from early and late flowering time mulberry trees. Data are presented as the mean ± S.E. (n = 3 independent RNA‐seq experiments). FPKM, fragments per kilobase per million reads. i) Flower time feature of MaERF110 transgenic Arabidopsis seedlings.

Figure 6

Sex determination region on chromosome 6 in M. atropurpurea. a) Manhattan plot of the results of GWAS for sex determination with the female genome “Tang 10” as a reference. The blue line indicates the significance of associated SNPs with a threshold value of −log₁₀(P) ≥ 8. b) Genetic differentiation (F _ST) between female and male trees across the genome based on the F _ST calculation of SNPs in 100‐kb sliding windows with a step size of 20 kb. The orange dashed column indicates the peaks the of sex divergence region overlapping with the peak of GWAS. c) Sex‐based association mapping located significant loci between 21 260 000 and 21 530 000 on chromosome 6. The squares represent SNPs (−log₁₀(P) ≥ 8) across the sex determination region, and nonsynonymous SNPs are highlighted in red. The candidate genes are indicated by purple (forward) and yellow (reverse) boxes. Functional annotation on these genes is shown in Table S17, Supporting Information. d) Gene expression levels of candidate genes located in SDR‐Y region. T1, stage of dormant bud sprouting; T2, stage of early axillary bud; T23, stage of late axillary bud. The star indicates the pentatricopeptide repeat gene (PPR). e) Whole‐genome alignments of each chromosome in the male reference against the female M. atropurpurea “Tang 10” reference. f) Whole‐genome alignments of chromosome A against the chromosome B in male M. atropurpurea “Huiqiu1.” g) Genome structure plot showing male‐to‐female differences in SDR and the distribution of male specifically genes across the Y‐SDR region. PPR indicates pentatricopeptide repeat gene. h) Statistics of the heterozygosity level of the 558 sex associated SNPs (p value < 1e‐8). Male individuals have higher heterozygosity than female.