Response of miR156-SPL Module during the Red Peel Coloration of Bagging-Treated Chinese Sand Pear (Pyrus pyrifolia Nakai)

MicroRNA156 is an evolutionarily highly conserved plant micro-RNA (miRNA) that controls an age-dependent flowering pathway. miR156 and its target SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) genes regulate anthocyanin accumulation in plants, but it is unknown whether this process is affected by light. Red Chinese sand pear (Pyrus pyrifolia) fruits exhibit a unique coloration pattern in response to bagging treatments, which makes them appropriate for studying the molecular mechanism underlying light-induced anthocyanin accumulation in fruit. Based on high-throughput miRNA and degradome sequencing data, we determined that miR156 was expressed in pear fruit peels, and targeted four SPL genes. Light-responsive elements were detected in the promoter regions of the miR156a and miR156ba precursors. We identified 19 SPL genes using the “Suli” pear (Pyrus pyrifolia Chinese White Pear Group) genome database, of which seven members were putative miR156 targets. The upregulated expression of anthocyanin biosynthetic and regulatory genes and downregulated expression of PpSPL2, PpSPL5, PpSPL7, PpSPL9, PpSPL10, PpSPL13, PpSPL16, PpSPL17, and PpSPL18 were observed in pear fruits after bags were removed from plants during the anthocyanin accumulation period. Additionally, miR156a/ba/g/s/sa abundance increased after bags were removed. Yeast two-hybrid results suggested that PpMYB10, PpbHLH, and PpWD40 could form a protein complex, probably involved in anthocyanin biosynthesis. Additionally, PpSPL10 and PpSPL13 interacted with PpMYB10. The results obtained in this study are helpful in understanding the possible role of miR156 and its target PpSPL genes in regulating light-induced red peel coloration and anthocyanin accumulation in pear.

