Two M. indica cultivars “Irwin” (IR, deep yellow-fleshed) and “Baixiangya” (BY, pale yellow-fleshed), widely cultivated in China, were chosen from the orchard of the South Subtropical Crops Research Institute, Zhanjiang (110°4′E, 21°12′N), Guangdong Province, China, for this study.

Genomic DNA was extracted from the fresh young leaves of mango cultivar “Irwin” according to the modified CTAB method, and used for both Illumina and PacBio sequencing. DNA quantity and quality were assessed using a NanoDrop 2000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and electrophoresis on a 0.8% agarose gel. For Illumina sequencing, we constructed a short-read (300 bp) library using the NEBNext Ultra DNA Library Prep Kit (New England BioLabs, Ipswich, MA, USA). For PacBio sequencing, 20-kb SMRTbell libraries were prepared according to the protocol released by PacBio, and then the library was sequenced on the PacBio Sequel System (Biomarker Technologies Corp., Beijing, China) with one SMRT cell, yielding more than 20 Gb of raw data.

A total of 18.87 Gb high-quality-filtered data were used to perform a 17-kmer analysis. The genome size, heterozygous rate, and repeat rate were estimated using the program Jellyfish (Marçais and Kingsford, 2011). MECAT (version 1.2) (Xiao et al., 2017) was used for the correction and assembly of PacBio data, and Polish (version 1.22) (Walker et al., 2014) was used for the assembly correction with Illumina short reads. To check the completeness of the assembly, the RNA-seq data from fruits were mapped to the genome assembly using BLASTN (version 2.6.0+) (Altschul et al., 1990); various sequence homology and coverage parameters were determined. Furthermore, we evaluated genome completeness using the Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis (Simão et al., 2015). The Hi-C sequencing of the chromosome-level genome assembly was performed according to the previously described method (Belton et al., 2012).

The present genome was annotated with three characteristics. (i) We first used de novo prediction and homology-based alignment to predict repetitive sequences across the mango genome. The repeat sequence database was constructed using the LTR_FINDER (Xu and Wang, 2007), PILER-DF (Edgar and Myers, 2005), and RepeatScout (Price et al., 2005). RepeatMasker and RepeatProteinMask software (Tamura, 1992) were used to identify repeat sequences. (ii) We then annotated the predicted protein-coding genes included in the databases, namely, National Center for Biotechnology Information non-redundant database, Gene Ontology, Kyoto Encyclopedia of Genes and Genomes, SwissProt, and EuKaryotic Orthologous Groups. (iii) Finally, we predicted miRNA and snRNA sequences using the INFERNAL software (Nawrocki and Eddy, 2013). The rRNA was identified using the BLAST program, and tRNAscan-SE was used to predict tRNA (Chan and Lower, 2019).

Healthy fruits of two cultivars (“Irwin” and “Baixiangya”) in three distinct developmental stages were sampled. The developmental stages were: S1, early development stage, 40 days after full boom (DAFB); S2, harvest stage, 86 DAFB; and S3, fully ripe stage, 96 DAFB. The flesh of the collected fruit was separated from the peel and frozen in liquid nitrogen for β-carotene measurement and RNA sequencing. All experiments were performed in triplicate.

Total RNA was extracted from the frozen flesh tissue of each sampling point and each cultivar using a Quick RNA Isolation Kit (Invitrogen, Carlsbad, CA, USA). The cDNA libraries were assembled according to the protocol of the manufacturer of the NEBNext Ultra RNA Library Prep Kit for Illumina (E7530) and NEBNext Multiplex Oligos for Illumina (E7500), and the libraries were sequenced on the Illumina HiSeq XTen platform (Biomarker Technologies Corp., Beijing, China). The DESeq2 software (Robinson et al., 2010) was employed to identify differentially expressed genes, which were filtered according to p ≤ 0.05, and expression value >0.

The β-carotene in the flesh was extracted as described by Ma et al. (2018). The carotenoids were analyzed using the Waters Ultra-high-performance liquid chromatography (Milford, MA, USA) with a photodiode array detector and a single quadrupole mass spectrometer detector in series (6120 Quadrupole, Agilent, Santa Clara, CA, USA). The column used was a YMC carotenoid (C30) column (Wilmington, NC, USA). The eluent phases were as follows: mobile phase A, acetonitrile:methanol = 3:1 (v/v); and phase B, 100% methyl tert-butyl ether (MTBE). Each eluent contained 0.01% butylated hydroxytoluene (BHT). Gradient elution was performed as follows: 0–2 min, 85:15 A:B; 2–4 min, 75:25 A:B; 4–7 min, 40:60 A:B; 7–10 min, 40:60 A:B; 10–13 min, 5:95 A:B; 13–23 min, 85:15 A:B. The flow rate was 0.8 ml/min, and the injection volume was 5 μl.

