text:
Effect of Tomato Pomace Addition on Chemical, Technological, Nutritional, and Sensorial Properties of Cream Crackers

The aim of this research was to determine the influence of tomato pomace (TP) addition on the chemical, nutritional, and technological characteristics of cream crackers made from wheat flour and 4%, 6%, 8%, and 10% TP. The TP-enriched cream crackers showed progressively increasing ash (from 0.69 of the control to 1.22 g/100 g dry matter of the 10% TP sample), fat (from 11.39 to 13.04 g/100 g), protein (from 13.53 to 15.60 g/100 g), total dietary fibre (from 4.08 to 7.80), carotenoids (from 0.55 to 8.56 mg/kg), tocols (from 57.59 to 71.63 mg/kg), free phenolic acids (from 100.08 to 277.37 mg/kg), free flavonoids (from 0.0 to 45.28 mg/kg), bound flavonoids (from 0.0 to 27.71 mg/kg), and fatty acids contents, antioxidant activity and dough viscosity. The colour coordinates increased via augmenting the amounts of TP. Thickness, volume, and specific volume decreased gradually with increasing TP; the enrichment reduced cracker hardness from 65.42 N (control) to 26.28 N (crackers with 10% TP), while the snapping force rose. Cream crackers with 8% TP showed the best sensory quality. Tomato pomace addition improves the nutritional quality of foods; furthermore, its recycling will help to solve the problems linked to the disposal of this industry waste.

Introduction
By-product management during fruit and vegetable processing is one of the biggest problems for the agro-food industry, because often large quantities of leftover pomace are thrown directly in landfills, causing serious environmental problems [1]. In a circular economy, waste materials are advantageously recycled as new raw materials. Presently, pomace is used as soil fertilizer [2] or as animal food [3], but it could be better employed as raw material to produce new functional products for human nutrition.
Tomato (Lycopersicon esculentum L.) is one of the most consumed fruits/vegetables in the world, with a 2020 annual production of 187 million tons [4]. Tomatoes are eaten fresh or are processed into tomato juice, tomato paste, etc. The waste, commonly called tomato pomace (TP), contains peels, seeds, and small pulp leftovers [5], and is around 4 million tons a year [6]. The recycling of this waste may be difficult because of rapid microbial decomposition; on the other hand, recycling is environmentally important and economically beneficial because TP contains important biologically active substances. The nutritional value of TP largely depends on its composition. In fact, tomatoes contain antioxidant compounds, such as tocols and carotenoids [7], phenolic compounds, vitamins C and E [8,9], as well as minerals (potassium, phosphorus, magnesium, calcium, and sodium) [10]; the peels have 2.5 times more lycopene than seeds and pulp [11], while the seeds are rich in protein, polyphenols, fibre, and lipids (mainly linoleic, oleic, and palmitic fatty acids) [12,13]. Due to the excellent nutritional components, after drying and grinding, TP can be used as raw material for the preparation of new products with nutritionally enhanced characteristics. So far, several studies have been conducted on the effect of tomato supplementation on soda crackers [14], bread and muffins [10], extruded barley snacks [15], biscuits [16], and cakes [17], but to the best of our knowledge, no information about the influence of tomato pomace on the antioxidant fate of bakery products, in particular cream crackers, is currently available. Hence, this study aims to evaluate the effects of the addition of increasing tomato pomace quantities (0%, 4%, 6%, 8%, and 10%) on the antioxidant properties (profile of carotenoids, tocols, free phenolics and bound phenolics, and antioxidant capacity) and on the physical, chemical, technological, and sensory characteristics of cream crackers.

Cream Crackers Preparation
The cream crackers were prepared according to Nakov and Ivanova [18], with minor changes. In short, to a dough made by blending wheat flour (250 g), margarine (37.5 g), salt (3.75 g), and sugar (3.125 g) with a MMC700W mixer (Gorenje, Velenje, Slovenia) for 3 min were added 70 mL water and a yeast suspension prepared with 6.25 g fresh yeast and 30 mL water at 35 °C. The resulting dough was mixed for 2 min at high speed and then rested for 24 h at 30 °C in a fermentation chamber. The enriched crackers were made by substituting wheat flour with 4%, 6%, 8%, and 10% TP, percentages chosen after thorough research of the scientific literature. The crackers without TP added (0%) served as a control sample. The fat dough was prepared by mixing for 2 min 62.5 g flour, 0.94 g salt, 1 g baking powder, and 40 g margarine with the above-mentioned MMC700W mixer.
The leavened dough was reduced to 6 mm thickness with a dough sheeter (Professional LMP500 laboratory laminator, Electrolux, Stockholm, Sweden). A 1 mm layer of fat dough was placed over the cracker dough, and the resulting double layer was folded positioning the fat layer inside two dough layers, and then it was laminated and folded again. After resting for 30 min at 4 °C in a polyethylene bag, four further laminations were performed, interspersed with three additional dough rests. Finally, the cream crackers were formed with a manual cutter, their surfaces were pierced to facilitate quick CO2 release, and they were baked at 240 °C for 7 min in a convection oven (Wiesheu Minimat Zibo, Wiesheu GmbH, Großbottwar, Germany). Five types of crackers were prepared, i.e., control crackers (100% wheat flour), and crackers enriched with 4%, 6%, 8%, and 10% TP. Two 40-cracker sets were produced from each flour/enriched mix.

Pasting Properties
The pasting properties of the wheat flour and the TP-enriched mixtures were tested using a Micro ViscoAmylo-Graph (Brabender OHG, Duisburg, Germany). The following parameters were determined: peak viscosity, breakdown, setback and final viscosity (in Brabenber Units, BU), and pasting temperature (in °C). The measurements were performed in duplicate.

Physical Characteristics
Mean thickness and width (mm) were measured on five crackers. The colour of dough and crackers in the CIE L*, a*, b* system was scored with a Chroma Meter CR-400 colorimeter (Konica Minolta, Tokyo, Japan) with the use of illuminant D65 on two sets of five random crackers. All data were presented as mean values of at least five replicates from each batch. The volume (cm 3 ) was determined on three samples using a Volscan Profiler (Stable Micro Systems Ltd., Godalming, UK). The spread factor was calculated as width/thickness multiplied by 10. The specific volume was calculated as volume divided by weight of the same cracker. The textural properties of the crackers (hardness and snapping force) were measured using a TA-XT2 Plus texture analyser (Stable Micro System Ltd. Godalming, UK). Images of cream crackers with different quantities of tomato pomace are shown in Figure 1.

Chemical Characteristics
Moisture and ash contents of wheat flour, TP, and crackers were determined according to methods ISO 665:2020 [19] and ISO 5984:2002 [20], respectively; protein content was assessed by the modified Lowry method, as described in Mae hre et al. [21]; lipid concentration was tested according to Soxhlet method ISO 136 [22]; and fibre content (soluble, insoluble, and total) was assessed as in method 32-07.01 [23]. Total carbohydrate content was computed by subtracting protein, lipid, and ash from the total. The results are reported as g/100 g dry matter (DM).
Identification of fatty acids in tomato pomace was performed using a Shimadzu-2010 gas chromatograph (Kyoto, Japan). The assay was performed with a CP7420 capillary column (100 × 0.25 mm, i.d. 0.2 μm, Varian Inc., Palo Alto, CA, USA) with carrier gas (hydrogen) and make-up gas (nitrogen). A gas chromatographic oven with five-step temperature program was used. All the analyses were performed in duplicate.

Tocols and Carotenoids
The analyses of tocols, trans-lutein, and trans-zeaxanthin were scored by normal phase HPLC, while lycopene and carotenoid isomers were analysed by reverse phase HPLC.
The tocol and carotenoid extracts were obtained from all the samples after saponification [24]. The NP-HPLC analysis was performed as described by Hidalgo and Brandolini [25]. The operating conditions for tocol identification were: Adamas ® Silica column 250 × 4.6 mm, 5 μm, and guard cartridge 10 × 4.6 mm, 5 μm (Sepachrom SRL, Rho, Italy); mobile phase, hexane:ethyl acetate:acetic acid (97.3:1.8:0.9, v/v/v); flow rate, 1.6 mL/min; pump L-2130 Elite LaChrom (VWR, Hitachi, Japan); fluorimetric detector Jasco 821 FP Intelligent Spectrofluorometer (Tokyo, Japan) at excitation-emission wavelengths of 290 nm and 330 nm, respectively; connected to a computer with the software Empower 2 (Waters Chromatography Division, Millipore, Milford) through the Waters e-SAT/IN module. Peaks were quantified using α-tocopherol, β-tocopherol, γ-tocopherol, and δ-tocopherol as external standards. The tocotrienols were quantified using the standard curves of their corresponding tocopherol. The operating conditions for lutein and zeaxanthin detection [26] were: Adamas ® Silica column 250 × 4.6 mm, 5 μm, and guard cartridge 10 × 4.6 mm, 5 μm (Sepachrom SRL, Rho, Italy); column oven at 20 °C L-2300 Elite LaChrom (VWR, Hitachi, Japan); mobile phase, hexane:isopropyl alcohol (5%); flow rate, 1.5 mL/min; pump L-2130 Elite LaChrom (Hitachi, Tokyo, Japan). Carotenoids were detected at 445 nm by a PhotoDiode Array Detector (Waters, Milford, MA, USA) set in the range of 200-650 nm. The HPLC system was controlled by the software Empower 2. The lycopene and carotenoid isomer quantification was carried out by reverse phase HPLC; the samples were prepared as described by Hidalgo et al. [27]. The carotenoid chromatographic analysis was carried out following Gupta et al. [28]. Exactly 20 μL volume of filtrate was tested under the following operating conditions: L-2130 Elite LaChrom pump (VWR, Hitachi, Japan), L-2300 Elite LaChrom column oven (VWR, Hitachi, Japan), diode Array Detector L2450 Elite LaChrom (Merck, Hitachi, Japan), at wavelength 445 nm, with EZChrom Client/Server software version 3.1.7. The analysis was carried out using an Adamas C30 column 5 μm-250 mm × 4.6 mm-HPLC Column (Sepachrom Srl, Rho, Italy), an Adamas pre-column thermostated at 20 °C, using mobile phases consisting of (A) methanol/water (98: 2, v/v), (B) methanol/water (95: 5, v/v), and (C) tert-methyl butyl ether. The gradient elution was 80% A, 20% C at 0 min, followed by a linear gradient at 60% A, 40% C at 2.0 min at a flow rate of 1.4 mL/min, at 2.01 min flow rate was changed to 1.00 mL/min with gradient changing to 60% B, 40% C followed by a linear gradient at 0% B, 100% C within 12 min and returning to initial conditions within 13 min. A rebalancing (7 min) was performed at initial concentrations of 80% A, 20% C. The column temperature was maintained at 20 °C. Elution peaks were monitored in a range of 250 to 700 nm using PDA. The quantification of the peaks was performed with the external standard method by preparing calibration curves of lycopene and β-carotene for all the other compounds. The analyses were carried out in duplicate and the results expressed in mg/kg dry matter (DM).