INTRODUCTION miR156 is one of the most abundant and evolutionarily conserved micro-RNAs (miRNAs) in plants, and affects diverse aspects of plant growth and development by regulating SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) genes (Poethig, 2010;Huijser and Schmid, 2011). miR156 and the SPL genes are down-and up-regulated as plants develop, respectively (Wu and Poethig, 2006). In Arabidopsis thaliana, 11 out of the 17 SPL genes are targeted by miR156, including SPL3, SPL4, and SPL5, which form one group. Constitutive overexpression of miR156-resistant SPL3 (rSPL3), rSPL4, or rSPL5 results in early flowering (Cardon et al., 1997;Wu and Poethig, 2006). Additionally, SPL3 is a direct upstream activator of three transcription factors (i.e., LFY, FUL, and AP1) that are floral meristem identity genes (Yamaguchi et al., 2009). A second group includes SPL9 and SPL15, which promote the juvenileto-adult phase transition (Schwarz et al., 2008;Wang et al., 2008). The expression of a miR156-resistant SPL9 (rSPL9) causes the extremely early flowering in plants Wu et al., 2009). SPL9 and SPL15 also influence the characteristic appearance of abaxial trichomes on adult leaves , and ensure plants are capable of flowering in response to an appropriate photoperiod (Schwarz et al., 2008). A third group of SPL genes includes SPL2, SPL10, and SPL11, which have a minor role in vegetative phase changes in addition to their major role in early embryonic patterning via a redundant regulation of the embryonic morphogenesis-to-maturation phase transition (Nodine and Bartel, 2010).
Besides endogenous developmental cues, miR156-SPL module could also perceive various external stimuli as input signals for developmental and morphological adaptive changes . In Arabidopsis, high level of CO2 accelerates flowering probably by decreasing and increasing the level of miR156 and SPLs, respectively (May et al., 2013). Ambient temperature is also an input signal for miR156-SPL model as overexpression of miR156 causes delayed flowering at lower ambient temperature (16 • C), associated with the downregulation of SPL3, FT, and FRUITFULL expression (Kim et al., 2012). Additionally, miR156-SPL module is regulated by abiotic stresses, including Phosphate Deficiency (Hsieh et al., 2009), salt stress (Ding et al., 2009), and drought stress (Sun et al., 2012).
The miR156-SPL module was recently reported to be involved in regulating anthocyanin biosynthesis in plants (Gou et al., 2011;Cui et al., 2014). Anthocyanin is the main pigment in flowers and fruits, and is responsible for characteristic reddish, bluish, and purple hues. Anthocyanin biosynthesis involves several structural genes, including PAL (phenylalanine ammonia lyase), CHS (chalcone synthase), CHI (chalcone isomerase), F3H (flavanone 3-hydroxylase), DFR (dihydroflavonol 4-reductase), ANS (anthocyanidin synthase), and UFGT (UDP-glucose:flavonoid 3-O-glucosyltransferase) (Winkel-Shirley, 2001). Anthocyanin biosynthesis is regulated by numerous transcription factors, including MYB, bHLH, and WD40, which form a MYB/bHLH/WD40 complex that mediates the expression of anthocyanin biosynthetic genes in diverse species including Arabidopsis (Broun, 2005), Chinese bayberry , and strawberry (Schaart et al., 2013). Regulation of anthoycanin accumulation has been well studied in fruit crops. In litchi, UFGT and LcMYB1 are the two key genes regulating anthocyanin biosynthesis and finally determining the pericarp color Lai et al., 2014). Delphinidin-derived anthocyanin is accumulated in the fruits of Lycium ruthenicum but not in the fruits of Lycium barbarum, which might be controlled by the expression patterns of both regulatory and structural genes and the transcriptional ratio of branch-node structural genes F3 ′ 5 ′ H/F3 ′ H (Zeng et al., 2014). In apple, a rearrangement in the upstream regulatory region of the MYB10 is responsible for increasing the level of anthocyanin throughout the plant to produce a phenotype with red foliage and red fruit flesh (Espley et al., 2009). Anthocyanin biosynthesis is also regulated by other factors, such as sugar. In Arabidopsis, sucrose treatment could strongly up-regulated the flavonoid and anthocyanin biosynthetic pathways, thus affecting flavonoid and anthocyanin content (Solfanelli et al., 2006). In radish, exogenous glucose, fructose, or sucrose promotes anthocyanin accumulation in hypocotyl by inducing expression of the biosynthetic genes including CHS and ANS (Hara et al., 2003). Analyses of wildtype, miR156-upregulated (Pro35S:MIR156), and miR156-downregulated (Pro35S:MIM156) A. thaliana lines revealed that increased miR156 abundance promotes anthocyanin accumulation. Additionally, at least one of the miR156 target genes (i.e., SPL9) encodes a protein which can negatively regulate anthocyanin biosynthesis through disrupting the MYB/bHLH/WD40 complex by competing with bHLH for the binding of MYB (Gou et al., 2011). Under salt and mannitol stress conditions, SPL9 decreases anthocyanin accumulation by repressing the expression of DFR (Cui et al., 2014). In the blood-fleshed peach cultivar "Dahongpao" and the whitefleshed cultivar "Baifeng", SPL1 represses the transactivation activity of the BL-PpNAC1 heterodimer, which subsequently inhibits PpMYB10 expression, leading to limited accumulation of anthocyanins (Zhou et al., 2015). Additionally, miR172 abundance is regulated by miR156 and its targets SPL genes Jung et al., 2011), and miR172 can be directly activated by SPL9 and SPL10 , indicating the possibility of miR172 participating in regulating anthocyanin biosynthesis in plants. These results reveal that the regulation of miR156/SPL modules during anthocyanin biosynthesis is affected by genotype and abiotic stress. There is currently little information regarding whether this regulatory pathway is influenced by light, which is one of the most important environmental factors required for anthocyanin accumulation in various plant species. COP1 and HY5 are the two key proteins involved in light-induced anthocyanin accumulation in fruit. COP1 is an ubiquitin E3 ligases, which is necessary for the ubiquitination and degradation of MdMYB1 protein in the dark and therefore negatively regulates light-induced anthocyanin biosynthesis in apple (Li et al., 2012). However, under UV-B irradiation, the expression of MdCOP1 is induced, which activates MdHY5 signaling by binding to the promoter regions of MdMYBs, and finally leads to the red coloration of apple peels (Peng et al., 2012). However, besides light signal transduction pathway, more other pathways which regulate light-induced anthocyanin biosynthesis in fruit need to be defined.
Pear (Pyrus spp.) is one of the most dominant temperate fruit crops worldwide. Anthocyanin accumulation in pear is considerably different than that of other important fruit species, including grape and apple. Grape rapidly softens at the beginning of the second growth cycle (the growth curve of grape is double-sigmoidal pattern), many weeks before ripeness, and the onset of rapid softening is called "veraison" by viticulturists. After varaison, red colored grape cultivars start accumulating anthocyanin, and the skin color is changed from green to red (Coombe, 1973). Apple shows two peaks of anthocyanin concentration during development: the fruitlet of apple is red, then the red color disappears, and anthocyanin is accumulated again and peaked toward maturity (Steyn et al., 2005). For pear, the fruit pigmentation patterns differ between the two main commercial pear species, Chinese sand pears (Pyrus pyrifolia Nakai) and European pears (Pyrus communis L.). Unlike red Chinese sand pears, in which anthocyanin accumulation peaks in mature fruit (Huang et al., 2009), European pears exhibit a decreased red coloration as fruits mature. The highest anthocyanin levels in European pears occur about midway between anthesis and harvest (Steyn et al., 2005). Chinese sand pears and European pears also respond differently to bagging treatments and postharvest UV-B/visible light irradiation. Fruit bagging and postharvest UV-B/visible light irradiation are considered attached and detached treatments to induce anthocyanin biosynthesis and improve fruit color. Fruit bagging has been used in apples (Arakawa, 1988;Bai et al., 2016) and Asian pears (Huang et al., 2009;Qian et al., 2013). Postharvest irradiation has been widely applied to apples (Ubi et al., 2006), Asian pears Sun et al., 2014), grapes (Li et al., 2009), and cherries (Kataoka et al., 2005). Ten days after bags are removed from plants or 10 days after the start of postharvest irradiation, Chinese sand pears appear red (i.e., accumulated anthocyanin), while anthocyanin is almost undetectable in European pears . The distinct pigmentation pattern in response to bagging treatments makes Chinese sand pear appropriate for studying the molecular mechanism regulating light-induced anthocyanin accumulation in fruits. However, it is currently unknown whether miR156/SPL modules affect this process in pear plants.
Different family members of miR156 and SPL genes have been identified (Cardon et al., 1999;Yu et al., 2015), and the abundance of miR156 exhibits tissue specificity (Xia et al., 2012;Yu et al., 2015). Pear miRNAs have been identified in fruit flesh and buds by high-throughput sequencing Niu et al., 2016). In this study, we used miRNA sequencing, degradome analysis, and a genome-wide identification of SPL gene family members, to select miR156 and SPL members that might participate in the light-induced anthocyanin biosynthesis in pear fruit peel. Furthermore, we searched for light-responsive elements in the promoter region of MIR156 genes (i.e., pre-miR156). We also examined the expression patterns of these potential miR156 and SPL genes after bagging treatments, and investigated the interaction between SPLs and the MYB/bHLH/WD40 transcription factor complex.