Total RNA was extracted from frozen flesh using TRIzol reagent (TaKaRa, Dalian, China) according to the instructions of the manufacturer. The primers used for real-time PCR were designed using Primer 5 and are listed in Supplementary Table 10. The PCR mixture (20 μl total volume) comprised 10 μl of LightCycler 480 SYBR Green I Master Mix (Roche, Mannheim, Germany), 2 μl of each primer, 2 μl of diluted cDNA, and 6 μl DEPC-H₂O. Reverse transcription (RT)-PCR) conditions were as follows: pre-incubation at 95°C for 5 min, followed by 45 cycles at 95°C for 5 s, 58°C for 15 s, and 72°C for 10 s. The expression levels of genes of interest were calculated according to the formula 2^−ΔΔCt.

Yeast one-hybrid (Y1H) assay was performed according to the Gold Yeast One-Hybrid Library Screening System User Manual (Clontech, Palo Alto, CA, USA). Either the MiPSY1 or MiPSY1-mut promoter was inserted into the Phis2 vector to generate the baits. MiPSY1-mut is a mutated version of MiPSY1 in which the core sequence (CACGTG) of two E-box elements has been altered and is “TCTAGC.” The full-length coding sequences (CDS) of MibZIP66 and MibHLH45 were cloned into the pGADT7 vector to generate the prey and were transformed into the yeast strain Y187 containing the bait. The TF–promoter interaction was determined according to the growth ability of the co-transformed yeast cells on synthetic defined -Trp-Leu-His media with 3-amino-1,2,4-triazole.

The full-length CDS of MibZIP66 and MibHLH45 was cloned into the vector pGreenII0029 62-SK as effectors. Meanwhile, the MiPSY1 promoter was recombined with pGreenII0800-LUC to create the reporter vector. All vectors were transformed into Agrobacterium tumefaciens (GV3101), and the cultures were adjusted to an optical density at OD 600 nm of 0.1 with an infiltration medium containing 100 μM acetosyringone, 0.5 M MES, and 10 mM MgCl₂. Agrobacterium with TFs or promoters was injected into Nicotiana benthamiana leaves using a 5-ml injection syringe. Transgenic plants were grown in a glasshouse condition for 3 days. Firefly luciferase (LUC) and Renilla luciferase (REN) were determined using the Dual-Luciferase Reporter Assay System (Promega Corp., Madison, WI, USA). Transcriptional activation was calculated as the ratio of LUC to REN.

The full-length coding sequences of MibZIP66 and MibHLH45 without the stop codon were amplified using PCR and cloned into the GFP vector, under the control of the CaMV 35S promoter. The vector constructs of 35S-MibZIP66-GFP and 35S-MibHLH45-GFP were transformed into Agrobacterium tumefaciens strain GV3101 through electroporation and then injected into tobacco leaves. MiPSY1 or the control vector (pYBA1132) plasmid was co-transformed with a plasmid coding for a mitochondrial marker red fluorescent protein (RFP) into Arabidopsis protoplasts by polyethylene glycol (PEG) transformation method according to Lister et al. (2007). Localizations of MiPSY1, MibZIP66, and MibHLH45 were examined using a fluorescence microscope (Zeiss Axioskop 2 Plus; Oberkochen, Germany).