Free and Bound Phenolic Compounds
The extraction of the soluble free and insoluble bound phenolic compounds was carried out according to the procedure reported by Brandolini et al. [29]. Approximately 1 g of flour was weighed into capped centrifuge tubes, and 15 mL of a methanol solution (80%) was added to each tube. After vortexing with a Mn1Minishaker (Ika, Staufen, Germany) for about 45-60 s, the samples were placed in an ultrasonic bath for 10 min, vortexed again for a few seconds, and centrifuged at 12,000 RPM for 10 min at 8-9 °C with a Centrikon K24 centrifuge (Kontron Instruments, Bletchley, UK). The supernatants, which contain the free phenolics, were recovered in a 250 mL flask, covered with foil to protect the content from light. The extraction was repeated two more times. The three extracts were merged, evaporated under vacuum using a Laborota 4000 rotavapor (Heidolph Instruments GmbH & Co. KG, Schwabach, Germany) for 40 min at 35 °C and completely dried by nitrogen flow for 1 min. The dry samples were resuspended with 2 mL of methanol:chromatographic water (8:2 v/v), vortexed, and filtered on a 0.45 μm PTFE membrane (Diana Beck Scientific, Angera, Italy).
The solid residues, which contained the insoluble bound phenols, were digested with 15 mL of 4M NaOH under nitrogen for 4 h at room temperature and constant agitation, brought to a pH of 1.5-2 with 6M HCl, and extracted twice with 20 mL of diethyl ether/ethyl acetate (1:1, v/v). After centrifugation (11,178 g, 10 min, 8 °C), the supernatants were clarified with sodium sulphate, filtered, evaporated as previously outlined, resuspended in 2 mL of methanol:water (1:1 v/v), and filtered with a 0.22 μm PTFE membrane (Millipore, Carrigtwohill County, Cork, Ireland). All extractions were performed under dim light, and the extraction tubes were wrapped with black paper to avoid sample degradation by photooxidation.
A 20 µ L volume of filtrate underwent reverse phase HPLC analysis using a column Adamas ® C18-AQ 5 μm 4.60 mm × 250 mm and a precolumn C18 5 μm 4.60 mm × 10.0 mm (Sepachrom SRL, Rho, Italy) thermostated at 30.0 °C; L-2130 pump, L-2300 column oven and L2450 Diode Array Detector Elite LaChrom (Hitachi, Tokyo, Japan). The identity of the compound was confirmed by the congruence of retention times and UV/VIS spectra with those of pure authentic standards. Thirty-three standards were injected; unidentified peaks were quantified using the calibration curve of the compound with a similar absorption spectrum and named as "phenolic derivative" as indicated in Brandolini et al. [29] and according to Dueñas et al. [30], Zalewski et al. [31], and Estivi et al. [32]; additionally, the derivatives of each compound were grouped together. The following calibration curves were constructed using 4-6 concentrations: apigenin 1-20 mg/L (detection limit

Antioxidant Capacity
The samples were extracted as suggested by Varas Condori et al. [33] with hexane (lipophilic fraction) and 80% methanol (hydrophilic fraction) in succession. Exactly 0.100 ± 0.010 g of sample was weighed in Eppendorf tubes and resuspended with 1 mL hexane solution. The suspension was vortexed with a WIZARD Advanced IR Vortex Mixer (Velp Scientifica Srl, Italy) and sonicated (Bandelin Sonorex, Berlin, Germany) for 10 min at 4 °C in the dark; subsequently, the samples were centrifuged (Centrifuge mod. 4224, ALC, Winchester, MA, USA) at 12,000 rpm for 10 min. The supernatant was thus recovered, and the extraction was repeated; the second supernatant was then combined with the previous one. For the hydrophilic fraction, 80% methanol was used on the sediment with a procedure similar to the previous one.
The antioxidant capacity of the lipophilic extracts was tested by the ABTS method [34]. To 150 µ L of extract, 3.0 mL of a radical solution consisting of 10 mL of aqueous 7 mM ABTS solution and 170 µ L of 140 mM aqueous potassium persulfate solution (absorbance value of 0.70 ± 0.02) were added and subjected to 6 min of reaction at 30 °C in a thermostated bath (EN.CO. Srl, Spinea, Italy). The absorbance of the extracts was measured at 734 nm using plastic cuvettes with a spectrophotometer (JASCO V-650 Spectrophotometer, Tokyo, Japan).
The antioxidant capacity of the hydrophilic fraction was measured by the ABTS method [34], as described above, and also by the FRAP method [35], slightly modified. Exactly 3.0 mL of the FRAP reagent were added to 150 μL of the extracted sample and, after 60 min of reaction at 37 ° C in a thermostated bath (Thermo Electron Corporation, EN.CO. Srl, Italy), the absorbance was measured at 593 nm using cuvettes with a spectrophotometer (JASCO V-650 Spectrophotometer, Tokyo, Japan). Both analyses were carried out in duplicate, and the antioxidant capacity was expressed in mmol of Trolox equivalent (TE)/kg DM.

Sensory Analysis
The sensory analysis of the five types of cream crackers with different TP contents was performed at the sub-department of Cereal Technology at the Faculty of Food Technology in Osijek, Croatia. Twenty trained people, after giving informed consent according to the guidelines on Ethics and Food Related research defined by the European Union [36], participated in the sensory tests. Appearance, texture, odour, aroma, and taste were scored from 1 to 5, where 1 was extreme dislike and 5 was extreme like; the overall quality was computed as averages of the five traits evaluated.

Statistical Data Analysis
The experimental data were analysed by one-way analysis of variance (ANOVA). When significant differences (p  0.05) were found, Fisher's least significant difference (LSD) at p ≤ 0.05 was computed. The analyses were performed using the Centurion XVI statistical program (Statgraphics Technologies Inc., The Plains, VA, USA). The mean and standard deviation values were calculated using the Office Excel program (Microsoft Corporation, Redmond, WA, USA). Table 1 reports the results of peak viscosity, breakdown, setback, final viscosity, and pasting temperature of the wheat flour and of the TP-enriched mixtures. The ANOVA showed that substituting wheat flour with TP significantly influenced (p ≤ 0.05) all these parameters. The progressive increase in final viscosity with the addition of TP may be related to the higher content of dietary fibre (cellulose, hemicellulose, and pectin), which acts as a hydrocolloid [37] and has a higher gelatinization capability. Due to the content of insoluble dietary fibre and protein (mostly coming from the tomato skin), TP can be a component for altering water sorption and regulating the rheological properties of food [38]. Nakov et al. [39] noticed a similar trend in the dough of cakes enriched with grape seed powder and attributed it to the high content of soluble dietary fibre, and/or to the high lipid content of TP, which may interact with other hydrophobic substances (e.g., gluten), thus increasing viscosity. On the other hand, Mironeasa and Gabriela [40] noticed a decrease in final viscosity, attributed to the increase in tomato pomace, a non-starch component, and to the concomitant reduction of wheat flour. Similarly, Sogi et al. [41] analysed the effect of tomato seed meal on the rheological characteristics of dough and noticed that replacing wheat flour with 10%, 20%, and 30% tomato seed meal did not change pasting temperature, but peak viscosity, breakdown, setback, and final viscosity decreased. Table 2 shows the colour of the dough and of the cream crackers, as well as the physical and textural characteristics of the cream crackers produced without (control) or with TP addition. In general, all these parameters varied according to the TP content of the cream crackers. The colour of dough and crackers became darker (lower L*) by increasing the TP amount, while a* and b* increased; that is, they became redder and yellower. The same results, due to the rich pigment content of TP, were reached by other authors [42,10].

Physical Characteristics
No significant width differences were detected, whereas the thickness and volume of the 8% and 10% TP crackers were significantly lower (р < 0.05), and the spread factor was higher than the control. It has to be stressed, however, that the incorporation of even the lowest TP amount (4%) caused a significant reduction (p < 0.05) in volume and specific volume. Bhat and Ahsan [42], Chung [43], and Bhat et al. [44] observed that adding more than 5% TP in biscuits significantly changed the diameter and thickness, probably because the higher TP fibre content diminished gas retention. A volume reduction was expected, because the wheat flour (rich in gluten, responsible for leavened products' volume) was replaced by the gluten-free TP. Additionally, bakery products' volume also depends on the ability to produce gases, which are linked to many factors such as growth agents, starch, fermentable sugars, pH, etc., but not to dietary fibre [45]. TP is rich in fibre ( Figure  2); therefore, this factor can lead to reduced cream cracker volumes. Nakov et al. [39] noticed that adding grape pomace to cakes led to width, thickness, and volume reductions proportional to the grape pomace content. TP addition significantly reduced cracker hardness from 65.42 N (control) to 26.28 N (crackers with 10% TP), and increased snapping force (the force required to break the cracker; Jukić et al. [46]) according to TP content. Chung [43] found a decrease in hardness from 60.35 kg/cm 2 (control cookies) to 38.98 kg/cm 2 (cookies with 10% TP), while Mehta et al. [10] had similar results in TP-enriched bread and muffins compared to the control samples. Probably, the hydrocolloid nature of TP inhibited the connection between starch chains and reduced starch retrogradation [37].

Chemical Characteristics
The ANOVA (not shown) proved the existence of significant differences (p < 0.05) between the control and the crackers with different amounts of TP for all parameters. Figure 1 shows the main chemical characteristics of white wheat flour and tomato pomace, as well as of crackers with and without TP.
The wheat flour moisture was far higher than that of the tomato pomace, as already evidenced by other authors [42,37]. Hence, its reduction with increasing TP percentage was observed, going from 9.46 g/100 g (control) to 5.64 g/100 g (crackers with 10% TP). Moisture reduction in biscuits with different amounts of TP was reported by Chung [43] and Mironeasa and Gabriela [40]. The amount of minerals (ash) in TP was more than 10 times higher than in wheat flour (5.46 vs. 0.51 g/100 g), and in the crackers increased from 0.69 g/100 g (control) to a maximum of 1.22 g/100 g (10% TP). A high ash amount in bakery products (soda crackers, biscuits, and cakes) with different TP enrichment was also described by Alazb et al. [17], Isik and Topkaya [14], and Sogi et al. [41].
The lipid content in wheat flour was around 15 times lower compared to tomato pomace, where the lipids come from the seeds. Mironeasa and Gabriela [40] found 19.15 g/100 g lipids in flour made only from tomato seeds, i.e., concentrations slightly higher than those of our TP (that contained both peels and seeds). The cream crackers were prepared with margarine, an important source of lipids, and as a consequence, the amount of lipids in the control crackers was 11.39 g/100 g; nevertheless, in the enriched samples they increased significantly and progressively, reaching 13.04 in the 10% TP crackers. Accordingly, Chung [43] and Bhat and Ahsan [42] found a lipid increase in biscuits due to TP presence in their composition.
The protein content of tomato pomace was 20.73 g/100 g and in wheat flour 11.83 g/100 g; the high concentration in TP was due to the presence of tomato seeds, a good source of protein. Consequently, the protein level of control crackers differed statistically from that of the TP-enriched crackers. Increasing the content of protein by adding TP has been recorded in bread and muffins [10], biscuits [42], and soda crackers [14].
Tomato pomace is considered an excellent source of dietary fibre, and particularly of insoluble fibre [47]. In fact, the high dietary fibre content of TP is a leading reason for its use in enriching bakery products [15]. In our experiments, the soluble, insoluble, and total dietary fibre recorded in TP were 4.62, 32.00, and 36.62 g/100 g DM, respectively, i.e., several times more than in wheat flour (1.32, 2.51, and 3.83 g/100 g DM). The ANOVA (not presented) and the results in Figure 1 show that the amount of soluble, insoluble, and total dietary fibre in the enriched crackers was statistically different and superior to that of the control crackers. Other authors reached similar conclusions [17,14,10].
Starch is the main component of wheat flour and is far more abundant than in TP (86.20 vs. 57.42 g/100 g). Hence, in the TP-enriched crackers, a progressive decrease of carbohydrates at increasing tomato pomace levels was scored, as already noticed by Alazb et al. [17].
As mentioned above, tomato seeds (a component of tomato pomace) are rich in fatty acids [47], especially monounsaturated (MUFA) and polyunsaturated (PUFA) [11]. The fatty acids observed in our crackers are shown in Table 3.
The ANOVA (not presented) highlighted significant differences (p ≤ 0.05) among samples for all fatty acids. The palmitic, stearic, oleic, and linoleic acids in the control were 20.28, 4.16, 31.09, and 38.60 g/100 g oil, respectively, and increased progressively according to the TP enrichment of crackers, up to 21.35, 4.92, 33.10, and 41.74 g/100 g oil, respectively, in crackers with 10% TP. When grouping the fatty acids by their saturation level (Table 3), it was evident that all the crackers contained high amounts of MUFA (again, the role of the margarine should be remembered), and that the content increase of different groups of fatty acids was directly proportional to TP richness.

Carotenoids and Tocols
The HPLC analyses recorded the presence of lutein (1.26 mg/kg DM) and zeaxanthin (0.06 mg/kg DM) in the wheat flour, and of lutein (4.14 mg/kg DM), zeaxanthin (0.36 mg/kg DM), 15-cis-neurosporene (20.59 mg/kg DM), δ-carotene-isomer (10.00 mg/kg DM), and trans-lycopene (11.89 mg/kg DM) in the tomato pomace ( Table 4). The total carotenoid content of TP ( Figure 2) was >35 times greater than that of wheat flour (1.32 vs. 46.89 mg/kg DM). The significant differences among samples for total carotenoids detected by the ANOVA (not shown) were clearly reflected in the increasing concentrations found in the TP-enriched cream crackers (Figure 2), ranging from 6.83 mg/kg DM (4% TP) to 8.56 mg/kg (10% TP). The replacement of wheat flour with tomato pomace significantly increased the lycopene content of bread and muffin [10]. Table 4. Carotenoids and tocols contents (mg/kg DM; mean ± standard deviation) of wheat flour (WF), tomato pomace (TP), and cream crackers with increasing TP percentages. Different letters indicate significant differences (p < 0.05) among cream crackers along a row. A more varied situation was observed for the tocols. In fact, the refined wheat flour contained only a limited quantity of β-tocotrienol (13.00 mg/kg DM), while the TP was rich in α-β-, γ-, and δ-tocopherols (55.30, 3.24, 116.65, and 3.47 mg/kg, respectively). Hence, the α-tocotrienols and δ-tocopherol found in the control crackers, as well as relevant quantities of the other tocopherols and tocotrienols, were due to the presence of margarine among the ingredients. Accordingly, the ANOVA demonstrated the existence of significant differences for total tocols contents between the control (44.38 mg/kg DM) and the four TP-enriched cream crackers (65.32-70.44 mg/kg DM), but not among the TP-enriched samples.

Phenolic Compounds
The refined wheat flour contained no free phenolics and only three insoluble bound phenolic acids (coumaric, ferulic, and ferulic derivative), while the TB contained several free and bound phenolic acids and flavonoids ( Table 5). As for tocols, the phenolic composition of the cream crackers was influenced by the presence of the margarine. For example, a free phenolic acid retrieved in the control and enriched cream biscuits (4-hydroxybenzoic acid derivative), not found in either the wheat flour or the TP, may stem from the margarine. In fact, the main ingredient of the margarine we used was oil from palm oil, whose fruits are very rich in 4-hydroxybenzoic acid [48], and the high temperatures during oil processing may induce a transformation towards a less polar compound. Additionally, the high and increasing concentrations of free gallic acid in the TP-enriched samples may be due to its release from different TP compounds during the long yeast fermentation step [49,50]. A less likely explanation could be the heat effect during baking, because in a yeast-less cake enriched with grape pomace, the baking step did not elicit an increase in gallic acid [39]. 98.14 ± 6.23 5.77 c ± 0.57 7.46 c ± 0.59 10.28 b ± 1.07 13.24 a ± 0.59 While the bread wheat flour contained only bound phenolic acids, the TP showed high concentrations of both free and bound phenolics ( Figure 3); furthermore, TP was very rich in flavonoids. Overall, the tomato pomace polyphenol content was more than 10 times (1211.4 vs. 73.2 mg/kg DM) that of the refined wheat flour; additionally, most of its compounds were of the free form (828.34 vs. 383.08 mg/kg DM), which has a better nutritional significance for superior bioaccessibility and bioavailability [26,51]; hence, adding TP, rich in free phenolics, will significantly improve the polyphenol composition of the cream crackers. Alazb et al. [17], using a different extraction protocol with alcohol, recovered 653.17 μg/g DM of phenolic compounds, mainly constituted (45%) by kaempferol.  The increasing TP concentrations in the enriched cream crackers led to a progressive augmentation of several soluble free phenolic acids (p-coumaric, p-coumaric der, ferulic, gallic) and flavonoids (naringenin, naringenin der, quercitin der, rutin), as well as of some insoluble bound phenolic acid (caffeic) and flavonoids (naringenin, naringenin der, quercitin der). Overall, the addition of 10% TP increased the free phenolics and bound phenolics contents by 322% and 133%, respectively, in comparison to the control, while the freeto-bound phenolics ratio jumped from 1.36 (control) to 3.28 (10% TP), significantly improving the nutritional characteristics of the cream crackers.
Isik and Topkaya [14] partially replaced wheat flour with 4%, 8%, and 12% dried tomato pomace meal in the preparation of crackers and noticed a concomitant total polyphenol content (TPC) increase from 52.52 to 127.58 mg GAE/100 g. Yagci et al. [52] concluded that increasing the levels of tomato pomace powder in extruded snacks significantly increased the content of individual free phenolics, including gallic acid, protocatechuic acid, 2,5-dihydroxybenzoic acid, chlorogenic acid, rutin, and quercetin.

Antioxidant Capacity
The antioxidant capacities of the lipophilic and hydrophilic extracts are depicted in Figure 4. The lipophilic extracts (carotenoids and tocols) showed limited antioxidant capacity in comparison to the hydrophilic extracts, possibly because of the much higher concentrations of phenolic compounds. The TP always displayed far higher values than the refined wheat flour. Accordingly, the TP-enriched cream cookies exhibited significantly increased antioxidant capacities when their TP content was augmented.

Sensory Analysis
The sensory analysis showed that the cream crackers with TP had good acceptability, insofar that those enriched with 10% TP displayed the best score for appearance, while those with 8% TP were the best for texture, aroma, taste, and odour (Figure 5a), as well as for overall acceptability, with a score of 4.4/5.0 points (Figure 5b). Mehta et al. [10] observed that adding TP in bakery products (bread and muffins) made their colour more appealing. The optimal 8% TP percentage for most of the parameters as well as for overall quality is consistent with Isik and Topkaya [14], who concluded that the addition of 12% or more TP in the production of soda crackers is not recommended, and with Bhat and Ahsan [42], who observed that biscuits with 5% TP received the best assessment. Similarly, Alazb et al. [17] did not record significant differences between a control cake and samples enriched with 2.5, 5, and 7.5% tomato pomace. The level of TP addition may be lower in high-rising leavened products such as bread. Majzoobi et al. [37] reported that bread with 5% and 7% TP was worse than the control sample, while Nour et al. [45] reported that bread with 6% TP was better than bread with 10% TP.

Conclusions
Tomato pomace is a by-product with excellent nutritional composition (high ash, fat, protein, and dietary fibre contents). Cream crackers enriched with different TP amounts (4%, 6%, 8%, and 10%) showed better physical and nutritional characteristics than the control, and the improvements were directly correlated to the TP content. The sensory analysis showed that crackers with 8% TP had the best acceptance among panellists. These results can help to promote the use of different nutrient-rich by-products from the food industry in the production of new functional foods.  Overall quality Tomato pomace content


provenance: AgriculturalAndFoodSciences,2022-2022-Agricultural And Food Sciences-2022.jsonl.gz.gz

text:
Comparative phylogenomic analyses and co-expression gene network reveal insights in flowering time and aborted meiosis in woody bamboo, Bambusa oldhamii ‘Xia Zao’ ZSX

Woody bamboos have peculiar flowering characteristics with intervals ranging from several years to more than 100 years. Elucidating flowering time and reproductive development in bamboo could be beneficial for both humans and wildlife. To identity the mechanisms responsible for flowering time and embryo abortion in Bambusa oldhamii ‘Xia Zao’ ZSX, a transcriptome sequencing project was initiated to characterize the genes involved in developing flowers in this bamboo species. Morphological studies showed that pollen abortion in this bamboo species was mainly caused by a delay in tapetum degradation and abnormal meiotic process. Differential expression (DE) and optimized hierarchical clustering analyses identified three of nine gene expression clusters with decreasing expression at the meiosis of flowering stages. Together with enriched Gene Ontology Biological Process terms for meiosis, this suggests that their expression pattern may be associated with aborted meiosis in B. oldhamii ‘Xia Zao’. Moreover, our large-scale phylogenomic analyses comparing meiosis-related transcripts of B. oldhamii ‘Xia Zao’ with well annotated genes in 22 representative angiosperms and sequence evolution analyses reveal two core meiotic genes NO EXINE FORMATION 1 (NFE1) and PMS1 with nonsense mutations in their coding regions, likely providing another line of evidence supporting embryo abortion in B. oldhamii ‘Xia Zao’. Similar analyses, however, reveal conserved sequence evolution in flowering pathways such as LEAFY (LFY) and FLOWERING LOCUS T (FT). Seventeen orthogroups associated with flowering were identified by DE analyses between nonflowering and flowering culm buds. Six regulators found primarily in several connected network nodes of the photoperiod pathway were confirmed by mapping to the flowering time network in rice, such as Heading date (Hd3a) and Rice FT-like 1 (RFT1) which integrate upstream signaling into the downstream effectors. This suggests the existence of an intact photoperiod pathway is likely the key regulators that switch on/off flowering in B. oldhamii ‘Xia Zao’.

Woody bamboos have peculiar flowering characteristics with intervals ranging from several years to more than 100 years. Elucidating flowering time and reproductive development in bamboo could be beneficial for both humans and wildlife. To identity the mechanisms responsible for flowering time and embryo abortion in Bambusa oldhamii 'Xia Zao' ZSX, a transcriptome sequencing project was initiated to characterize the genes involved in developing flowers in this bamboo species. Morphological studies showed that pollen abortion in this bamboo species was mainly caused by a delay in tapetum degradation and abnormal meiotic process. Differential expression (DE) and optimized hierarchical clustering analyses identified three of nine gene expression clusters with decreasing expression at the meiosis of flowering stages. Together with enriched Gene Ontology Biological Process terms for meiosis, this suggests that their expression pattern may be associated with aborted meiosis in B. oldhamii 'Xia Zao'. Moreover, our large-scale phylogenomic analyses comparing meiosis-related transcripts of B. oldhamii 'Xia Zao' with well annotated genes in 22 representative angiosperms and sequence evolution analyses reveal two core meiotic genes NO EXINE FORMATION 1 (NFE1) and PMS1 with nonsense mutations in their coding regions, likely providing another line of evidence supporting embryo abortion in B. oldhamii 'Xia Zao'. Similar analyses, however, reveal conserved sequence evolution in flowering pathways such as LEAFY (LFY) and FLOWERING LOCUS T (FT). Seventeen orthogroups associated with flowering were identified by DE analyses between nonflowering and flowering culm buds. Six regulators found primarily in several connected network nodes of the photoperiod pathway were confirmed by mapping to the flowering time network in rice, such as Heading date (Hd3a) and Rice FT-like 1 (RFT1) which integrate upstream signaling into the downstream effectors. This suggests the existence of an

Introduction
The transition from vegetative growth to reproductive growth is an important developmental phase in the life cycle of plants. Floral initiation is commonly regulated by various environmental and physiological cues. Bamboos (Bambusoideae) are a lineage of perennial forest grasses native to all continents except Europe and Antarctica. Bambusoideae are classified into 116 genera and approximately 1439 species (Bamboo Phylogeny Group, 2012;Clark et al., 2015;Lin et al., 2018). Bamboos consist of woody and herbaceous species. Woody bamboo species are characterized with lignified culms, complex branching, bisexual flowers, and a prolonged vegetative phase lasting from several to more than a hundred years before flowering (McClure, 1993). To biologists, the most captivating phenomenon is the extended synchronous flowering of woody bamboos which is related to the long blooming cycle (Janzen, 1976). Most woody bamboos blossom only once in their whole lifetime and usually die after flowering. Bamboo flowering is unpredictable, long-periodic, gregarious, and uncontrollable, making it unwieldy and devastating to bamboo forest resource and farmers. Once bamboos flower, the subsequent die-off can result in sudden ecological consequences-ecological environment of the bamboo forests is changed, and the economic interests of bamboo farmers and the food and living environment of pandas are no longer guaranteed. However, until now the underlying mechanism of bamboo flowering remains unclear. Therefore, it is important to ascertain the unique flowering pathways and genes involved in bamboo flowering time and floral development.
Flowering depends on many factors, including environment, nutrition, climate, and physiological status. Several hypotheses have been proposed for the flowering mechanisms in bamboo (Gielis et al., 1999;Franklin, 2004). Some MADS-box genes have been found to affect the flowering time in angiosperms (Wysocki et al., 2016). For example, SOC1 gene in Arabidopsis and its homolog OsMADS56 in Oryza sativa (rice) induce floral development (Ryu et al., 2009;Lee and Lee, 2010). And as antagonists of SOC1 gene, SVP in Arabidopsis and OsMADS22/ OsMADS55 in rice control flowering (Hartmann et al., 2000;Lee et al., 2012). In certain bamboos species, several important flowering-related genes have been identified (Tian et al., 2005;Tian et al., 2006;Lin et al., 2010;Xu et al., 2010;Hisamoto and Kobayashi, 2013), and environmental factors and different hormones induce flowering in vitro (Nadgauda et al., 1990). Moreover, the gene expression studies provide insights into the mechanisms for bamboo floral initiation and development. In 2013, the first draft of nuclear genome for Phyllostachys edulis was published, which provides a reference genome assembly for comparisons with other bamboo species . Next generation sequencing with RNA-Seq method (Wang et al., 2009) has allowed gene expression to be examined at a large scale for non-model organisms. This method was first used by Zhang et al. (2012) to study the floral transcriptome in Dendrocalamus latiflorus (a bamboo species), and then a following study on different tissues such as leaf, stem, shoot, and root of the same species by Liu et al. (2012).  presented transcriptome results from six different heights of moso bamboo fast growing shoots using the Illumina HiSeqTM 2000 platform. An active Dof-Hd3a-MADS-flowering pathway in moso bamboo was conducted by Gao et al. (2014). In Gao et al. (2015), suggested that a large number of differentially expressed microRNAs and their targets participated in diverse primary biological pathways and played significant regulatory roles in flowering of moso bamboo. In Ge et al. (2017), proposed that PheDof1, PheMADS14 and six miRNAs may play vital regulatory roles in flower development and floral transition in moso bamboo.
In Hou et al. (2021), identified a SOC1-like gene, BoMADS50, from Bambusa oldhamii, whose conserved and specific SNP function may be responsible for the existence of synchronous and sporadic flowering in bamboo. Nevertheless, until now few studies have reported the anatomy of bamboo flowers and embryo development due to so long phase of vegetative growth. Meanwhile, most bamboo species have relatively low maturing rate (Zhou, 1998;Lin et al., 2018), thus only a few studies have explored the mechanism of low seed production in bamboos (Hu et al., 1994;Huang et al., 1999;Wang et al., 2006;Lin and Ding, 2012;Lin and Ding, 2013;Lin et al., 2015;Wang et al., 2015). B. oldhamii 'Xia Zao' is an edible shoot and woody species, with great ecological and economic values. It is preferentially distributed in subtropical regions such as Fujian province, China. Current research of B. oldhamii are mainly focused on plant physiology, morphology, ecology and material science. However, there have been few molecular studies on flowering of this species (Lin et al., 2010). In 2013, B. oldhamii 'Xia Zao' flowered in the Xiapu bamboo garden, China, and there are no reports that the seeds or fruits have ever been found after flowering.
In this present study, in order to investigate the mechanism of flowering time and pollen abortion in B. oldhamii 'Xia Zao', we performed a comprehensive transcriptome sequencing for gene expression profiling and explored the cytological meiotic process with paraffin section technique. Here, we collected 27 flower and nonflower samples from B. oldhamii 'Xia Zao' representing different developmental phases: culm bud, underground shoot, above-ground shoot, and spikelet with different development phases. Using a large-scale transcriptome profiling, we interrogated the expressed transcriptomes including flower tissues at different stages. These transcriptomes were analyzed for two complementary objectives: (1) identify the key regulators that switch on/off flowering in B. oldhamii 'Xia Zao' through phylogenomic and comparative sequence analyses using transcriptomic data; (2) identify the key causal floral-related genes responsible for meiotic failure by examining the abnormal microstructure. By using a de novo transcriptome assembly, we generated a rich resource of sequence information including a list of flower and meiotic genes. Our cytological studies revealed abnormalities associated with meiotic failure, which is likely linked to nonfunctional sequence evolution of two key genes PMS1 and NEF1 with known roles in pollen and tapetum development (Ariizumi et al., 2004;. Moreover, we identified the retention of conserved photoperiod pathway in B. oldhamii 'Xia Zao' that may contribute to bamboo flowering.

Tissues, libraries, and sequence data
The samples of B. oldhamii 'Xia Zao' were collected on November 26, 2013 from the Xiapu breeding base, Ningde city, Fujian province (119°59′14″E, 26°48′10″N), China. A total number of 27 tissues were collected from three clumps: two flowering clumps (clump 1 and clump 2) and one nonflowering clump (clump 3). As B. oldhamii 'Xia Zao' has a lifestyle dominated by clonal reproduction, these three clumps were developed from the same asexual clone and were genetically identical. Six tissues including culm buds (here defined as F1_Cb, Figure 1), underground shoot buds (here defined as F1_ShU, Figure 1), above-ground shoot buds (here defined as F1_ShA, Figure 1), stage-1 flowers (spikelet length < 1cm) (here defined as F1_Fl1, Figure 1), stage-2 flowers (spikelet length between 1 and 1.5 cm) (here defined as F1_Fl2, Figure 1), and stage-3 flowers (spikelet length > 1.5 cm) (here defined as F1_Fl3, Figure 1) were collected from flowering clump 1. The buds of nonflowering culm base from clump 3 were collected as a control (here defined as NF3_Cb, Figure 1). Three biological replicates were included for above each tissue. In addition, to get an idea of the variation within different clumps, three flowering tissues of clump 2 with 2 biological replicates from each were generated (here defined as F2_Fl1, F2_Fl2, F2_Fl3). An illustration of the samples is shown in Figure 1. The samples were immediately frozen in liquid nitrogen and stored at -80°C. RNA extraction, library preparation, and Illumina paired-end sequencing were performed to generate RNA-Seq data with the read length of 2×100 bp and an insert size from 200 bp to 350 bp on the Hiseq 2000 platform. Library preparation and Illumina paired-end sequence followed the procedure described in Illumina Truseq ™ RNA sample prep kit.

Assembly, cleaning, and annotation
A total of more than 1.5 billion Illumina paired-end RNA-Seq reads were generated from 27 bamboo cDNA libraries (Table 1). To decrease the amount of data for further processing, duplicate reads were removed with Prinseq (Schmieder and Edwards, 2011). Adapter sequences were detected using fastQC and Trimmomatic (Bolger et al., 2014) was used to remove adapter sequences and low-quality bases. A de novo transcriptome assembly was subsequently performed with Trinity release (Grabherr et al., 2011) using cleaned reads from all the 27 tissues. The resulting assembly was further cleaned by removing contigs without coding regions (Iseli et al., 1999) and redundant transcripts using the Assembly Post Processing pipeline in PlantTribes (https://github.com/ dePamphilis/PlantTribes). Further, a representative sequence was selected to reduce the complexity of the data for downstream analyses. The procedure to select representatives from the transcriptome followed the procedure described by . Contigs were classified into 22-genome orthogroup classification using HMMScan, and one longest sequence was selected from a set of unigenes that fall within the same Trinity component and the same orthogroup. A final cleaned transcriptome assembly was generated to contain 167,058 representative sequences that retain sequences with open reading frames (ORF) greater than 200 bp. The coding sequences (CDS) and peptide sequences (PEP) were predicted using ESTScan. Sequences were further classified into the 22genome orthogroup classification  using the Gene Family Classifier pipeline in PlantTribes (https://github. com/dePamphilis/PlantTribes). The putative function of each sequence was analyzed with InterproScan, Pfam analyses, and annotation information (including the best Arabidopsis hit, GO information) from its assigned orthogroup.

Read mapping and quantification
Filtered clean reads from each library were mapped onto the final cleaned assembly using bowtie2 (Langmead and Salzberg, 2012) with default parameters. A gff file with the length of each transcript was generated, combined with the mapping output BAM file, expression quantification was performed with featureCounts (in SubREAD (Liao et al., 2014). Reads were combined by including paired-end reads and single reads using featureCounts (Liao et al., 2014). A gene expression file was generated to contain read counts for each gene from each library. The gene expression file was then filtered to remove contigs with mean read count across all the libraries less than 1.
Read counts mapped in each library were then normalized by fragments by kilobase per million mapped reads (FPKM) (Mortazavi et al., 2008) that accounts for bias in contig length as well as the sequencing depth from each library.

Sample clustering with pvclust and principal component analyses
Problematic libraries identified using pvclust (Suzuki and Shimodaira, 2006) with the assumption that replicates within one sample should group with each other. Bootstrapped sampling clustering was performed by the pvclust (Suzuki and Shimodaira, 2006) command in R using the gene expression matrix that quantifies the expression of each gene in log2FPKM across each library. Correlation distance and complete clustering were adopted. Replicated libraries that didn't group with the other replicates within the same sample were discarded. PCA was performed with the rlog transformed read counts by rlogTransofrmation in DESeq2 and prcomp function in R. The code for PCA analyses is also available on GitHub.

Identification of differentially expressed genes
DE analyses were performed between any two samples using DESeq2 (Love et al., 2014) where the rlog transformation (rlogTransofrmation function in R) significantly increases the power of DE identification. Read counts of all genes from each two tissues were used for DE analysis. A P-value cutoff of 0.01 and a fold-change cutoff of 2 were required for DE genes.

Clustering analyses
We utilized an optimized hierarchical clustering (Bar-Joseph et al., 2001) using a customized R function for clustering analyses. The optimized hierarchical clustering is preferred as it adopts an improved measurement of gene-and-gene distancecorrelation distance, instead of Euclidean distance (Bar-Joseph et al., 2001). The advantage is that it identifies more tightly linked cluster of genes with similar expression pattern across all stages regardless of their absolute level of expression. The gene expression matrix in log2FPKM was used as the input of hierarchical clustering. Then the top 4000 genes with the median absolute deviation (mad function in R) were used for optimized hierarchical clustering. Overall gene expression clustering was then examined, and the number of clusters was selected based on the transcriptional pattern. Each cluster of interest was identified by applying function "hclust" and "cuttree" function.

Identification of known floral and meiosis genes
Using a comprehensive literature search from the published functional characterization studies on flowering process in Arabidopsis and rice, we curated a list of genes involved in floral development regulatory network (Supplementary Table  S3), as well as in meiosis (Supplementary Table S2). We also expanded the list of meiosis genes by adding a set of genes with tissue specific expression in pollen mother cells from an RNA-Seq study by Tang et al. (2010), which we term as PMC genes.

Identification of unannotated genes in Phyllostachys
Several orthogroups miss Phyllostachys genes in the annotated gene models, such as orthogroup 7224 (LEAFY) (Weigel et al., 1992), orthogroup 2086 (FT), the absence of which is likely due to poor annotation. To recover the unannotated genes in Phyllostachys, we performed a TBLASTN search against the Phyllostachys genomic scaffolds using the predicted peptide sequences either in B. oldhamii 'Xia Zao' or rice. We then identified the genomic hit in Phyllostachys ) (by creating a gff file from the aligned boundary and used "fastaFromBed" function to get the genomic sequence), which was then aligned with the B. oldhamii 'Xia Zao' query CDS sequences to identify possible exon-intron boundaries. Alternatively, if this method failed to identify clear boundary and if the gene of interest is quite conserved, we used the full-length CDS sequence of rice to identify the genomic hits in Phyllostachys. Then the exonic regions of Phyllostachys sequences were used to predict possible peptide sequences with all possible six frames; the correct frame was determined as the one that well aligned with the predicted peptide sequences in B. oldhamii 'Xia Zao' or Oryza.
Sample fixation, paraffin section preparation and sample staining Sample fixation and preparation followed the procedure described by Li (1991). Samples were taken and directly fixed in formaldehyde acetic acid (FAA: 35-40% formaldehyde, 70% alcohol and 100% acetic acid in a proportion of 5:90:5 v/v) at 4°C until preparation. At the start of preparation, the samples were treated with Ehrfich's haematoxylin as a whole, then washed in running water. After the samples became blue, 50%, 70%, 85%, 95% and 100% ethanol were used successively (2×1.5-2 h) and then kept for 1h in solution of xylol and ethanol in proportions 1:1, next treated with pure xylol for 1 h (2×), and then kept 2-3 h in solution of xylol and paraffin in proportions 1:3, 1:1 and 3:1 successively, placed in a mixture of xylol-paraffin and kept at 58°C, stored overnight allowing xylol to evaporate. Then the samples were kept in pure paraffin for 2-3 h (3×). Finally, the samples were embedded in paraffin. Transverse sections (7mm) were cut using a rotary microtome (Leica RM2255) and double stained with 1% alcoholic Safranin O (Sigma S-2255) (in 50% ethanol) and 1% Fast green (Fluka 05500). The sections were permanently mounted in Canada balsam. Throughout the observations, images were captured with a microscope (Leica DM5000B). At least 3 replications were caried out for each experiment and micrograph.
RNA extraction, cDNA synthesis and reverse transcription quantitative real-time PCR Total RNA was extracted from the collected samples according to the manufacturer's instructions (RNA pure Ultrapure Total RNA Rapid Extraction Kit, Yuanpinghao Bio). The purity/concentration of RNA samples was determined using a NanoDrop 2000c spectrophotometer (NanoDrop, Thermo Scientific). Samples with concentrations greater than 100 ng.ml -1 and an optical density absorption ratio A260/A280 greater than 1.8 were used for cDNA synthesis and the integrity of RNA bands were detected by 1.0% agarose gel electrophoresis. Individual RNA samples were stored at -70°C and then 1 mg of total RNA were used as template in RT reactions with the PrimeScript RTase reverse transcriptase (TaKaRa), according to the manufacturer's instructions.
The primers for 12 genes from B. oldhamii 'Xia Zao' were designed by using the Primer 5.0 software, referring to the Supplementary Table S5. (http://www.genome.wi.mit.edu/cgibin/primer/primer5.cgi). All primer pairs were initially tested via standard RT-PCR using the Premix Ex Taq (TaKaRa) and the presence of a single amplification product of the expected size for each gene was verified by electrophoresis on a 2.5% agarose gel. cDNA products for qPCR were diluted 10-fold prior to use in real-time PCR. To ensure the absence of contamination or primer-dimer formation for each primer pair, a template-free control reaction was run. qRT-PCR reactions were carried out in 96-well blocks with an Applied Biosystems StepOne Plus Real-Time PCR system using TransStart Tip Green qPCR SuperMix (TransGen Bio). The reaction conditions were recommended by the manufacturer (30 s at 95°C, 40 cycles of 95°C for 5 s, and 55-60°C for 30 s). The dissociation curve was obtained by heating the amplicon from 60 to 95°C. All qRT-PCR reactions were carried out in triplicates. Each reaction was repeated three times, and the expression levels were calculated using the 2 −DCT method (Livak and Schmittgen, 2001). Finally, RNA-seqderived expression levels were compared with those derived from qRT-PCR.

Results
Microspore and male gametophyte development and anther cytological abortion of B. oldhamii 'Xia Zao' Each floret of B. oldhamii 'Xia Zao' has six anthers, and each anther has 4 sporangiates. At the stage when the spikelet length was less than 0.5 cm, the anther was still in an undifferentiated state, a bulky differentiated sporogonium can be seen beneath the epidermis cells at the corner of the anther on the transverse section ( Figure 2). Porogonium divides periclinally to produce a primary parietal cell and a primary sporogenous cell ( Figure 2). During the early developmental stage of secondary sporogenous cells, there were only a few secondary sporogenous cells identified in loculus ( Figure 2). The young anther walls usually consisted of three layers: from the outer to the inner layers, epidermis, endothecium and secondary parietal cell layer (Figure 2), at this time the spikelet length was < 1.0 cm (stage FI1). The structure of anther wall was fully developed in the microsporocyte stage. Outside-in layers of the anther wall, the epidermis, endothecium, middle layer and tapetum, can be seen in specimen (Figure 2). At this point, the length of spikelet grows to 1.5 cm (stage FI2). When the spikelet length was > 1.5 cm (stage FI3), microsporocytes would go into a period of meiosis. The abnormal meiotic behaviors such as a denser cytoplasm, abnormal spindle fiber and unformed equatorial plate were frequently observed (Figure 2). At the stage of microspore development, degeneration was still observed in the anther such as unclear nucleus in pollen grains and distorted tapetum (Figure 2). At the vacuolated uninucleate microspore stage, some microspores shrank seriously, resulting in deformation and degeneration of these microspores (Figure 2). Even when the 2-cell pollen grains were formed, their pollen walls were shrunk and aberrant, and the residual tapetum could still be observed (Figure 2). In some anthers, almost all microspores deformed, and no uninucleate microspores could successfully develop into pollen grains.

Sample collection and assembly statistics
More than 1.5 billion reads were generated from six major developing tissues during bamboo floral development (Figure 1) and resulted in 942051 contigs in the Trinity assembly. These tissues include culm base buds from flowering (F1_Cb) and nonflowering plants (NF3_Cb), underground (F1_ShU) and above-ground shoot buds (F1_ShA), and florets of different size (F1_Fl1, F1_Fl2, F1_Fl3) representing different stages of flower development. After filtering and selecting representatives (following Method "Assembly, cleaning, and annotation"), 167,058 contigs with an average unigene length of 766 bp and an N50 unigene length of 380 bp were represented in the final cleaned non-redundant assembly (Table 1). To examine if shorter contigs were attributed to possibly artificial unigenes constructed from traversing the de brujin graph path in Trinity assembly, we only considered 78,568 contigs that are assigned to an orthogroup (Table 1). The N50 statistics showed a two-fold increase of 657 bp, and the average sequence length was 1056 bp (Table 1). To evaluate the completeness of the transcriptome assembly, we examined the capture rate of the ultra-conserved orthologs (UCOs) gene set. The gene capture rate of 99.1% (Table 2)  ( Suzuki and Shimodaira, 2006). Because they are not grouped with the other replicate libraries of the same tissue, it represented some inconsistency within-sample replicates that could be due to sample harvest, RNA extraction, or sequencing process, and was valid for removal. The global transcriptional profile by the PCA plot reveals that the biggest variation between the samples lies between clump 1/3 and clump 2 tissues, as the first principal component separating clump 1 (together with clump 3 culm buds) and clump 2 floral tissues explains 70.4% of the variation (Supplementary Figure S1). The second principal component shows that floral stage 2 and floral stage 3 are most distinct from the remaining vegetative tissues and floral stage 1 tissues from clump 1 and nonflowering clump 1 culm buds. Interestingly, culm buds' tissue from flowering clump 1 and nonflowering clump 3 are more similar to each other, indicating a similar transcriptional profile within the same culm buds tissue despite from different clones (Figure 3 and Supplementary Figure S1). The unexpected grouping of F1_Fl1_R3 and F2_Fl1_R1 (Figure 3) again support these two samples are possible outliers. The three replicates of flower 2 (blue) and flower 3 (light green) in clump 1 show tight correlation. Conversely, the clump 2 flower tissues show greater within-sample correlation than between sample variation. In addition, there is a higher between-clone variation in reproductive tissues (between clump 1 and clump 2) than in vegetative culm buds' tissues (between clump 1 and clump 3). Because of this high level of variation from clump 2, samples from clump 1 (6 tissues that start with "F1") and clump 3 (one tissue that start with "NF") were focused for downstream differential expression and clustering analyses in flowering evolution of B. oldhamii 'Xia Zao'. PCA plot of the 27 samples of B. oldhamii 'Xia Zao' generated for this transcriptome study. The samples include one culm buds tissue of three biological replicates from nonflowering clump 3labeled as "NF3_Cb", six tissues each with three biological replicates from flowering clump 1labeled as "F1", as well as three tissues each with two biological replicates from flowering clump 2labeled as "F2". "Cb"culm buds, "ShU"underground shoots, "ShA"aboveground shoots, "Fl1" -flower stage 1 with each spikelet length less than 1 cm, "Fl2" -flower stage 2 with each spikelet length between 1cm and 1.5 cm, "Fl3" -flower stage 3 with each spikelet length greater than 1.5 cm. spikelet.

Validation of digital gene expression by qRT-PCR
To confirm the genes in different development samplings, ten genes were randomly selected for qRT-PCR. Expression of four genes (Hd3a, MADS56, E class gene, MADS15) by qRT-PCR was fitted well with the pattern of Tag-seq analysis (Figures 4A-D). And another eight genes (B class, C class, E class, ZIP4, TDR, PTC1, CYP703A3, G class) similarly behaved between qRT-PCR and deep-sequencing methods (Supplementary Figure S2).

Differential gene expression analyses
Pairwise DE analyses were performed between any two samples of the seven tissues from the nonflowering clump 3 and flowering clump 1. As a result, 21 DE analyses were generated that identified a number of DE genes between any of the seven tissues. The biggest variation is revealed almost between clump 3 and clump 1 tissues, except between the culm buds' tissues from the two (Supplementary Figure S3). It suggests that the variation between clump is smaller than the variation between tissue from the same clump. Within clump 1, the variation is greatest between culm buds and stage-3 floral tissues (12,587 DE unigenes). Significant difference was also revealed between stage-1 floral tissues and stage-3 floral tissues (7712 DE genes).

Abundant floral tissue expression accounts for the major transcriptional profile
Our optimized clustering identifies a total of nine global transcriptional clusters (Supplementary Figure S4 and Table  S1) based on the top 4000 genes that show the highest expression deviation from the median. Cluster 1 shows a dominant signature of flower gene expression with upreuglated expression in all three flower stages, but quite low level of expression in vegetative tissues (culm buds, undeground and above ground shoot buds). Cluster 2, however, shows an opposite pattern of vegetative tissue specific expression. Cluster 3 appears to capture genes specifically expressed in shoot, likely contributed by genes involved in rapid cell expansion. Each cluster is represented by a different set of genes with the numbers vary across clusters. The major expression pattern is characterized with abundant expression in flower stage 2 and stage 3 in cluster 8 ( Figure 5A), indicating a large number of genes contribute to the floral transcriptome. qRT-PCR validation of the expression for four selected genes in different developmental stages. The transcript levels were normalized to that of PP2A and UBC18.
To identify if unigenes from cluster 8 encode the genes involved in meiosis development, we did an enrichment analysis using this set of genes. Protein sequences were used to blast search against the rice peptide sequences to obtain the best hits, which were used for Gene Ontology enrichment analysis using the webtool (http://www.ricearray.org/analysis/overview.shtml). As expected, several GO Biological process (GO BP) categories related to meiotic development were significantly enriched among this set of genes, which include inflorescence development, pollen development, anther morphogenesis, anther development, tapetal cell differentiation, and tapetal layer formation. These categories were mapped onto six genes with a known role in flower regulation and development in rice, which include two genes (LOC_Os02g02820 (AMS in Arabidopsis) (Sorensen et al., 2003) and LOC_Os07g36460) encoding a helixloop-helix DNA binding transcription factor (TF)) and one encoding a FT-like 2 homologous to Flowering locus T gene (LOC_Os06g06320, Hd3a from orthogroup 2086) ( Figure 5).

Misregulated transcriptional profile likely associated with meiotic failure
As morphology of B. oldhamii 'Xia Zao' meiotic processes indicates evidences of meiotic abortion (Figure 2), we attempted to look for any abnormal patterns related to expression. From the nine overall transcriptional clusters (Supplementary Figure S4), we identified three clusterscluster 5, cluster 7, and cluster 9 that could be associated with possible meiotic failure (Supplementary Figure S5). All these three clusters display high expression in at least one flower stages, yet show decreased expression towards flower 2 or flower 3 stage (Supplementary Figures S5A-C), indicating a failure to complete the meiotic development. However, the number of genes from these three clusters (682) accounts for only a small number of genes, especially when it comes to comparison with the 1115 unigenes from the major transcriptional cluster 8. Nevertheless, this indicates a possibility of misregulation of gene expression which could be associated with meiotic failure. Gene enrichment analysis identified a distinct set of functional enriched terms from the major cluster 8 (Supplementary Figure S5D), suggesting genes involved in processes such as meiosis, megasporogenesis, and tapetal layer development may be misregulated to cause possible meiotic failure, whereas genes related to enriched functional terms from cluster 8 such as pollen development don't appear to be responsible for meiotic failure of B. oldhamii 'Xia Zao'.
Conserved AMS sequence evolution in B. oldhamii 'Xia Zao' As one of the enriched flower associated functional categories from cluster 8 is anther development and AMS (LOC_Os02g02820 in Figure 5) is a key gene involved in this process (Sorensen et al., 2003), we examined if AMS shows any nonconserved sequence evolution in B. oldhamii 'Xia Zao'. Orthogroup peptide multiple sequence alignments (MSA) shows three fragmented AMS contigs, two align on the N-termiinal, and one aligns on the C-terminal, and a gap between the three sequences indicate failure to assemble the full-length AMS (Figure 6). To examine whether this was due to assembly artifact, we aligned the rice AMS CDS with these three B. oldhamii 'Xia Zao' AMS contigs. By also searching aginst the original trinity assembly, and also by filling the gaps (gap mainly exists in exon 5) using raw read sequence data (Figure 6), we managed to construct the full-length AMS sequence which contains seven exons in the CDS (Figure 6). In the end, the full-length CDS of AMS encodes a quite conserved protein with a length of 547 AA. The protein sequences of AMS in B. oldhamii 'Xia Zao' and rice show 64.3% identity (Figure 6 and Supplementary Table S4), and CDS sequences show 73% identity ( Figure 6 and Supplementary Table S4). The correct placement of B. oldhamii 'Xia Zao' AMS gene as a sister to its ortholog in Phyllostachys and rice indicates conserved protein sequence evolution for AMS gene on the sequence level (Supplementary Figure S6A). This indicates high level of complexity in transcriptome assembly, and for conserved genes, manual curation is possible to create the full-length coding sequences. In conclusion, AMS is not missing in B. oldhamii 'Xia Zao' even though evidences suggest meiosis don't function properly in B. oldhamii 'Xia Zao' (Figure 2). However, downregulation of AMS in B. oldhamii 'Xia Zao' flower 3 stage ( Figure 6) is different from nonreduced expression of rice AMS (Figure 6), suggesting that a noncanonical gene regulation might be associated with late stage of meiosis.
Nonconserved sequence evolution for PMS gene might be associated with meiotic failure In examining sequence evolution for all meiotic orthogroups, we identified an 8-bp insertion in exon 3 of B. oldhamii 'Xia Zao' PMS1 gene (orthogroup 8661), the fulllength CDS of which was constructed by comparing transcript comp149596_c0_seq1 with rice CDS. Alignment between B. oldhamii 'Xia Zao' PMS1 transcript with rice PMS CDS reveals a nonspliced 695bp-intron8 in B. oldhamii 'Xia Zao' CDS ( Figure 7 and Supplementary Figure S7). The constructed full-length contains 12 exons in the CDS and high level of sequence similarity of 84.7% identity with rice CDS. Close examination of the two sequences reveals an 8-bp insertion in exon 3 of PMS, which is supported by reads, and has disrupted the frame by 2 bp and causes a premature stop codon resulting in a much truncated protein. It is believed that this small insertion may be associated with meiotic failure in B. oldhamii 'Xia Zao', as this gene plays a role in DNA mismatch repair and mutation in Arabidopsis shows significant decrease in seed The major transcriptional profile includes genes highly expressed in flower tissues and functional enrichment shows the enriched Gene Ontology Biological Process terms. (A, B) show the expression profile of genes in this cluster, and (C) shows the enriched functional categories and representative gene in rice. Zhao et al. 10.3389/fpls.2022.1023240 Frontiers in Plant Science frontiersin.org production , which is consistent with the reduced seed production in B. oldhamii 'Xia Zao'.
Identify floral initiation genes in early stages of B. oldhamii 'Xia Zao' development One objective of B. oldhamii 'Xia Zao' floral transcriptome is to investigate if it shares the same floral initiation pathway as other flowering plants such as Arabidopsis and rice. To identify if floral initiation genes were already expressed in early developmental stages such as culm buds, we performed differentially expressed (DE) analyses between nonflowering culm buds and flowering culm buds and showed that a large number of selected genes known to function in pollen developmental and meiosis were significantly enriched at early flowering stage (F_FI1) (Supplementary Figure S8). A total number of 3080 contigs from 1953 orthogroups defined by 22 sequenced plant genomes ( Figure 8) were identified, which we cross-checked with our curated list of floral genes (Supplementary Table S3). We claim a B. oldhamii gene an ortholog of known floral gene when it is assigned to the same orthogroups for known floral genes from rice or Arabidopsis. This identified 17 known floral orthogroups from the DE transcripts between nonflowering and flowering culm buds (Figure 8). Only four contigs showed a reduced expression in flowering culm buds compared to nonflowering culm buds, with A schema illustrates the process of full-length CDS construction. Exon 1, exon 2, exon 3, exon 4, exon 6, exon 7 were constructed by high sequence similarity in assembled contigs with rice CDS. The left region of exon 5 up to 784 bp of rice AMS CDS and the right region of exon 5 from 911 to 939 were constructed by comparing B. oldhamii 'Xia Zao' transcriptome assemblies to rice CDS. The gap in the middle of exon 5 (orange region labeled as "extended-by-reads") were filled by using read sequence data. In this process, three seeds were used to perform "grep" command in the reads to extend the right regions, and three iterations ("grep2extend", "grep2extend-iter2", and "grep2extend-iter3") were able to fill the gap in exon 5. (E, F) are the pairwise alignment (from Geneious) of full-length AMS CDS and peptide sequence between B. oldhamii 'Xia Zao' and rice.
a majority showing an upregulated expression in flowering culm buds ( Figure 8). These floral regulators that were expressed as early as culm buds include genes such as LEAFY (orthogrup 7224), AP1 which interacts with LEAFY (Liljegren et al., 1999), VIN3 (Kim and Sung, 2013) and GA1 (Wilson et al., 1992) ( Figure 8), a group of genes with known roles in integrating various environmental factors to regulate flowering. By using Phyllostachys genome and conserved sequences in rice, we i d e n t ifi ed the full-leng th coding sequence of FT (Supplementary Table S4), LEAFY (Supplementary Table S4), two genes with well-known roles in flower initiation. Surprisingly, instead of relaxed sequence constraint (Supplementary Figures S6B, C and Table S4), they all retain high level of sequence conservation. This suggests the retention of key processes involved in control of flowering time in at least two sequenced bamboo species, B. oldhamii 'Xia Zao' and Ph. edulis. Our gene expression network analysis was based on the sequence assembled by the de novo transcriptome, and we didn't really rely on for the association analysis. We mainly focused on the overall expression pattern of the gene. Additionally, we conducted a large number of comparative studies with reference to the known flowering time network of rice, and the reason we selected rice over bamboo is that the flowering pathway research mechanism in rice is relatively unambiguous. We pay attention to several important ways to regulate flowering time, such as photoperiod pathway, autonomous pathway, vernalization pathway and so on. We identified the differentially expressed genes in flowering and nonflowering tissues by blast to find their best hit in rice, and then compared these flowering pathways are conserved in bamboo. The results showed that only genes in the photoperiod pathway were found to be differently expressed in flowering and nonflowering tissues. The significance of this work is that we have demonstrated that the flowering pathway of bamboo may be complete, and the reason it is so difficult to bloom is due to genetic regulation. Our sampling and analysis also show that the photoperiod pathway may play a major role in controlling the response of bamboo to flowering or not (Supplementary Figure S9).

Discussion
There are no seeds or fruits have ever been found after flowering in B. oldhamii 'Xia Zao'. In this study, various stages of anther development were discerned with microscopic observation of transverse sections. There are several types of abortions observed in the development phases of anthers in B. oldhamii 'Xia Zao'. For example, abnormal diad cells were observed during meiosis. Similar phenomenon was also reported in other bamboo species, such as Menstruocalamus sichuanensis  and Shibateae chinensis (Lin and Ding, 2012).  also reported that the cease of development led to the deformation of uninucleate pollen in one cultivars of sterile chrysanthemum. Besides, at the stage of microspore, distorted tapetum was still observed in the anther. Therefore, it could be concluded that the aberrant tapetal behaviors are associated with male sterility Meric et al., 2003;Shi et al., 2010). Holford et al. (1991) summarized the following three types of abnormal tapetal behavior in a male sterile onion: the premature breakdown of the tapetum at tetrad stage; the hypertrophy of the tapetum after the diad stage followed by its premature autolysis; the tapetum remaining in good condition but for an abnormal long period of time. Buyukkartal et al. (2005) considered that persistent tapetum blocked nutrient transportation and then led to low pollen fertility and affected the proportion of seed setting and grain formation. Therefore, abnormal tapetal cells caused anthers abortion, which is a common phenomenon in B A FIGURE 7 Gene structure schema of PMS1 (A) and NEF1 (B) gene from B. oldhamii 'Xia Zao' and rice reveals nonsense mutations causing premature stop codons. Insertion is marked byblue arrow and orange label, stop codon is indicated by a blue arrow. Zhao et al. 10.3389/fpls.2022.1023240 Frontiers in Plant Science frontiersin.org plants. In the study, the anther began to abort at stage FI3 in B. oldhamii 'Xia Zao', accompanied by unsuccessful meiosis of microspore mother cells and abnormal development of the tapetum, resulting in pollen abortion. Even when the 2-cell pollen grains were formed, it could also be seen that microspores and tapetal cells degenerated together in B. oldhamii 'Xia Zao'. This phenomenon provides a clue that the aberrant development of meiosis, tapetal cells, and microspores is associated with low seed setting rate in bamboo plants.
In this study, we identified a group of genes involved in meiosis, megasporogenesis, and tapetum development, which could be responsible for the abortion phenomenon during meiosis and pollen maturation. For example, a number of known genes involved in meiosis (Supplementary Table S2) were also identified in the study, such as TFs AtTGA10 (orthogrouop 4636) and AtSCC2 (orthogroup 5234), genes encoding mismatch repair proteins, AtMSH4 (orthogroup 5792) and AtRAD50 (orthogroup 9771), a meiotic specific protein AtZIP4 (orthogroup 9504), and a gene in structure maintenance of chromosome AtMSC4 (orthogroup 4970). The fact that many of these meiotic genes maintain high sequence conservation suggests they play a similar functional role in bamboo plants. Significantly, AMS is identified as a key gene involved in meiosis. AMS encodes a basic helix-loop-helix (bHLH) TF, which plays a pivotal role in tapetum development and microspore formation. It was reported that mutations in this gene can lead to tapetum abnormalities and microspore degeneration (Sorensen et al., 2003). For instance, AMS mutations downregulated the expression of the genes involved in meiosis in Arabidopsis, and many genes related to tapetum metabolism were also directly regulated by AMS. In addition, AMS may change mitochondrial permeability and induce tapetum PCD by inducing the expression of LTP family genes. For example, rice TdR, an AMS homologous gene, can directly regulate cysteine proteases related to PCD (Na et al., 2006). In this study, the AMS gene from B. oldhamii 'Xia Zao' has a high sequence similarity with that from rice ( Figure 6). However, the down-regulation of AMS in flowering stage of B. oldhamii 'Xia Zao' (Figure 6) is different from the non-reduced expression of AMS in rice (Figure 6), suggesting that a noncanonical gene regulation might occur in the late stage of meiosis.
Previous studies have shown that MUTL proteins are involved in DNA mismatch repair and play an important role in maintaining the stability of the meiotic genome (Kolodner and Marsischky, 1999;Hoffmann and Borts, 2004). Low seed setting rate is a widespread characteristic of MutL homologues mutants in eukaryotes. Evidence from Arabidopsis confirm that the deletion of AtMLH1 leads to a significant reduction in pollen collapse and fruit set (Dion et al., 2007). Jackson et al. (2006) also reported that two base insertions from mutant alleles of AtMLH3 resulted in truncation of the AtMLH3 protein, leading to significant decreases in fruit set in Arabidopsis. Kim et al. (2015) reported that AcPMS1 may play a restorative role in CMS onion. In yeast, PMS1, MLH1, MLH2, and MLH3 have been reported to be MUTL protein homologs, with PMS1 and MLH1 forming the major dimer involved in the repair of single mismatched bases (Bray and West, 2005). In this study, our comparative sequencing analysis revealed high sequence similarity of PMS CDS between B. oldhamii 'Xia Zao' and rice ( Figure 7A). However, the insertion of eight bases in exon 3 of PMS in B. oldhamii 'Xia Zao' caused a shift mutation that resulted in a scrambled frame and led to a premature stop codon at alignment position 606 ( Figure 7A). Therefore, we speculate that the premature termination of PMS1 contributed to the failure of meiosis in B. oldhamii 'Xia Zao'.
In addition to PMS, we are concerned that another gene that plays a major role in the regulation of pollen wall formation, NEF1, also has a nonsense mutation in B. oldhamii 'Xia Zao' ( Figure 7B). Pollen exine is an essential component of male gametophyte protection, defense, and female gametophyte recognition NEF1 encodes a membrane protein present in the plastids of the tapetum that may be involved in the transport of polysaccharide required for the synthesis of the primexine (Ariizumi et al., 2004;Wang., 2013). In Arabidopsis, NEF1 was reported to have a significant role in regulating the formation of pollen exine, and the absence of this gene will prompt to a deficit in primary exine synthesis, sporopollen accumulation failures, and the formation of fragmented microspores (Ariizumi et al., 2004;Dong et al., 2005;Chang et al., 2012). Besides, NEF1 has a high degree of evolutionary conservation, with many plant genomes containing a single copy of the NEF1 homolog, and this conserved structure ensures similarity in gene function (Wang, 2013). Only a single copy of NEF1 was found in B. oldhamii 'Xia Zao', and its sequence is highly homologous to rice ( Figure 7B). However, the presence of a base insertion behavior in the CDS region of NEF1 in B. oldhamii 'Xia Zao' that distinguishes it from the rice CDS may be responsible for the distortion and rupture of the outer pollen wall (Figures 2H, I and Supplementary Figure S10).
Compared to other flowering plants, bamboo has a long nutritional growth and an irregular flowering period (Zheng et al., 2020), and we know little about the molecular mechanisms that regulate flowering time in bamboo. In this study, we have therefore analysed the transcriptome differences between nonflowering and flowering tissues of B. oldhamii 'Xia Zao' by high-throughput sequencing, leading to the identification of many genes involved in candidate flowering pathways in B. oldhamii 'Xia Zao'. The evidence has revealed that only the genes involved in photoperiodic pathway are differentially expressed in the different tissues (Supplementary Figure S9). In Arabidopsis, three key genes GIGANTEA (GI), CONSTANS (CO) and FT in photoperiod pathway were reported to form a regulatory pathway to control photoperiod response and promote flowering in Arabidopsis (Suaŕez-Loṕez et al., 2001;Yanovsky and Kay, 2002;Putterill et al., 2004;Turck et al., 2008). There is a similar photoperiod pathway in rice, which is composed of GIGANTEA (OsGI), Heading date 1 (Hd1) and Heading date 3a (Hd3a), which are homologous to GI, CO and FT, respectively (Yano et al., 2000;Shoko et al., 2002;. In this study, we screened the photoperiod flowering network of B. oldhamii 'Xia Zao' and detected six photoperiod pathway node genes homologous to rice: LEAFY, OsMADS56, Hd3a/RFT1 and OsMADS14/OsMADS15 (Figure 9). Except for OsMADS56, the other five genes were highly expressed in flowering culm buds of B. oldhamii 'Xia Zao'. FT gene acts as a key floral pathway integrator, essential for determining flowering time. It acts as an important flower-forming element that is transported to the SAM and thus induces flowering in plants (Corbesier et al., 2007;Tsuji et al., 2013). Hd3a and RFT1, which are highly homologous to FT as florigenic genes under short-and long-day conditions, respectively, are the main flowering inducers in rice (Shoko et al., 2002;Komiya et al., 2008). OsMADS14 and OsMADS15, two genes AP1-like genes, located downstream of flower-forming elements, have been shown to be positive regulators of flowering in rice (Kobayashi et al., 2012). LEAFY, a transcriptional regulator that determines the formation of floral meristem, is highly conserved in different plant species (Weigel et al., 1992). Most known LEAFY genes are expressed in the vegetative and floral meristem (Weigel et al., 1992;Carmona, 2002;Wang et al., 2012;Yang et al., 2017;Hu et al., 2020). For example, homologous gene of LEAFY, RFL was expressed in only young panicles in rice (Kyozuka et al., 1998). LiLFY1 was expressed in only young flower buds and apical meristems in lily (Wang et al., 2008). In this study, we assembled the full-length coding sequences of FT and LFY using the genome of B. oldhamii 'Xia Zao', which show a high degree of conservation with rice. However, the expression pattern of the genes displays a completely opposite trend in the two plant species. This suggests a degree of mis regulation of flowering-related genes may be related to the peculiar flowering characteristics of bamboo plants.

Conclusion
Bamboo species represents a good evolutionary model to study flowering evolution, due to its predominant asexual reproduction using clonal propagation. We thus hypothesize Mapping of DE genes in flowering and nonflowering culm buds in B. oldhamii 'Xia Zao' reveals the conserved evolution of photoperiod pathway. The yellow color represents integrated node gene. that its flowering pathway must have relaxed to show nonconserved sequence evolution. This signature identified a cluster of genes including FT and AMS, two known key regulators of flowering and meiotic development. We examined in great detail of the sequence and expression evolution of these two genes, as well as another floral regulator, LEAFY, by using the genome of moso bamboo and conserved sequences in rice. Surprisingly, instead of relaxed sequence constraint, the full-length coding sequence of these three genes all retain high level of sequence conservation. LEAFY also shows its induced expression in flowering culm buds compared to the culm buds of nonflowering bamboo, indicating the presence of a functional flower and meiotic pathway for at least these three genes. In addition, by comparing to the documented transcriptional profile in rice, all three genes show an opposite expression trend. These lines of evidence indicate that gene regulation instead of protein sequence evolution may underlie the reduced flowering and meiotic failure in B. oldhamii 'Xia Zao'. In addition to these master regulators, our analyses of sequence evolution revealed quite conserved sequence evolution for additional meiotic genes including TFs (AtTAG10, AtZIP4) and genes involved in mismatch repair (AtRAD50, AtMSH4). However, two genes -PMS1 and NEF1 show INDEL structural variation that disrupt the frame, and the high level of sequence conservation in the CDS with rice suggests that they are rather recent events. The processes these two genes contribute to and the associated meiotic abortion indicate a potential link between them and reduced seed production in B. oldhamii 'Xia Zao'.

Data availability statement
The data presented in the study are deposited in the NCBI SRA repository, accession number PRJNA876039.

Author contributions
SL designed the research. WZ, LZ, and WY carried out the experiment. ZY and CG performed the bioinformatics analyses of the transcriptome data. WZ, CG, and WY wrote the manuscript. ZY and YD revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding
This work was financially supported by the National Key Research & Development Program of China (2021YFD2200503), the National Natural Science Foundation of China (Grant No. 31870595;No. 32001292) and the Priority Academic Program Development of Jiangsu Higher Education Institutions.