Plant Materials, Treatments, and Total RNA Extraction
Red Chinese sand pear (cultivar "Meirensu") fruits were obtained from a commercial orchard in Kunming, Yunnan, China. Ten mature trees that were similar in size and uniformly exposed to sunlight were selected. We bagged 80 fruitlets (the vertical diameter of fruit was around 2 cm) per tree 40 days after full bloom with double-layered yellow-black paper bags [Kobayashi (Qingdao) Co., Ltd., China]. Half of the fruits were re-exposed to sunlight after the bags were removed 10 days before harvest (the total bagging period lasted for about 3 months). The fruits that remained bagged served as control samples. Fruits were harvested 0, 1, 3, and 6 days after bag removal (DABR), stored on dry ice, and transported to the lab as quickly as possible. The skin from 10 individual fruits were pooled, instantly frozen in liquid nitrogen, and stored at −80 • C for assays.

Small-RNA Library Construction, Sequencing, and Data Analysis
Total RNA from re-exposed and control fruit peels collected at 0, 1, 3, and 6 DABR was isolated using a CTAB method (Qian et al., 2014a). After ligation of 3 ′ -and 5 ′ -adapters, and RT-PCR amplification, small RNA (sRNA) library was size fractioned by PAGE gel. A mixed sRNA library was constructed and sequenced using the Illumina Genome Analyzer IIx (Illumina, Inc., Santa Clara, CA, USA). After filtering and adapter cutting, the clean tags were used by BLASTall to search the GenBank and Rfam 10.0 databases (Kozomara and Griffiths-Jones, 2014) to detect rRNAs, scRNAs, snoRNAs, snRNAs, and tRNAs, which were subsequently removed. The remaining sRNA tags were aligned with mRNA sequences to identify any degraded mRNA fragments, which were eliminated (http://peargenome.njau.edu. cn/default.asp?d=1&m=1) . The potential miRNAs were then mapped to the pear genome sequence (http:// peargenome.njau.edu.cn:8004/default.asp?d=4&m=2)  using SOAP 2.20 (http://soap.genomics.org.cn/soapsplice. html), and their distribution on the genome was analyzed. Only sRNA tags that formed good stem-loop structures and had a miRNA/miRNA * pair were considered as potential miRNAs. The secondary structure was predicted by RNAfold, and novel miRNAs were identified by MIREAP. The criteria used to identify candidate miRNAs were as previously described (Niu et al., 2013). Analyses of the high-throughput sequencing profiles were based on the number of reads. The miRNA expression levels were converted to transcripts per million normalized values as follows: normalized expression = actual miRNA count/ total number of clean reads × 1,000,000. Raw data from sRNA sequencing was submitted to Sequence Read Archive (SRA ID: SRR5074534).

Degradome Library Construction, Data Analysis, and Target Identification
Equal amounts of RNA from seven independent pear fruit peel libraries (re-exposed and control fruit peels collected at 0, 1, 3, and 6 DABR) were pooled to construct a degradome library. After 5 ′ adapter contaminants were eliminated and genome mapping was completed as described for the sRNA data, the degradome sequencing data were analyzed using the CleaveLand (version 3.0) (Addo-Quaye et al., 2009) program. The alignment score threshold was set to 4.5 for conserved and moderately conserved miRNAs and to 5 for novel and candidate miRNAs (Xia et al., 2012). The apple consensus gene set from AppleGFDB and the annotation information for miRNA target genes were retrieved from the Genome Database for Rosaceae . Degradome data were normalized to transcripts per million values. Based on the number of degradome sequences and their abundance values, the cleavage sites of miRNA targets were classified into 5 categories as previous study (Xia et al., 2012): Category 0: > 1 raw read at the position, read abundance at position is equal to the maximum read abundance on the transcript, and there is only one maximum on the transcript. Category 1: > 1 raw read at the position, read abundance at position is equal to the maximum read abundance on the transcript, and there is more than one maximum position on the transcript. Category 2: > 1 raw read at the position, read abundance at position is less than the maximum but higher than the median read abundance for the transcript. Category 3: > 1 raw read at the position, read abundance at position is equal to or less than the median read abundance for the transcript. Category 4: only 1 raw read at the position. Raw data from degradome sequencing was submitted to Sequence Read Archive (SRA ID: SRR5074535).

′ -Rapid Amplification of cDNA Ends
We completed 5 ′ -rapid amplification of cDNA ends (RACE) using the SMARTer RACE cDNA Amplification Kit (Clontech, Palo Alto, CA, USA) to verify the miR156-mediated cleavage events identified during degradome sequencing. A 2-µg sample of mixed total RNA isolated from seven independent pear fruit peels (re-exposed and control fruit peels collected at 0, 1, 3, and 6 DABR) was ligated with 5 ′ RNA adapters at room temperature. All polymerase chain reaction (PCR) products were inserted into the pMD18-T vector (Takara, Dalian, China) and sequenced. Specific primers were designed for nested PCR (Supplementary Table S1 in Supplementary File 1).

Gene Ontology Analysis
An additional analysis involving the Gene Ontology (GO) databases was completed to more thoroughly characterize the functions of miRNA targets. Blast2GO was used to store and manage information from the GO databases (https:// www.blast2go.com). All targets were identified using BLASTx searches against the GO protein database. The GO analyses were completed using a query search of the GO databases.

Identification and Annotation of Pear SQUAMOSA Promoter Binding Protein (SBP)-Box Genes
The hidden Markov model (HMM) profile of the SBP domain (Accession no. PF03110) was downloaded from the Pfam database (http://www.sanger.ac.uk). This domain was then used as a query for BLASTP searches of the GenBank non-redundant protein database and the Pear Genome Database (http:// peargenome.njau.edu.cn:8004/default.asp?d=4&m=2) . All hits with an E < 0.01 were collected. The nonredundant protein sequences encoded by putative PpSPL genes were manually checked for the presence of the SBP domain, and genes without a complete SBP domain were eliminated. All putative PpSPL genes were confirmed again by cloning and sequencing, using a mixed cDNA template obtained from different tissues of "Meirensu" plants, including young fruits, young leaves, petals, young stems, and the skin of mature fruit. Primer details are provided in Supplementary Table S1 in Supplementary File 1. The PpSPL genes were named according to their gene ID order in the pear genome.

Sequence Alignments and Phylogenetic Analyses
Multiple sequences were aligned using DNAMAN software (version 5.2). The sequence logo was obtained using the online WebLogo platform (http://weblogo.berkeley.edu). Phylogenetic trees were constructed using the neighbor-joining method of the MEGA 4.0 program (Tamura et al., 2007), and the bootstrap test was replicated 1,000 times. The miR156 target sites were identified using sequence alignments and manual analyses.

Quantitative Real-Time Polymerase Chain Reaction Analysis
Total RNA was extracted using a modified CTAB method (Qian et al., 2014a). The total RNA sample was treated with DNase I to remove any contaminating genomic DNA, and then quantified. First-strand cDNA was synthesized from 4 µg DNAfree RNA using a Revert Aid TM First-Strand cDNA Synthesis Kit (Fermentas, Glen Burnie, MD, USA). A 2-µL sample of 10fold diluted cDNA was used as the template for gene cloning and quantitative real-time polymerase chain reaction (qPCR) analyses.
The qPCR mixture contained 10.0 µL SYBR Premix Ex Taq TM (Takara, Ohtsu, Japan), 0.4 µL each primer (10 µM), 2 µL cDNA, and 7.2 µL RNase-free water in a total volume of 20 µL. The PCR was completed using a LightCycler 1.5 instrument (Roche, Germany) with the following program: 95 • C for 30 s; 40 cycles of 95 • C for 5 s and 60 • C for 20 s. A template-free control was included for each primer pair. The qPCR primers (Supplementary Tables S2, S3 in Supplementary File 1) were designed using Primer 3 software (http://frodo.wi.mit.edu/cgibin/primer3/primer3_www.cgi). All qPCR data were normalized using the threshold cycle value for the Pyrus species actin gene (PyActin, JN684184). Each sample was analyzed three times.

Quantitative Real-Time Polymerase Chain Reaction Analysis of miR156 in Pear
Total RNA was extracted using the pBiozol Total RNA Extraction Reagent (BioFlux, China). According to the instructions provided for the miRNA cDNA Synthesis Kit (Takara), 1 µg total RNA was polyadenylated by poly(A) polymerase. The poly(A)-tailed total RNA was reverse-transcribed by PrimeScript R RTase with a universal adapter primer (containing oligo-dT). The qPCR analysis was completed using SYBR R Premix Ex Taq II (Perfect Real Time) (Takara) and a LightCycler 1.5 instrument. Each sample was analyzed three times, and 5S rRNA was used as the internal control gene Niu et al., 2016). The qPCR primers are listed in Supplementary Table S4 in Supplementary File 1.

Yeast Two-Hybrid Assay
A yeast two-hybrid assay was conducted using the Matchmaker TM Gold Yeast Two-Hybrid System (Clontech). We used the pGADT7 vector, which contains the GAL4 activation domain, and the pGBKT7 vector, which carries the GAL4 DNA binding domain. The gene open reading frames were inserted into the multiple cloning sites of pGBKT7 and pGADT7 using the In-Fusion R HD Cloning Kit (Clontech) to obtain the gene-BD and gene-AD plasmids, respectively. Details of the primers used during the yeast two-hybrid assay to obtain the open reading frames are listed in Supplementary Table S5 in Supplementary File 1. The gene-BD plasmid was inserted alone or together with the gene-AD plasmid into the Y2HGold yeast strain according to a polyethylene glycol/lithium acetate method. The autoactivation of the colonies transformed with the gene-BD plasmid was tested on synthetically defined (SD) medium lacking tryptophan, but supplemented with 5-bromo-4-chloro-3-indolyl α-D-galactopyranoside (X-α-Gal) and aureobasidin A (AbA). The co-transformed colonies were selected on SD medium lacking leucine and tryptophan (i.e., double dropout medium), and screened for growth on quadruple dropout SD medium (i.e., lacking adenine, histidine, leucine, and tryptophan) supplemented with X-α-Gal and AbA to verify positive interactions.

Statistical Analysis
After running ANOVA, LSDs (α = 0.05) were calculated for mean separations of anthocyanin concentration using the Data Processing System (DPS, version 3.01, Zhejiang University, Hangzhou, China). Differences of gene expression were statistically evaluated by the analysis of variance and Tukey's test using SPSS 19.0 (SPSS, Chicago, IL, USA). Probability values of <0.05 were considered statistically significant.

Genome-Wide Identification of Pear miRNAs and Their Targets, and Analysis of Promoter Elements of miR156 Precursors
A total of 13,954,538 raw reads were acquired through highthroughput sequencing of sRNA from pear fruit peels of bagged plants. Corrupted adapter sequences (e.g., missing 3 ′ adapter), sequences shorter than 17 bases or 25 bases after removing the 3 ′ adapter, and junk reads were eliminated. The remaining 9,915,560 clean reads were used for predicting miRNAs (Supplementary Table S6 in Supplementary File 1). A total of 167 conserved, 34 moderately conserved, and 21 novel miRNAs were identified in the sRNA library (Supplementary Tables S7, S8 in Supplementary File 1). Among the detected miRNAs, 26 miR156 members were identified, including some that were highly expressed in fruit peels such as miR156c (424 reads). Some members were detected at very low levels, including miR156sd (1 read) and miR156bd-3p (1 read). In addition, all the novel miRNAs were highly expressed in fruit peel, with the raw reads from 125 to 1,866 (Supplementary Tables S8 in Supplementary  File 1).
Due to the results from sRNA and degradome sequencing, miR156a and miR156ba were chosen as two examples of miR156s, and the genomic region 1,500 bp upstream of pre-miR156a and pre-miR156ba was considered the promoter region and was analyzed using the PlantCARE online program (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). We identified 16 light-responsive elements in the pre-miR156a promoter (i.e., two Box 4 elements, six G-boxes, one GA-motif, one GAG-motif, two GT1-motifs, one I-box, one MRE, and two sp1 elements) ( Table 1). The pre-miR156ba promoter contained one AE-box, two Box 4 elements, seven G-boxes, one GAG-motif, three GT1-motifs, one TCCC, and two sp1 elements ( Table 1).

Genome-Wide Identification of PpSPL Genes
We isolated 19 pear PpSPL genes, which were named based on the gene ID order in the pear genome (i.e., PpSPL1-19) (Supplementary Table S11 in Supplementary File 1). All PpSPL family members shared a conserved SBP domain, with two zinc finger-like structures (i.e., Zn-1 and Zn-2) and a highly conserved bipartite nuclear localization signal (NLS; Figure 1). To further study the evolutionary relationships among PpSPL genes and SPL genes from other plant species, we prepared a dataset of 95 SPL sequences from apple, grape, rice, tomato, and A. thaliana. The amino acid sequences of the SBP domains (75 amino acids) encoded by the SPL genes were used to construct a phylogenetic tree with the neighbor-joining algorithm (Figure 2). The plant SBP domains were clustered into seven distinct groups (i.e., G1-G7). The PpSPL family members were distributed over all seven cluster groups, but were most closely associated with apple genes. Additionally, PpSPL9 and PpSPL13 were homologs of A. thaliana AtSPL9, which negatively regulates anthocyanin biosynthesis.
Of the miR156 members, miR156a/g/s/sa and miR156ba shared highly similar sequences. Therefore, a conserved primer was designed for analyzing the total expression of miR156a/g/s/sa and miR156ba. After bags were removed from plants, miR156a/ba/g/s/sa transcript abundance increased, peaking at 1-3 DABR (i.e., >3-fold increase) followed by a gradual decrease, but the expression level was still relatively high at 6 DABR (Figure 6). In contrast, miR156a/ba/g/s/sa transcript abundance in the control samples was low at all time points.

Interactions among Transcription Factors Related to Anthocyanin Biosynthesis
Because anthocyanin biosynthesis involves three transcription factors (i.e., MYB, bHLH, and WD40), the interaction among these proteins was tested with a yeast two-hybrid assay. The PpMYB10, PpbHLH3, PpbHLH33, and PpWD40 open reading frames were inserted into the pGBKT7 and pGADT7 vectors to obtain BD-bait and AD-prey protein complexes. We used SD medium lacking tryptophan, but supplemented with X-α-Gal and AbA, to verify the autoactivation of BD-bait protein complexes. PpMYB10-BD, PpbHLH3-BD, and PpWD40-BD were not autoactivated on medium containing 125 ng/ml AbA. However, PpbHLH33-BD was autoactivated on medium with 1,000 ng/ml AbA (data not shown), so it could not be analyzed by the yeast two-hybrid assay. Our results indicated that PpMYB10 and PpWD40 did not interact with each other, while both transcription factors interacted with PpbHLH3 and PpbHLH33 (Figure 8). These findings suggested that a MYB10/bHLH/WD40 protein complex exists in pear.

Interactions among PpSPL Proteins and Transcription Factors Related to Anthocyanin Biosynthesis
In this study, DFR expression was slightly upregulated (<2-fold increase) during anthocyanin accumulation (Figure 5), and NAC genes were not considered. Thus, proteins encoded by genes that were downregulated after bags were removed from plants (e.g., PpSPL2, PpSPL5, PpSPL7, PpSPL9, PpSPL10, PpSPL13, PpSPL16, and PpSPL17) were fused to BD, while PpMYB10, PpbHLH3, and PpbHLH33 were fused to AD to test potential interactions in a yeast two-hybrid assay. Unfortunately, PpSPL2-BD, PpSPL9-BD, and PpSPL18-BD were autoactivated (data not shown). Our results indicated that PpSPL10 and PpSPL13 interacted with PpMYB10 (Figure 9), while there was no interaction between PpSPL and PpbHLH (Supplementary Figure 1).

DISCUSSION 11 PpSPLs Were Putative miR156 Target Genes in Pear
Different family members of miR156 and SPL genes have been isolated in plants (Cardon et al., 1999;Yu et al., 2015). There are eight miR156 genes in the A. thaliana genome (i.e., miR156A, miR156B, miR156C, miR156D, miR156E, miR156F, miR156G, and miR156H). The expression patterns of these genes were highly tissue-specific. For example, miR156A was highly expressed in the first two pairs of leaves in seedlings, miR156B was specifically expressed in the central cylinder of the primary root primordial cells, and miR156D was highly expressed in seedlings (Yu et al., 2015). This variability in expression patterns suggests that the miR156 genes have diverse roles influencing plant growth and development.
Because of their crucial contribution to the regulation of a broad range of developmental processes, SPL genes have been subjected to genome-wide identification and analysis in diverse plant species. Consequently,17,27,19,19, and 15 SPL gene family members have been identified in A. thaliana (Cardon et al., 1999), apple , grapevine , rice (Xie et al., 2006), and tomato (Salinas et al., 2012), respectively. Of the known SPL genes, only SPL9 and SPL10 in A. thaliana and SPL1 in peach have been shown to participate in regulating anthocyanin biosynthesis (Gou et al., 2011;Zhou et al., 2015). There is currently no information regarding if specific miR156 and SPL gene family members are involved in regulating anthocyanin accumulation in pear.
In this study, we identified 26 miR156 members (Supplementary Tables S7 in Supplementary File 1), and 19 pear PpSPL genes (Supplementary Table S11 in Supplementary File 1). According to the phylogenetic tree, the SBP domains were clustered into seven distinct groups (i.e., G1-G7), which was similar to the results of a previous study on apple . The PpSPL family members were more closely related to the SPL genes from eudicots, including apple, grape, tomato, and A. thaliana, than the rice SPL genes (i.e., monocotyledonous species). These findings imply that SPL genes in eudicots diverged more recently from a common ancestor than the SPL genes of monocots. In addition, 11 PpSPL genes are potentially regulated by miR156 (Supplementary Tables S10 in Supplementary File 1;  Figure 3), which was similar to apple, where more than half of the MdSBP genes (15 out of 27) are targeted by miR156 .
Increased Expression Levels of PpMYB10, PpCHS, and PpUFGT Were Related to Light-Induced Anthocyanin Accumulation in Pear During the accumulation of anthocyanin after bags were removed from plants, the expression of most anthocyanin biosynthetic genes were upregulated. PpCHS and PpUFGT expression was especially upregulated at 3 DABR, with about 200-and 600-fold increases over expression levels at 0 DABR ( Figure 5). Chalcone synthase catalyzes the first committed step of flavonoid and anthocyanin biosynthesis, and the CHS gene is considered to be a UV-B-responsive gene. Additionally, it is targeted by HY5, which is a bZIP transcription factor that mediates UV-B-induced gene expression changes downstream of UVR8 (Binkert et al., 2014). In our previous studies, CHS expression was induced by light  and methyl jasmonate (Qian et al., 2014b). Similar results were reported for other species, including apple (Ubi et al., 2006), Chinese bay berry (Niu et al., 2010), and Petunia hybrida (Koes et al., 1989). The UFGT enzyme catalyzes the production of stable anthocyanin by adding sugar moieties to the anthocyanin aglycone. The transcription of UFGT is one of the most rate-limiting steps in the anthocyanin biosynthesis pathway (Montefiori et al., 2011). In white grapes, the loss of fruit skin coloration is due to the absence of UFGT transcription (Boss et al., 1996). In litchi (Litchi chinensis), of the six genes encoding proteins associated FIGURE 2 | Phylogenetic analysis of SPLs from pear and other plant species. A phylogenetic tree was constructed using SBP domain protein sequences from pear (PpSPL, marked with black dots), apple (MdSBP, marked with white dots), grape (VvSBP, marked with white squares), rice (OsSPL, marked with white rhombuses), tomato (CNR or SlySBP, marked with white triangles), and Arabidopsis thaliana (AtSPL, marked with black triangles). The symbols for Arabidopsis thaliana and tomato genes are inverted in the right and left halves of the figure. PpSPL1, PpSPL4,PpSPL5,PpSPL6,PpSPL9,PpSPL10,PpSPL12,PpSPL13,PpSPL14,PpSPL15,and PpSPL18 are putative miR156 targets. The SBP domain sequences used to construct the phylogenetic tree as well as the accession numbers or locus IDs and data sources of the SBP-box genes from different species are listed in Supplementary Table S12 in Supplementary File 1.
with anthocyanin biosynthesis, the UFGT genes were the only ones whose transcript levels were significantly correlated with the pericarp anthocyanin content (Wei et al., 2011). These observations suggest that UFGT is required for anthocyanin accumulation.
Among transcription factors, the PpMYB10 and PpWD40 expression levels increased significantly, with peaks at 1 DABR ( Figure 5). Additionally, the R2R3 MYB class of transcription factors is the dominant regulator of anthocyanin biosynthesis in several plant species, including pear (Feng et al., 2010), apple (Espley et al., 2009), grape (Kobayashi et al., 2004), and Chinese bayberry (Niu et al., 2010). Specific R2R3 MYBs were reported to activate the promoter regions of anthocyanin biosynthetic genes, such as DFR (Niu et al., 2010) and UFGT , thereby promoting anthocyanin accumulation. In Chinese sand pear, PpMYB10 expression was considerably upregulated during anthocyanin biosynthesis in different developmental stages (Feng et al., 2010), following methyl jasmonate treatments (Qian et al., 2014b), and after postharvest irradiation . These findings, together with the results of the present study, indicate that upregulated PpMYB10 expression is related to light-induced anthocyanin accumulation in pear. Interestingly, almost no anthocyanin concentration is detected in European pear under bagging treatment, and the expression of PpMYB10 is not up-regulated as well , indicating that bagging treatment doesn't activate the light transduction pathway in European pear, or probably the regulation of PpMYB10 expression from the down-stream transcription factors of light transduction pathway (such as HY5) is blocked due to certain factors. Thus, the transcription of PpMYB10 is not triggered and the anthocyanin biosynthesis pathway is not activated.  FIGURE 5 | Expression of anthocyanin biosynthetic genes and regulatory genes at different stages after bags were removed from plants. Each value represents the mean ± standard error of three biological replicates.

miR156-SPL Module Also Responded during the Light-Induced Anthocyanin Biosynthesis in Pear
During the biosynthesis of anthocyanin after bags were removed from plants, miR156 a/ba/g/s/sa abundance increased and peaked at 1-3 DABR (Figure 6). In contrast, we observed a decrease in the transcription levels of PpSPL2, PpSPL5, PpSPL7, PpSPL9, PpSPL10, PpSPL13, PpSPL16, PpSPL17, and PpSPL18 (Figure 7). Additionally, PpSPL9 and PpSPL13 were determined to be homologs of A. thaliana AtSPL9, which regulates anthocyanin biosynthesis (Gou et al., 2011). In A. thaliana, increased miR156 activity promotes anthocyanin accumulation, whereas SPL activity is negatively correlated with anthocyanin content in wide-type plants and transgenic lines expressing Pro35S:MIR156 and Pro35S:MIM156 (Gou et al., 2011). Furthermore, during anthocyanin accumulation in A. thaliana plants treated with salt and mannitol, miR156 expression levels increased, while SPL expression was downregulated (Cui et al., 2014). In FIGURE 6 | Expression of miR156a/ba/g/s/sa at different stages after bags were removed from plants. Each value represents the mean ± standard error of three biological replicates.
the red-fleshed peach cultivar "Dahongpao", SPL1 was highly expressed in white-fleshed young fruits, but expressed at low levels in red-fleshed ripening fruits (Zhou et al., 2015). These FIGURE 7 | Expression of PpSPL genes at different stages after bags were removed from plants. Each value represents the mean ± standard error of three biological replicates.
FIGURE 8 | Yeast two-hybrid analysis of the physical interactions among transcription factors related to pear anthocyanin biosynthesis. DDO: double dropout synthetically defined (SD) medium (i.e., lacking leucine and tryptophan). QDO/X/A: quadruple dropout SD medium containing X-α-Gal and aureobasidin A (AbA) (i.e., SD medium lacking adenine, histidine, leucine, and tryptophan, but supplemented with X-α-Gal and AbA). The empty pGADT7 vector was used as a negative control. results indicate that miR156 expression is closely related to anthocyanin accumulation, while the expression of the target SPL genes is negatively correlated with anthocyanin biosynthesis. Concerned about the patterns anthocyanin accumulation and gene expression, we hypothesized the procedure of light-induced anthocyanin biosynthesis from up-stream to down-stream: Exposure to light increases the abundance of miR156 (peaked at 1 DABR), which degrades SPL genes (down-regulated at 1 DABR), resulting in the up-regulation of PpMYB10 (peaked at 1 DABR) and formation of the MYB/bHLH/WD40 complex, subsequently activates the expression of most structural genes including PpCHS, PpCHI, PpF3H, PpANS, and PpUFGT (peaked at 3 DABR). More proteins encoded by these genes are translated, catalyze the reactions of anthocyanin biosynthesis pathway, and finally lead to the accumulation of anthocyanin from 3 DABR to 6 DABR.

PpMYB, PpbHLH, and PpWD40 Could form a Protein Complex
The regulation of anthocyanin biosynthesis involves numerous transcription factors, including MYB, bHLH, and WD40, which form the MYB/bHLH/WD40 complex that mediates the expression of anthocyanin biosynthetic genes in diverse species, including A. thaliana (Broun, 2005), Chinese bayberry , and strawberry (Schaart et al., 2013). Among these transcription factors, bHLH acts as a bridge interacting with MYB and WD40 to form the complex. An earlier study involving apple revealed that there is no direct interaction between MYB and WD40 proteins . To the best of our knowledge, the current study is the first to examine the MYB/bHLH/WD40 complex in pear. Our yeast two-hybrid results indicated that PpMYB10 and PpWD40 interacted with PpbHLH3 and PpbHLH33, but not with each other. PpbHLH3 and PpbHLH33 were observed to interact with each other (Figure 8). Our findings imply that the MYB10/bHLH/WD40 transcription factor complex is conserved among plant species.

PpSPL10 and PpSPL13 Could Form a Heterodimer with PpMYB10
Our results revealed that PpSPL10 and PpSPL13 (orthologs of A. thaliana AtSPL9) interacted with PpMYB10 (Figure 9). We also observed that DFR expression was slightly upregulated (<2-fold increase) during anthocyanin accumulation (Figure 5). Therefore, the miR156-SPL9-DFR pathway, which responds to abiotic stresses in A. thaliana (Cui et al., 2014), may not be responsible for the light-induced accumulation of anthocyanin in pear. According to the previous results in A. thaliana (Gou et al., 2011), we hypothesize the involvement of miR156-SPL module in light-induced pear peel coloration and anthocyanin accumulation probably through the interaction between SPL and MYB10, which is the crucial part of MYB/bHLH/WD40 complex.
In conclusion, after bag removal, during the accumulation of anthocyanin in pear and up-regulation of genes related to anthocyanin biosynthesis, the expression of miR156 and its target PpSPL genes, including PpSPL2, PpSPL5, PpSPL7, PpSPL9, PpSPL10, PpSPL13, PpSPL16, PpSPL17, and PpSPL18 were increased and decreased, respectively. In addition, three transcription factors MYB10, bHLH, and WD40 could form a protein complex, and PpSPL10 and PpSPL13 could also interact with PpMYB10. Our results provide some clues about the miR156-SPL module probably participating in light-induced red peel coloration in pear.

AUTHOR CONTRIBUTIONS
MQ, JN, YT, and DZ designed the study; MQ, JN, and QN conducted the experiments; MQ, JN, QN, SB, LB, JL, YS, YT, and DZ analyzed the data; MQ, YT, and DZ wrote the article.