The genome of the mango cultivar “Irwin” was initially estimated to be 327 Mb with a 1.37% heterozygosity rate and 47.28% of the repeat sequence ratio based on the 17-mer sequence (Supplementary Figure 1 and Supplementary Table 2). In total, 20.26 Gb of sequencing data were produced with an average coverage depth of 61.96× using the PacBio single molecular long-read method (Supplementary Table 1). After heterozygosis was accounted for, a draft genome of 396 Mb was assembled, with N50 of 1.03 Mb (Table 1). The assembly of our genome was comparable to that of the recently published genome of the traditional Indian cultivar “Alphonso” (392.9 Mb) (Wang et al., 2020), but better than that of the Southeast Asian variety “Hong Xiang Ya” (~372 Mb) (Li et al., 2020). Then, Hi-C contact frequency derived from Hi-C sequencing was used to order and orient the polished contigs into Hi-C scaffolds. As a result, a total of 97.19% (375 Mb) sequences in the mango genome were anchored to 20 chromosomes, and the chromosome size ranged from 12.8 to 29.8 Mb (Supplementary Tables 3, 4, and Figure 1A). These results are consistent with those of studies by Li et al. (2020) and Wang et al. (2020). Two independent methods were used to evaluate assembly quality. The alignment rates of the RNA-seq reads from the three fruit developmental stages of the two cultivars were ~90.29 and 95.87% (Supplementary Table 5). In contrast, the BUSCO assay showed that ~94% of the core eukaryotic genes (1,991) were retrieved in the assembly, further confirming the continuity and high quality of the assembled genome (Supplementary Table 6). Furthermore, we performed the whole-genome alignment of our assembly and previous mango assembly. One-to-one syntenic blocks showed ~79.61% of our genome sequence matched 84.75% of the “Hong Xiang Ya” genome sequence (Figure 1B), indicating our genome sequence was substantially more complete.

We identified a total of 138.29 Mb of repeated sequences, representing 34.92% of the assembled mango genome. Among these repetitive elements, 12.77% were unknown repetitive sequences. The most abundant repeat type was long terminal repeats (LTRs), accounting for 15.49% of the genome. Additionally, we also found long interspersed nuclear elements (LINEs, 0.57%) and DNA transposons (6.25%) (Supplementary Table 7). Combined de novo prediction, protein-based homology searches, and RNA-seq analysis identified a total of 36,756 protein-coding genes with an average coding sequence length of 1,142 bp and 5.5 exons per gene (Supplementary Table 8). Moreover, a total of 109 rRNAs, 75 snRNAs, 483 tRNAs, 279 snoRNA, and 131 miRNAs were identified.

Studying the changes in β-carotene and the regulation of its biosynthesis will provide insight into the genetic basis of mango fruit flesh color. During fruit development and ripening, β-carotene content increased in “Irwin” mangoes, but remained low in “Baixiangya” mangoes. At the fully ripe stage, β-carotene concentrations in “Irwin” flesh was up to 50.06 μg/g of fresh weight, but was only 2.64 μg of fresh weight in “Baixiangya” flesh (Figure 2B). This difference in β-carotene concentration was also evident visually by the distinct color of the fruits (Figure 2A).

Comparative RNA-seq analysis was performed between “Irwin” and “Baixiangya” mangoes to study the expression profiles of β-carotene biosynthesis-related genes. A total of 52 genes were found to be associated with the carotenoid biosynthesis pathway, and 39 genes were differentially expressed between three developmental stages in both cultivars (Figure 3). Among these differentially expressed genes, the expression levels of MiDXS (Mango_gene17478 and Mango_gene23777, encoding 1-deoxy-d-xylulose-5-phosphatesynthase) and MiPSY (Mango_gene16452 and Mango_gene22766) were relatively low during the early fruit development stage but increased greatly with continued fruit ripening in both cultivars. LCYB (Mango_gene10477) showed a considerably higher expression level in “Irwin” than in “Baixiangya” mangoes at the early stage of fruit development. During the fruit ripening period, its expression level remained nearly unchanged in “Irwin,” but decreased in “Baixiangya” mangoes. Similarly, Mango_gene22766 on chromosome 7 showed a markedly high level of expression in “Irwin” during fruit development and ripening but a very low level in “Baixiangya” mangoes. At the fully ripened stage (S3), the Mango_gene22766 expression level was 14-fold higher in “Irwin” than in “Baixiangya” based on qRT-PCR data (Figure 4). Hence, the high expression level of these genes is likely to lead to the accumulation of β-carotene in the Mango fruit. In contrast, the transcript levels of MiCCD (Mango_gene31763, Mango_gene29128, and Mango_indica_newGene_586) exhibited a significant decline during the fruit maturation phase, indicating a likely negative role in the development of Mango fruit flesh color. According to the annotation of the genome, Mango_gene22766 plays an important role in driving the metabolic flux toward the β-carotene formation. Therefore, it can be inferred that Mango_gene22766 may affect color formation by regulating β-carotene production in the fruit flesh. We sequenced two coding sequences of Mango_gene22766 from the flesh of “Irwin” and “Baixiangya,” and found they were both 1,299 bp long and encoded a peptide of 432 amino acids long (Supplementary Tables 11, 12). Except for the 400th and 417th amino acid, the rest of deduced amino acids were identical between the two cultivars (Supplementary Figure 2).

