Plant materials

Three tetraploid wheat (Triticum turgidum L.) NILs with different spike phenotypes, including bh50 with NSs, bh51 with FRSs, and bh53 with RSs, were used for proteomic analysis (Fig. 1). The NILs were planted in a greenhouse, and immature spike tissue was harvested from the main tiller of each plant at the spikelet initiation stage (following the double ridge stage, Fig. 1).

Protein extraction

Total protein from immature spikelet tissue was extracted when the spikelets were approximately 3 mm long using a Plant Total Protein Extraction Kit (Sigma-Aldrich, St. Louis, MO, USA), in accordance with the manufacturer’s instructions. In total, 250 mg of plant tissue was ground in liquid nitrogen and extracted using phenol, in accordance with the method described by Isaacson et al. (2006) with minor modifications [18]. The mixture was purified using acetone, and the purified proteins were dissolved in lysis solution at 30 °C for 1 h. After centrifugation, the concentration of extracted protein was determined using the Bradford method (Bio-Rad, Hercules, CA, USA, and proteins were stored at − 80 °C for iTRAQ analysis [19].

Protein quantitation and iTRAQ labeling

Samples containing approximately 10 μg of protein were run on 12% SDS-PAGE gels and were visualized using Coomassie brilliant blue stain, in accordance with Candiano’s protocol [20]. The stained gels were scanned using an Image Scanner (GE Healthcare, Chicago, IL, USA) at a resolution of 300 dots per inch. Aliquots containing 100 μg were obtained from each sample by adding 5 volumes of cold acetone. A total of 50 μL of dissolution buffer (4% SDS, 100 mM DTT, 150 mM Tris-HCl, pH 8.0) and 4 μL of a solution containing reducing agents (8 M urea, 150 mM Tris-HCl, pH 8.0) were added, and the mixture was incubated at 60 °C for 1 h. Subsequently, 2 μL of cysteine-blocking reagents were added at room temperature to collect the peptide. The resulting peptide mixture was labeled using the 8-plex iTRAQ reagent, in accordance with the manufacturer’s instructions (Applied Biosystems, Inc., Foster City, CA, US). The bh50 samples were labeled 113 and 114; bh51 samples were labeled 116, 117, and 118; and bh53 samples were labeled 118, 119, and 121 (Fig. 1). All labeled samples were then multiplexed and vacuum-dried for further identification.

2D-LC separation and RPLC-MS/MS analysis

Prior to LC-MS/MS analysis, excess labeling reagent was removed by strong cation exchange chromatography. Peptides were separated using an Agilent 1200 HPLC System (Poly-SEA 5 μL 300 Å 2.0 × 150 mm with 215 nm and 280 nm UV detection; Agilent, Santa Clara, CA, USA) at 0.3 ml/min using a nonlinear binary gradient starting with buffer A and transitioning to buffer B (10 mM KH₂PO₄, pH 3.0, 500 mM KCl, 25% acetonitrile). The first fraction was collected at 0–5 min. Thereafter, fractions were collected at 4 min intervals from 6 to 44 min, and the final fraction was collected at 45–50 min to generate a total of 12 fractions. Each fraction was dried in a vacuum freeze dryer for LC-MS/MS analysis. Peptide samples were resuspended in Nano-RPLC buffer A (0.1% formic acid, 2% acetonitrile). Online Nano-RPLC was performed using an Eksigent nanoLC-Ultra 2D system (AB SCIEX, Framingham, MA, USA). Proteins were identified using two technical replicates for each biological replicate. Samples were loaded onto a C₁₈ nanoLC trap column (100 × 3 cm, C₁₈, 150 Å) and washed with buffer A at 2 μL/min for 10 min. An elution gradient of 5%–3.5% acetonitrile (0.1% formic acid) over 70 min was used on an analytical ChromXP C₁₈ column (75 μm × 15 cm, C₁₈, 3 μm 120 Å) with a spray tip. Data were acquired using a Triple TOF 5600 system (AB SCIEX) fitted with a Nanospray III source and a pulled quartz tip as the emitter. For information-dependent acquisition, survey scans were acquired over 250 ms, and as many as 35 product ion scans were collected if these exceeded a threshold of 150 counts per second with a charge of 2⁺ to 5⁺.

iTRAQ data analysis

MS/MS spectra were processed with Protein Pilot software (Protein Pilot 4.0; AB SCIEX) against the Triticum aestivum database using the Paragon algorithm [21]. An automatic decoy database search strategy was employed to estimate the false discovery rate (FDR) using the PSPEP (Proteomics System Performance Evaluation Pipeline Software, integrated in the Protein Pilot software), which was calculated as the false positive matches divided by the total matches. Proteins were identified using the following parameters: sample type = iTRAQ 8-plex (peptide-labeled), Cys; alkylation = iodoacetamide; digestion = trypsin; instrument = TripleTOF5600 (AB SCIEX); database = Triticum aestivum.fasta. The iTRAQ 8-plex was chosen for protein quantification, and unique peptides/proteins with a global FDR of < 1% were considered for further analysis.

Bioinformatic analysis of differential abundance proteins (DAPs)

To identify DAPs from the protein data, we used the following criteria: number of unique peptides ≥2; fold change (FC) threshold for upregulation/downregulation = 1.5/0.67; and maximum allowed FC = 100. The website-based tools agriGO (http://bioinfo.cau.edu.cn/agriGO/analysis.php) and KOBAS (ver. 3.0; https://www.biostars.org/p/200126/) were used to analyze the DAPs identified for Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. An online tool for generating Venn diagrams (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to visualize overlapping DAPs between groups (bh51 compared with bh50 and bh53 compared with bh50).

RT-PCR assays

Total RNA from immature spike tissue was isolated using Trizol reagent (Invitrogen, Carlsbad, CA, USA) in accordance with the manufacturer’s instructions. Approximately 2 μg of DNA-free total RNA and a PrimeScript 1st Strand cDNA synthesis kit (TaKaRa Bio, Inc., Dalian, China) were used for first-strand cDNA synthesis. Quantitative real-time PCR reactions were performed on a Mastercycler ep realplex real-time PCR system (Eppendorf, Hamburg, Germany) with SYBR Premix Ex Taq (TaKaRa Bio, Inc.). The cDNA was amplified using specific primers (Additional file 1: Table S1). Primer pairs were designed for qRT-PCR using the local AlleleID software (ver. 6.0; http://www.softpedia.com/get/Science-CAD/AlleleID.shtml). The expression levels of the corresponding genes are presented relative to corresponding control samples under the indicated conditions, with normalization against the Tubulin internal control gene (F: GACGCATGTCCATGAAGGAG; R: CCAATGCAAGAAAGCCTTGC). Relative gene expression levels were calculated using the 2^–△△CT method [22]. Correlation coefficients and histograms were calculated and plotted using R software for statistical computing (ver. R-3.0.0; http://www.r-project.org).

The developmental characteristics of RSs and FRSs during the spikelet initiation phase

To determine the developmental characteristics of spikes in the NILs, seedling morphology and spike architecture were examined. The bh51 and bh53 lines had classic SS traits, and both had more spikelets and grains per spike than NS phenotype (bh50, wild-type) plants (Fig. 1) [2]. Immature spikes from the main tiller of the stem were collected at the seedling stage for further examination by scanning electron microscopy. The results showed that lateral meristems emerged simultaneously between the glume and lemma, indicating that this could be a key stage in the transformation of the spikelet meristem to a floret meristem. Tissues from all three NIL genotypes were harvested for 8-plex iTRAQ analysis.

The principal component (PCA) and cluster analysis

An iTRAQ-based shotgun quantitation approach was used to compare overall differences in the proteomes of the three experimental groups (bh50, bh51, and bh53). A total of 6249 proteins were identified with FDR < 1%, and 3834 changes in protein abundance were assessed (Additional file 2: Table S2). As a starting point for the DAP analysis, expression data were used to determine the global relationship between the different genotypes. PCA was performed using the individual replicates from each genotype. The results demonstrated that although the first two principal components could only explain 41.58% of the variation, the data points in the graph formed tight clusters, especially the controls (Fig. 2a). Additionally, we were able to separate the replicates from each genotype into three discrete groups (Fig. 2b), which is a further indication that the data are reliable.

Expression patterns of DAPs

Applying the cut-off threshold of a 1.5-fold change for increased accumulation and a 0.67-fold change for decreased accumulation, together with the number of unique peptides ≥2, totals of 70 and 74 proteins showed differential accumulation in the bh51 versus bh50 and bh53 versus bh50 groups, respectively. Among these DAPs, 25 proteins were upregulated and 45 were downregulated in bh51 compared with wild type (bh50). In contrast, the expression of 27 proteins was increased and 47 decreased in bh53 compared with bh50 (Additional file 3: Table S3). When a Venn diagram was used to determine common DAPs among the three experimental groups, 27 proteins were unique to the bh51/bh50 group, 24 were unique to the bh53/bh50 group, and 38 were shared between the two groups (Fig. 3). Additionally, almost all the shared DAPs exhibited the same fold-change trends. For example, 97.37% of the proteins were either upregulated or downregulated when bh50 protein expression was compared with bh51 and bh53 expression, whereas only one protein (W5DNF1) changed the direction of its expression between the two experimental groups (Additional file 3: Table S3).

Functional annotation and categories of DAPs

The proteins that overlapped in the two experimental groups (38 DAPs) were analyzed using bioinformatics to understand the molecular pathways and processes involved. A total of 12 proteins were annotated in the biological process (BP) category, and 20 terms were significantly enriched. In contrast, 26 proteins were annotated in the cell component (CC) category, and 20 terms were significantly enriched. A total of 20 proteins were annotated in the molecular function (MF) category, including 10 significantly enriched terms (Fig. 4). In the BP analysis, “post-embryonic development” and “organelle organization” accounted for the largest proportion of DAPs (adjusted p-value = 4.63E^− 02), followed by biological processes that included “chromosome organization”, “peptidyl-amino acid modification”, and “histone modification” (Additional file 4: Table S4). The CC analysis demonstrated that most of the annotated proteins were intracellular. Binding activity was the dominant MF among the GO assignments.

KEGG analysis provides more information on spikelet development

Temporal and spatial expression of candidate genes

To look for correlations between mRNA levels and protein abundance, transcriptional analysis of six proteins was performed using qRT-PCR (Additional file 1: Table S1). Transcript levels of the genes and abundance of the corresponding proteins displayed similar trends for W5A3B9 (protein argonaute 4B-like), A0A077RWA8 (BEL1-like homeodomain protein 9), W5GU54 (peptidyl-prolyl cis-trans isomerase), W5ASV5 (D-3-phosphoglycerate dehydrogenase), and W5FH53 (putative DNA repair protein RAD23). The only exception was W5BLU3 (involved in histone methylation) for the bh51–bh50 group comparison. In general, W5A3B9, A0A077RWA8, and W5ASV5 were downregulated, whereas W5GU54 and W5FH53 were upregulated. Individually, the expression and abundance of W5BLU3 only increased in the bh53–bh50 group comparison, whereas there was a discrepancy between the proteomic and transcript data in the bh51–bh50 group comparison (Fig. 5). Quantitative analysis demonstrated that transcript levels of W5ASV5 measured by qRT-PCR in the bh51–bh50 and bh53–bh50 group comparisons changed 0.505- and 0.456-fold, respectively, whereas the corresponding protein levels were altered 0.476- and 0.531-fold. The linear regression equation (R² = 0.4907) and correlation coefficient (0.7005) are also shown in Fig. 5.

Putative GO and KEGG analysis networks

To examine some of the regulatory networks more closely, DAPs of particular interest were submitted to the STRING database for interactive pathway analysis based on data from Arabidopsis (Additional file 5: Table S5) [24]. The pathways/functions include regulation of transcription factor activity (BLH8), regulation of transcription (ADA2B, AGO4), glycolysis (GAPC1), removal of reactive oxygen species (APX1), heat shock protein-related (Hsp90.1, WPP2) and probable histone-arginine methyltransferase (PRMT4B). Because the number of DAPs was limited, independent spots, such as protein folding (CYP40) and even serine biosynthesis (PGDH1), are not shown in the network in Fig. 6. This information provides insight into the mechanism of lateral meristem growth, and possible interactions between these proteins could be tested by yeast two-hybrid experiments in the future.

Wfzp is a member of a gene family characterized by two AP2 domains. Other members include APETALA 2 and proteins that participate in inflorescence meristem development [25, 26]. Unlike ERF-type proteins, these AP2-type transcription factors usually bind to the promoters of genes containing the consensus motif: 5’-gCAC(A/N)N(A/T)TcCC(a/g)ANG(c/t)-3′ [27]. Additionally, all of them contain an miR172 target site, including TARGET OF EAT1 (TOE1), TOE2, TOE3, SCHLAFMUTZE (SMZ), and SCHNARCHZAPFEN (SNZ) in Arabidopsis [28, 29]. Many of these regulatory pathways have been investigated in model plants. Zhu (2011) proposed a conserved regulatory network with key roles for miR156 and miR172 in triggering changes in the developmental phase and in flower development in Arabidopsis and monocotyledons [30]. When over-expressed in plants, the 21-nucleotide non-coding microRNA miR172 was able to downregulate AP2 target genes, including IDS1 [31], SID1 [32], and SNB [33]. In this study, we also identified downregulated proteins (e.g., W5A3B9) that bound to miR172 and were involved in targeting AP2 mRNAs and mediating the cleavage of miRNA targets. Additionally, A0A077RWA8, which acts synergistically with BLH8 in Arabidopsis, is shown in the network to interact with AP1-like proteins, including LFY [34] and AGL8 [35]. The fact that our results are consistent with previous findings confirm their reliability. We further investigated some of the key proteins in the network that are associated with major biological processes to examine any links between their functions and inflorescence differentiation, although the precise nature of any interaction between these proteins has yet to be established.

Histone modification and post-transcriptional gene regulation

The W5A3B9 protein and its Arabidopsis homolog ARGONAUTE 4 (AGO4) have important roles in RNA-mediated silencing systems, including post-transcriptional gene silencing in plants [36]. Additionally, 24-nt miRNAs, including miR172, are recruited by AGO4 to induce methylation and silence target genes [30, 37, 38]. The expression level of AGO4 protein decreased in both of our group comparisons, and similar results were observed in our qRT-PCR experiments (Fig. 5). Although the expression level of miR172 in the NILs is unclear, we predict that the target AP2-like genes could be upregulated. Previous results demonstrated that growth inhibition in axillary meristems was affected by single-nucleotide polymorphism mutations in the ORF of WFZP. Therefore, mutations in FZP may stimulate ectopic development at the spikelet meristem and influence organ identity at the floral meristem [4, 6, 7].

Histone methyltransferases are histone-modifying enzymes (e.g., histone-lysine N-methyltransferases and histone-arginine N-methyltransferases) that catalyze the transfer of methyl groups to histone lysine and arginine residues [39]. Methylated histones can either repress or activate transcription [40]. W5BLU3 is orthologous to Arabidopsis PRMT4B and may be recruited to promoters upon histone H3 methylation, activating transcription via chromatin remodeling [41]. In our study, the expression of this protein was increased in both group comparisons and in the qRT-PCR results comparing bh53 with bh50, which suggests that it may activate the expression of downstream genes. Additionally, FRIGIDA-like protein 4a (FRL4A) in Arabidopsis not only regulates flower development but also participates in cell proliferation [42]. In rice, orthologs of PRMT5 appear to be key mitotic regulators and are essential for the accurate pre-mRNA splicing that regulates flowering in plants [43].

In Arabidopsis, the transcriptional adapter ADA2b (a homolog of A0A1D5YPJ3) is required to stimulate acetyltransferase activity of GCN5 on free histones or nucleosomes by opening up the promoter region [44]. It may also mediate auxin and cytokinin signaling to control cell proliferation and play a role in repressing freezing tolerance pathways at warmer temperatures [45]. Because the A0A1D5YPJ3 protein was upregulated in both our bh51 and bh53 proteomic group comparisons, this gene may promote transcriptional activity in lateral meristem development.

DAPs in the ubiquitination pathway could affect protein degradation

The ubiquitin pathway plays important roles in regulating cell cycle progression via 26S-proteasome-mediated protein degradation and involves two different E3 ligase complexes called SCF (skp-cullin-F-box protein) and APC/C (anaphase-promoting complex). Many studies suggest that this pathway is also involved in regulating spike morphogenesis [46]. Ikeda et al. (2007) reported that Aberrant Panicle Organization 1 (OsAPO1) encodes an F-box protein that regulates spikelet and floral identity and determinacy in rice and is also an ortholog of Unusual Floral Organ (UFO) in Arabidopsis [47]. Additionally, the Arabidopsis AtDA1 protein can act as a ubiquitin receptor and regulate cell proliferation, producing larger floral organs [48]. Therefore, this pathway may also include other proteins relevant for our study. The ubiquitin receptor RAD23d associates with the 26S proteasome docking subunit RPN10 to facilitate the recognition of ubiquitinated substrates for ubiquitin/26S proteasome-mediated proteolysis [49]. RPN10 plays a role in balancing cell expansion with cell proliferation rates during shoot development [50]. The expression levels of RPN10 in both the proteomic and qRT-PCR results were elevated, which may indicate that protein degradation is accelerated. Therefore, the function and network associated with the RAD23d receptor gene should be investigated more thoroughly in future studies.

Transcription factors involved in shoot development

BEL1-like homeodomain protein 8 (BLH8) is homologous to A0A077RWA8 and functions as a transcription factor maintaining stem cell fate in the shoot apical meristem, specifying floral primordia, and establishing early internode patterning events during inflorescence development [51]. It also represses the transcription of AGAMOUS (AG) in floral and inflorescence meristems [52]. Eleven BLH8 regulation network interactors were identified using the STRING database annotations (http://string-db.org/cgi/input.pl?UserId=oWlcBNoivY9Y&sessionId=8Pmofu57T0vk&input_page_show_search=on). These included AP1 (A-class MADS box gene) and a member of the KNOTTED class of homeodomain proteins (KNAT1, KNAT2, KNAT5, and KNAT6) encoded by the STM gene, which has a previously demonstrated role controlling development of the inflorescence meristem [53]. Additionally, a close correlation was observed between the proteomic studies of A0A077RWA8, which was downregulated approximately 2-fold in bh51 compared with that in bh50 (0.591) and in bh53 compared with that in bh50 (0.463). Because A0A077RWA8 represses AG (a C-class MADS box protein), it may play a role in floral meristem development by an unknown molecular mechanism.

DAPs involved in carbohydrate and nitrogen metabolism

The Arabidopsis homolog of W5ASV5, D-3-phosphoglycerate dehydrogenase 1 (PGDH1), is required for the phosphorylated pathway of serine biosynthesis and for mature pollen development [54]. Phosphorylation of histone Ser residues is important for regulating cell division and development [55]. The proteomic and mRNA results from both of our group comparisons showed an approximately 2-fold decrease in the expression of W5ASV5. Because this enzyme is involved in an energy release process in glycolysis, the protein may act as a negative regulator in lateral meristem developmental processes. We will investigate this gene further in our next study.

The Arabidopsis homolog of W5HB91, glyceraldehyde-3-phosphate dehydrogenase (GAPC1), is a key enzyme in glycolysis and plays essential roles in the maintenance of cellular ATP levels and carbohydrate metabolism [56]. It is also involved in fruit and seed development because the gapc-1 null mutant has aborted/empty embryonic sacs in its basal and apical siliques. The relationship between the RS and FRS phenotypes and decreased enzyme activity merits further study.

Upregulated HSPs may be involved in developmental processes

Peptidyl-prolyl cis-trans isomerase (CYP40), which is orthologous to W5GU54, catalyzes the cis-trans isomerization of proline imidic peptide bonds in oligopeptides and is involved in promoting the juvenile phase of vegetative development, regulating the position of floral buds, floral morphogenesis, and the expression of HSPs [57]. In our study, the proteomic and qRT-PCR expression patterns for W5GU54 were consistent, demonstrating an upregulation in both the RS and FRS genotypes. Another vital protein, W5DZ32, which is homologous to WPP2 in Arabidopsis, was developmentally associated with the nuclear envelope and promoting cell division [58]. WPP2 could act together with HSP70–1 to facilitate WIT1 targeting at the nuclear envelope [59]. The upregulated expression from the proteomic assays suggests that these genes may promote a transition from vegetative to reproductive development.