We named the cloned full-length sequence of Mango_gene22766 as MiPSY1 (Supplementary Table 11). MiPSY1 was suggested to be distributed in the mitochondria by the WoLF PSORT Prediction online software (http://www.genscript.com/wolf-psort.html). To examine its subcellular location, MiPSY1 was fused to RFP and transiently expressed in Arabidopsis protoplasts. Fluorescence microscopy analysis showed that the MiPSY1 protein was localized in the mitochondria (Figure 5). In summary, MiPSY1 might be a critical candidate gene controlling β-carotene formation in the mango fruit flesh.

To further investigate the potential TFs regulating MiPSY1 expression, a weighted gene co-expression network analysis of differentially expressed genes was performed (Supplementary Figure 3A). A total of 14 modules were identified (Supplementary Figure 3B). Analysis of module–trait relationships showed that the “MElightcyan” module was significantly correlated with β-carotene content (r = 0.86, p = 0.00) (Supplementary Figure 3B). Moreover, 49 TFs (including bHLH, bZIP, and NAC) were gathered in this module and were predicted to potentially interact with MiPSY1 (Supplementary Table 9). According to PLACE and Plant-CARE databases, seven TFs, namely MibZIP36 (Mango_gene12953), MibZIP66 (Mango_gene35447), MibHLH105 (Mango_gene17414), MibHLH31 (Mango_gene17464), MibHLH64 (Mango_gene20675), MibHLH45 (Mango_gene30145), and MibHLH130 (Mango_gene30794), were predicted to be regulatory cis-acting elements of the MiPSY1 promoter (Supplementary Table 13). Thus, we cloned these seven TFs and analyzed their sequences. Further phylogenetic analysis indicated that the seven identified TFs were grouped with Arabidopsis thaliana bZIP or bHLH proteins. To confirm the RNA-seq results, the seven TFs identified were further analyzed using qRT-PCR, which validated the RNA-seq results. The transcript levels of the seven TFs were significantly correlated with MiPSY1 expression and β-carotene content during fruit development and ripening (Figure 4). In summary, these results suggest that these seven TFs might share the regulatory function of A. thaliana bZIP or bHLH and act as candidate regulators of β-carotene biosynthesis.

Dual-luciferase transient expression assays were performed to test the potential regulatory effects of the key candidate TFs on the MiPSY1 promoter. The LUC/REN assays indicated that MibZIP66 and MibHLH45 could significantly activate the MiPSY1 promoter, with transcription activation increments of 2.7- and 2.3-fold, respectively. The effect of the other five TFs on the MiPSY1 promoter was limited. The action of MibZIP66 and MibHLH45 in combination had no significant additive effect on the MiPSY1 promoter activity (Figure 6). MibZIP66 and MibHLH45 were localized exclusively in the nucleus (Figure 7).

The PlantTFDB analysis suggested that MibZIP66 and MibHLH45 could bind to the CACGTG sequence in the MiPSY1 promoter, which was consistent with previous studies reporting that many bZIP and bHLH TFs specifically bind to the ACGT-containing elements, such as G-boxes (Pires and Dolan, 2010; Wolfgang et al., 2018). In our study, two CACGTG sequences were found to be located between the positions 852 and 800 bp, upstream of the transcriptional start codon in MiPSY1 (Supplementary Table 13). YIH experiments and cis-element mutagenesis were performed. As shown in Figure 8, all of the yeast cells could grow normally on two dropout minimal media (SD/-Leu/-Trp). Nevertheless, only the yeast cells co-transformed with the positive control, pGADT7-MibZIP66+MibPSY1-WT promoter, or pGADT7-MibHLH45+MibPSY1-WT promoter grew well on selective medium SD/-His/Leu/-Trp supplemented with 60 or 90 mM 3-amino-1,2,4-triazole. The pGADT7-MibZIP66+MibPSY1-Mut promoter and pGADT7-MibHLH45+MibPSY1-Mut promoter could not grow on the selective medium when the two G-box elements CACGTG of the MiPSY1 promoter were both mutated to TCTAGC (Figure 8). These results suggest that MibZIP66 and MibHLH45 could specifically bind to the G-box motif (CACGTG) of the MiPSY1 promoter to regulate the expression of MiPSY1 in the mango flesh.