text:
Sustainable Applications for the Valorization of Cereal Processing By-Products

This review article revises the sustainable practices and applications to valorize valuable components recovered from cereal processing by-products. After introducing cereal processing by-products, their healthy compounds, and corresponding functional properties, the article explores reutilization opportunities of by-products emphasizing specific sources (e.g., oat and wheat bran, distillers’ dried grains, etc.) and the biorefinery approach. Proteins and soluble dietary fibers such as arabinoxylans are of particular interest due to their content in the cereal processing by-products and their easy extraction based on conventional technologies such as enzyme-assisted extraction and membrane filtration. Non-thermal technologies have also been suggested to improve sustainability recovery approaches. Finally, the article discusses the different applications for the recovered high-added value compounds that span across biotechnology, foods, and bakery products.

Introduction
Cereals comprise an essential source in the human diet and a significant part of livestock feed for thousands of years, while their processing represents a substantial asset to the food production chain [1]. Epidemiological studies have indicated that the consumption of whole cereal grains is correlated with a lower risk of developing cardiovascular and chronic diseases such as cancers, diabetes, and obesity [2]. The health benefits of cereal grains have been attributed to their content in high amounts of lipids, proteins, dietary fiber, tocopherols, Vitamins B, and E [3]. On the other hand, these compounds are concentrated in the hull, bran, and cereal germ that are removed from foods during cereal manufacturing.
Cereal manufacturing includes different processes such as dry and wet milling, malting, and pearling, which produce by-products of varying nature and chemical compositions, e.g., wheat and maize bran, rice bran and corn germ oil, and distiller's dried grains [4]. Nevertheless, all of them contain valuable nutritional components (similar to whole grains) that could be converted to biofuels, bioplastics and biopolymers. Alternatively, they could be recaptured and reused in the food chain, finding innovative nutraceutical and pharmaceutical applications [5] as well as applications in fermentative for the production of bioactive microbial metabolites, enzymes, single cell proteins and oils [6]. For example, alkylresorcinols (existing only in rye and wheat bran) could reduce cholesterol absorption, and reduce the risks of chronic diseases such as obesity and diabetes [7]. Besides, the market of functional food ingredients supporting the immune system will grow in the post-pandemic era [8,9]. At the same time, there is an urgent need to valorize food processing by-products in order to improve food security in the next decades [10]. Cereal processing by-products are rich in β-glucan that has been proposed as a bioactive food compound against SARS-CoV-2 infection [11].
Besides, the current managing practices that discharge cereal processing by-products to the environment are not sustainable. This article revises the valorization opportunities  Table 2 summarizes the functional components of cereal processing by-products. Rice bran high amounts of carbohydrates (~50%), oil (~20%), proteins (~13%), and dietary fiber (~11%) that is comprised of β-glucan, pectin, and gum. It also contains significant amounts of minerals such as magnesium, phosphorous, and iron, as well as bioactive phytochemicals such as ferulic acid, phytic acid, squalene, polycosanols, phytosterols, oryzanols, and tocotrienols [14,15]. On the other hand, corn bran is mainly composed of insoluble dietary fibers such as hemicelluloses (ca. 700 g/kg), cellulose (ca. 240 g/kg), and lignin fractions (ca. 10 g/kg) [16]. Compared to other cereals bran, corn bran has the highest dietary fiber content, tocopherols, and polyphenolic compounds that have well-known antioxidant properties and can be used in different applications. For example, they can be used as bioactive compounds in cosmetics or natural substitutes (antioxidant preservatives, stabilizers, emulsifiers, and coloring agents) in foods preventing the potential adverse effects associated with the consumption of the artificial ingredients [17][18][19][20]. Corn bran is a rich source of ferulic acid compared to other vegetables, fruits, and cereals, while sorghum bran contains high amounts of 3-deoxyanthocyanidins. The latest is a rare class of flavonoids with intense cytotoxic activities. Besides, the oil recovered from rice bran is rich in tocopherols, phytosterols, and tocotrienols and particularly γ-oryzanols, which show ten-fold higher antioxidant properties compared to tocopherols [20][21][22]. Table 2. Functional components of cereal processing by-products and health properties (remade from [5]).

Cereal
By-Product Target  The processing by-products of oat comprise the most typical examples of underestimated substrates with high valorization potential. Oats are rich in proteins (11-20 g/100 g) and β-glucan (2.2-7.8 g/100 g) [28,29], which have been correlated with the control of normal cholesterol's level in the blood [30]. This property of β-glucan and proteins as well as their viscoelastic characteristics allow considering them as ideal additives in foods and confectionary [31][32][33]. For instance, hydrolysis of oat's carbohydrates (e.g., starch, maltose, and β-glucan) has been implemented for the industrial development of nondairy formulations such as milk alternatives, cream, ice cream, and other products consumed by lactose-intolerant individuals [34][35][36]. The biodegradation process is conducted using multi enzymatic cultures, and the viscoelastic properties of the hydrolysates are optimized using kinetics modeling [37]. The by-product of this process (oat mill waste) is typically utilized as livestock feed, although it is still rich in β-glucan and proteins. Thus it could be used in different food applications, e.g., to mimic and replace fat of cheese and yogurt [38][39][40]. Sibakov et al. [41] applied a dry fractionation to remove lipids from oat grain and extract β-glucan from oat grains. This recovery approach produced another by-product rich in unfolded globulin that aggregates due to the acidic conditions [42,43]. Aggregated proteins show poor functionality, e.g., decreased solubility, viscosity, emulsification, and foaming ability [44,45]. On the other hand, the aggregation of proteins could be utilized for the structure formation of lactic acid fermented products. Thereby, it is essential not only to recover target compounds, but also to develop tailor-made applications [46].
Wheat bran is rich in dietary fiber and particularly arabinoxylans, which are non-starch polysaccharides with high ratios of arabinose/xylose β-glucans, fructans, cellulose, and lignin [25,47]. Arabinoxylans can be found in the cell walls of all cereal grains interlinked with ferulic and other phenolic acids [47]. They present different structures and properties, namely water retention, viscoelastic properties, high nutritional characteristics, and healthpromoting impacts, e.g., reduce the glucose levels in the blood [18]. Wheat bran's fiber contributes to several beneficial physiological effects (e.g., laxation, acceleration of intestinal transit, increase in fecal bulk, and blood cholesterol attenuation) [5,48]. More specifically, arabinoxylans can lower blood sugar levels and are widely used in food (e.g., as sweetener, medicine, and cosmetic industries) [49,50]. On the other hand, phenolic acids show high oxidation resistance and exert in vitro anti-inflammatory and anticancer effects [51].
Rye bran is also rich in dietary fiber-containing higher contents of β-glucan and fructan. Still, lower amounts of arabinoxylans than wheat (5%, 7%, and 23% compared to 2.5%, 3%, and 26%, respectively) [52]. Arabinoxylans have prebiotic properties (e.g., bifidogenic enhancement) and can be used as fat replacers and functional components in meat, dairy and bakery products [18,38,53,54]. The recovery of arabinoxylans from wheat, rye, and other cereals bran is typically conducted using several enzymatic, chemical, and hydrothermal treatments that lead to arabino-oligosaccharides with different lengths. Subsequently, the challenge during extraction is to match the structure and size of the recovered arabinoxylans with the required properties of the final products. For instance, the application of alkali or acid solutions causes the breakage of ferulic acid, leading to the extraction of oligosaccharides with a low arabinose substitution ratio [18].
Barley is the primary raw material for the production of beer. During malting and brewing, the main leftovers of the process (~85%) are the well-known brewers' spent grains (BSG). BSG contains high amounts of lignin (28%) and non-cellulosic polysaccharides (28%), including β-glucan and arabinoxylans, and cellulose (17%), as well as has considerable amounts of proteins and polyphenols [18]. BSG is merchandized as livestock feed, although it comprises an excellent source for recovering their food-grade components. The proteins of barley, wheat, corn, and rice processing by-products have better functional properties compared to those of the endosperm due to their content in essential amino acids and their valorization potential in the market of extruded and cereal-based baked products is high [18,55]. The fermentation of BSG with different lactic acid bacteria (e.g., with Leuconostoc pseudomesenteroides DSM20193 and Weissella confusa A16 as suggested by Koirola et al. [56]) has been examined in different applications such as yoghurt [57], pasta [58] and bread [59]. Table 3 presents different approaches for the recovery of valuable compounds from different cereal processing by-products. As shown above, cereal bran, germ, and husk are rich in functional macromolecules such as proteins, β-glucan, arabinoxylans, pectin, and smaller molecules such as ferulic acid and vitamins. At first, each of these by-products should be pre-treated using grinding, wet milling, or drying [46]. The second step requires the separation of macro-and micro-molecules using different techniques such as isoelectric solubilization or alcohol precipitation. Isoelectric solubilization/precipitation can separate proteins more selectively (due to their charge) compared to alcohol's addition that precipitates proteins together with soluble (β-glucan, arabinoxylans, and pectin) and insoluble (e.g., lignin and cellulose) dietary fiber [46]. Another approach for the separation of macro-and macro-molecules in cereal processing by-products is the utilization of membrane processes. For example, Patsioura et al. [32] used a simple ultrafiltration step and a 100 kDa-polysulfone membrane to concentrate β-glucan from oat processing waste in a cross-flow module and remove smaller molecules in the permeate. The proposed process operated at low transmembrane pressure (≤2 bar), but it was not selective enough to separate β-glucan from proteins. Table 3. Different approaches and technologies for the recovery of valuable compounds from different cereal processing by-products.

Target Compound Recovery Method Results References
Oat mill waste β-Glucan Polysulfone membrane was applied in the pilot cross-flow module for the ultrafiltration of β-glucan containing feeds (<600 mg/L) recovered from the industrial oat mill waste Two thirds (~67%) of β-glucan had been recovered [32] Rice by-products Peptides Combined treatment with proteases and sequential fractionation with one microfiltration (0.2 µm) and three polyethersulfone ultrafiltration (8, 5, and 4 kDa, respectively) steps The retentate with the highest peptide content was the one recovered with the 8 kDa membrane [60] Wheat bran Proteins The adherent endosperm was removed from wheat bran with brushing before being extracted using the alkaline treatment and isoelectric precipitation High nitrogen dispensability with minimum and maximum protein solubility pH values of 5.5 and 11.5, respectively.
[61] Extraction of arabinoxylans with a lower ratio of arabinose/xylose as a factor of alkali's power. [63] Rice bran Proteins Coupling alkaline extraction with microwave-assisted extraction The optimum conditions were found to be a power of 1000 W, a water solid-to-liquid ratio of 0.89 g bran/10 mL, and an extraction time of 90 s. Microwaves enhanced by 1.5 fold the protein yield of alkaline extraction.
[ 64] In order to increase selectivity, many researchers suggested the application of sequential ultrafiltration systems for the valorization of by-products derived from different cereal processing by-products (e.g., from corn, rice, and wheat). The advantage of combined membrane systems is the synergy resulting from their integration, with overall significant benefits in plant compactness, product quality, energy efficiency, and environmental impact [18]. For instance, Ferri et al. [60] used a combined treatment with proteases (Protamex, alcalase, and neutrase) and sequential membrane processes for this purpose (Figure 1a). Rice by-products were first subjected to enzymatic hydrolysis (ph = 7 and 60 • C for 2 h) before recovering the digestate with centrifugation and its sequential fractionation with one microfiltration (0.2 µm) and three polyethersulfone ultrafiltration (8,5, and 4 kDa, respectively) steps. The retentate with the highest peptide content was the one recovered with the 8 kDa membrane. Compared to polysulfone, the polyethersulfone membrane is less hydrophobic, allowing the rapid absorption of polar compounds (e.g., polyphenols) and the release of soluble dietary fiber in the permeate [40]. In another effort, Castro-Muñoz and Yáñez-Fernández [65] used a combination of a hollow fiber microfiltration (0.2 µm) and an ultrafiltration (100 kDa, UFP-100-E-4A) membrane (Figure 1b). The first retentate was rich in suspended grain particles and macromolecules, allowing its application as a carbon source to generate biogas, bioethanol, and relevant biotechnological products. The second retentate contained carbohydrates that could be further processed with a 1 kDa ultrafiltration membrane for the recovery of soluble calcium that can be reused in the maize's nixtamalization process and a polyphenol-rich permeate for food, pharmaceutical, and cosmetic applications. cereal processing by-products (e.g., from corn, rice, and wheat). The advantage of combined membrane systems is the synergy resulting from their integration, with overall significant benefits in plant compactness, product quality, energy efficiency, and environmental impact [18]. For instance, Ferri et al. [60] used a combined treatment with proteases (Protamex, alcalase, and neutrase) and sequential membrane processes for this purpose (Figure 1a). Rice by-products were first subjected to enzymatic hydrolysis (ph = 7 and 60 °C for 2 h) before recovering the digestate with centrifugation and its sequential fractionation with one microfiltration (0.2 μm) and three polyethersulfone ultrafiltration (8,5, and 4 kDa, respectively) steps. The retentate with the highest peptide content was the one recovered with the 8 kDa membrane. Compared to polysulfone, the polyethersulfone membrane is less hydrophobic, allowing the rapid absorption of polar compounds (e.g., polyphenols) and the release of soluble dietary fiber in the permeate [40]. In another effort, Castro-Muñoz and Yáñez-Fernández [65] used a combination of a hollow fiber microfiltration (0.2 μm) and an ultrafiltration (100 kDa, UFP-100-E-4A) membrane (Figure 1b). The first retentate was rich in suspended grain particles and macromolecules, allowing its application as a carbon source to generate biogas, bioethanol, and relevant biotechnological products. The second retentate contained carbohydrates that could be further processed with a 1 kDa ultrafiltration membrane for the recovery of soluble calcium that can be reused in the maize's nixtamalization process and a polyphenol-rich permeate for food, pharmaceutical, and cosmetic applications.  The extraction of proteins from precipitates or concentrates are typically conducted using the alkaline solution (e.g., 0.1 M NaOH at 60 • C for 60 min for proteins extraction from BSG) and subsequent filtering (180 µm) of the supernatant, and acid precipitation with 2.0 M citric acid addition (pH modification to 4.0) [66], or by adding Na 2 HPO 4 in sodium dodecyl sulfate solution [67]. For example, Idris et al. [61] investigated the recovery of proteins from wheat bran. At first, the adherent endosperm was removed from wheat bran with brushing before being extracted using the alkaline treatment and isoelectric precipitation. This approach led to high nitrogen dispensability with minimum and maximum protein solubility pH values of 5.5 and 11.5, respectively.
Alkaline extraction has also be used for the recovery and fraction of arabinoxylans from cereal processing by-products since under these conditions, and the hydroxyl ions accelerate cellulose's swelling and intermolecular hydrogen bonds' disruption [68]. Alkaline also hydrolyses the carboxyl and acetyl groups of benzyl groups and uronic acids of lignin. The ester linkages between arabinose residues and the ferulic acid, leading to the solubilization of hemicellulose from cell insoluble matrix. For example, Mandalari et al. [63] reported that the sequential fractionation of BSG and wheat bran using KOH of increasing strength (0.5, 1.0, and 4.0 M) of growing strength and 50 mM Na 2 CO 3 allowed the extraction of arabinoxylans with a lower ratio of arabinose/xylose as a factor of alkali's power.
However, there are concerns about the detrimental effects of this technique on proteins' and arabinoxylans' functionality and nutritional value, restricting their applications as additives in food products [18]. For example, the protein extracted using alkaline is poorly soluble, while the increased alkali strength results in the recovery of arabinoxylans with lower molecular [69]. The foaming and emulsification properties and the solubility of the isolated proteins could be increased using enzymes, e.g., specific proteases that modify proteins smoothly without destroying any amino acids [70]. This process generates a mixture of peptides owing different molecular weights and lengths as well as physicochemical and functional properties compared to non-hydrolyzed proteins, e.g., protein hydrolysates extracted from BSG exhibit immunomodulatory properties that could help the body against inflammatory diseases [71]. In order to increase proteins' extractability, relevant biotechnological approaches have been suggesting the co-generation of other products, and the low yields of target peptides have limited their application [18].
For instance, Tang et al. [62] applied ultrasounds to recover proteins from barley spent grain, referring to a yield of 104 mg/g under the optimum conditions (solid-to-liquid ratio of 2 g/100 mL and intensity of 88.2 W/100 mL of extractant). This yield (145.6 mg/g) can be further enhanced by coupling ultrasonic-assisted extraction with ultrafiltration using a 30-kDa membrane. Phongthai et al. [81] investigated the extraction of proteins from rice bran, referring to a yield of~4.7% under the optimum conditions of 0.99 g/10 mL solid-to-liquid ratio, 76% sonication amplitude, and 18 min extraction time. Moreover, the proteins were hydrolyzed using three enzymatic preparations (Subtilisin A, Actinase E, and Neutrase 0.8 L), resulting in different hydrolysis degrees (e.g., 20%, 14%, and 6%, respectively). By increasing the hydrolysis degree, the foam capacity decreased, and the antioxidant capacity increased, resulting in peptides with different properties and potential applications. In another study, Phongthai et al. [64] used microwaves to extract proteins from rice bran. The optimum conditions were found to be a power of 1000 W, a water solid-to-liquid ratio of 0.89 g bran/10 mL, and an extraction time of 90 s. Microwaves are enhanced by 1.5 fold the protein yield of alkaline extraction. By applying a lower alcalase hydrolysis degree of 5%, the emulsification and foaming capacities increased, and the antioxidant capacity decreased.

The Biorefinery Concept
The modern bioeconomy requires not only the application of sustainable technologies, but also their integration in the biorefinery concept for the conversion of biomass to a range of biobased products for food, energy, textile, and other industrial applications [82]. Biobased products could replace products manufactured from fossil fuel feedstocks taking into account that renewable biomass prices have slowly and steadily decreased. Over the last years, the utilization of commonly used biomass materials (e.g., potato, wheat, sugar beet, and corn) as a source of glucose for the production of biofuels and biobased chemicals has been partially commercialized [18]. Besides, the feasibility of bioethanol production and appropriate processing steps from corn has been demonstrated by performing calculations that scale-up benchtop or pilot plan operations [83].
The first generation of bioethanol production is based on transforming the whole grain kernel by grinding and mixing with enzymes and water to progress the degradation of starch [84]. At this process, kernel's components (e.g., proteins, dietary fiber, and germ) are not fermented, but are concentrated into a by-product, the so-called distillers' dried grains with solubles (DDGS). It is estimated that 100 kg of grains can generate~40 L of ethanol,~32 kg of DDGS, and~3 kg of CO 2 [85]. Another sector contributing to the increasing amounts of DDGS is whiskey's distilleries that produce potable ethanol using blended grains of barley, rye, wheat, and maize. This substrate is more attractive compared to the DDGS derived from distilleries, as it possesses a food-grade nature. Thus, the corresponding recovered products could be used directly in food formulations [86].
DDGS is a heterogeneous material whose rich composition in protein, lipids, carbohydrates, and other valuable ingredients varies depending on the initially blended cereals [87,88]. DDGS is mainly used as livestock feed, but its composition allows considering it as a substrate that can generate numerous high added-value biobased products within the biorefinery concept [18]. To this line, two approaches could be used. In the first case, the grains are subjected to dry milling or the quick germ technique that starts by soaking the whole grain in water up to 3 h at 60 • C. Thereafter, components such as germ meal, oil, and arabinoxylans are recovered from the residual grains before the starchy endosperm fermentation, saccharification, and bioethanol production [89,90]. In the second case, the DDGS and its intermediate products are directly fractionated into valuable components such as xylooligosaccharides, protein, phenolic acids, oils, and phytosterols [91]. The selected process should not be affected by the variability and the composition of the feedstock and should be easily incorporated into existing production processes. From an economic perspective, the thermal processing of intermediate products is the most expensive process of bioethanol production [83].
The valorization of wheat bran also fits well with the biorefinery concept. Wheat bran includes approximately 15% of wheat kernel constituents and is rich in non-starch carbohydrates (55-60%, dry matter based), starch (14-25%), and protein (13-18%). Other minor constituents include fat (3-4%), minerals (3-8%), and other components like lignans, flavonoids, phenolic compounds, polyols, amino acids, and organic acids [92,93]. The non-starch carbohydrates of the wheat bran are mainly composed of soluble dietary fiber and particularly arabinoxylan (52-70%) [48]. In order to convert wheat bran to biobased products, several pretreatment methods are required such as acid hydrolysis, solubilization of lignin and hemicelluloses with organic solvents, enzymatic depolymerization of arabinoxylans, wet alkaline oxidation, and steam explosion, among others [94,95]. After the hydrolysis of cellulose and hemicelluloses, the released sugars are converted to bioethanol, glycerol, butanol, organic acids (e.g., levulinic, acetic, and formic and), and other products through submerged fermentation [96]. The products can be further transformed into other biobased products, e.g., polyethylene (used in packaging) by dehydrating ethanol to ethylene, polymerizing it, or polybutylene succinate produced by the esterification of succinic acid and butane-1,4-diol [95].

Applications
As shown above, primary and secondary by-products of cereals processing contain several high added-value compounds that could be converted into different commodities for different sectors. Besides, numerous applications have today been commercialized, e.g., Promitor Soluble Fiber from Tate & Lyle, Stabilized Rice Bran from NutraBio, Gama Oryzanol from Swanson Health Products, Oat Fiber Plus Tablets from Now Foods, and Life Extension NK Cell Activator from Swanson Health Products [18].
Cereal processing by-products are rich in nutrients that could replace the typical carbon sources in media preparations used in industrial enzyme production and microbial processes [97]. For instance, they could be utilized for the production of microbial enzymes (e.g., xylanase, cellulases, proteases, and amylases) and single-cell protein at low cost and high production yields [98][99][100][101]. The industrial production of enzymes is based on the cultivation of certain bacterial and fungal strains (mainly Aspergillus oryzea, Aspergillus niger, and Bacillus subtilis) that produce them during their metabolism. For instance, Gomathi et al. [102] investigated carboxymethyl cellulose production by Aspergillus flavus using submerged fermentation and wheat bran as a substrate. Hashemi et al. [103] studied the production of α-amylase using Bacillus sp. KR-8104 as a starting culture in a submerged fermentation system. The addition of BSG (by 5%, w/v) and removing dextrin from the culture medium enhanced α-amylase production by five times, making it more economically sustainable. Solid-state fermentation using different agricultural by-products (e.g., rice husk, gram bran, maize bran, wheat bran, and straw) has also been tested to produce xylanase by Bacillus subtilis ASH. Among the different substrates, wheat bran provided the highest xylanase yields because it contained high amounts of carbohydrates and proteins [104,105]. In a similar effort, Tanasković et al. [106] proposed the utilization of Bacillus sp. TMF-2 for the solid-state fermentation of wheat bread, as this strain triplicated the soluble phenolic content of wheat bran, accelerating a significant increase of antioxidant capacity and free radical scavenging activity.
The secondary cereal processing by-products (e.g., BSG and DDGS) have been investigated in different food applications, e.g., to fortify bakery products, snacks, cookies, and flavors. To this line, several pretreatment methodologies (e.g., chemical, biochemical, hydrothermal, enzymatic, or physical like pre-soaking, milling, and extrusion) have been applied for the retardation of rancidity and the optimization of end-products' quality. Moreover, safety assays should be conducted together with studies dealing with the rheological, functional, and sensory properties of the developed food products before their usage in commercialized commodities.
The most popular application of cereals' bran is their usage in baked products in order to increase their content in dietary fiber. The percentage of cereals' bran in bread may vary according to the application and the needed properties, e.g., how much we want to increase bread's fiber content, how much we want to decrease its glycemic index, and what kind of sensory properties are needed. The addition of bran in baked products may affect the quality negatively and sensory characteristics of baked products, e.g., causing porosity and elasticity changes, increasing bitterness due to the contained phenolic compounds (pinoresinol and syringic acid), reducing bread's volume, or reducing nutrients absorption due to their increased content in phytic acid [18,107,108]. Likewise, cereals' bran contains high amounts of oil that are susceptible to oxidation, causing off-flavors and in the baked products. Therefore, oil's removal from bran is necessary before application.
For example, the incorporation of wheat bran in a high percentage of 15-20% has been referred to as increased water immobilization during dough making and decreases its crumb textural quality and volume [109]. Coda et al. [110] referred that a 160 µm-particle size of wheat bran resulted in high volume bread. Le Bleis et al. [111] prepared a French bread dough by fortifying it with 20% wheat bran of coarse (1.8 mm) and fine (18 µm) particle sizes and noted many differences in bread's properties, e.g., loss of dough stability and increase of its viscosity, bread density and porosity during proofing, firmer crumb and stiffer crust as well as a decrease in mixer's mechanical energy. Boita et al. [112] noted that the increased incorporation of wheat bran in the flour (from 25% to 100%) increased water absorption. On the other hand, it decreased stability, dough development time, extensibility, and viscosity. These results were attributed to the thinned and weakened gluten network formed by incorporating wheat bran in the dough. In order to reduce the negative impact of bran on the functional properties of the dough and the increased bioavailability of the minerals, Sanz Penella et al. [113] combined the addition of wheat bran with phytase and α-amylase.
Cereal processing by-products have also been used for the development of extruded snacks. For instance, Nascimento et al. [114] combined BSG (up to 30% in the dough) with rice flour in order to prepare extruded snacks with desirable characteristics compared to those prepared only with rice flour. Stojceska et al. [115] used a similar percentage (30%) of BSG to develop extruded ready-to-eat snacks with increased dietary fiber content (at least by 10%). On the other hand, the incorporation of wheat bran in extruded food products faces similar problems with bread application in spite of sensory characteristics and texture, e.g., by increasing the amount of bran in the mixture, the hardness and density of the products grew, while the crispness and the expansion volume are decreased [116,117]. It is thus vital to optimize the amount of added wheat bran.
Cereals bran has also been suggested to fortify spaghetti. However, their increased incorporation is known to affect the cooking quality of pasta negatively, e.g., higher cooking losses and decreased water absorption [118]. Aravind et al. [119] prepared spaghetti from durum semolina substituted with different amounts (10-30%) of wheat bran and referred that the increased inclusion of bran led to higher cooking losses, decreased firmness, and greater surface roughness. On the other hand, the dietary fiber and antioxidants content was increased compared to the control pasta sample, while starch digestibility was not affected. On the other hand, Padalino et al. [120] investigated gluten-free spaghetti development using oat bran-rich (22%) in β-glucans and maize flour. According to their findings, the addition of the noted hydrocolloids improved the firmness and elasticity of the pasta, leading to low bulkiness and adhesiveness.
Finally, other cereal processing by-products such as rice distilling lees have been proposed as a potential substrate for the production of seasonings and flavors or to fortify cookies [121]. Ertaş [122] applied different processing techniques such as autoclave, microwave, and oven stabilization to improve the nutritional content of cookies. According to the results, the fortification with 30% microwave-treated bran provided cookies with high mineral content, while the autoclave-treated cookies showed the highest loss of phytic acid. The best sensory properties (e.g., appearance, flavor, color, and taste) were obtained for the cookies fortified with 10% of oven-treated bran.

Conclusions
The conventional utilization of cereal processing by-products for livestock feed and composting is a low added-value solution for the cereal sector that seeks more sustainable solutions within the stressing post-pandemic era and bioeconomy frame of our times. On the other hand, these by-products have considerable potential to be used as substrates for the production of different products for biotechnology, food, and pharmaceutical applications. The most sustainable strategy for the development of such supply chains is the biorefinery concept that is ideal for the valorization of side streams and the integration of recovery processes. Most of the research and market implementation studies conducted in the field deal with the valorization of wheat bran by-products and DDGS. After the recovery of valuable compounds such as proteins and arabinoxylans, a further valorization of the residual materials for the production of biofuels and other biobased products could increase the energy and material efficiency as well as the profit obtained from the downstream processing of the cereal processing by-products. To this line, innovative recovery approaches, non-thermal technologies, and more integral strategies should be further investigated together with tailor-made applications in foods and other sectors. Several compounds extracted from cereals, bran, and other waste streams have been referred to possess beneficial health claims for consumers, but more in vivo, human clinical trials regarding their digestibility, intake, and absorption in the body should be performed to validate these claims. The applications of cereal processing by-products or the corresponding recovered high added-value compounds in foods should also comply with the safety and quality regulations for human consumption, and thereby more efforts are needed in this direction. The most popular application is the fortification of bakery and other food products with cereals bran that brings not only health benefits and improved functional properties, but also problems in bread's quality and organoleptic character. Therefore, more insights and more in-depth investigations are needed to address these issues prior to the commercialization of respective applications.


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

text:
Multi-Trait Genomic Prediction Models Enhance the Predictive Ability of Grain Trace Elements in Rice

Multi-trait (MT) genomic prediction models enable breeders to save phenotyping resources and increase the prediction accuracy of unobserved target traits by exploiting available information from non-target or auxiliary traits. Our study evaluated different MT models using 250 rice accessions from Asian countries genotyped and phenotyped for grain content of zinc (Zn), iron (Fe), copper (Cu), manganese (Mn), and cadmium (Cd). The predictive performance of MT models compared to a traditional single trait (ST) model was assessed by 1) applying different cross-validation strategies (CV1, CV2, and CV3) inferring varied phenotyping patterns and budgets; 2) accounting for local epistatic effects along with the main additive effect in MT models; and 3) using a selective marker panel composed of trait-associated SNPs in MT models. MT models were not statistically significantly (p < 0.05) superior to ST model under CV1, where no phenotypic information was available for the accessions in the test set. After including phenotypes from auxiliary traits in both training and test sets (MT-CV2) or simply in the test set (MT-CV3), MT models significantly (p < 0.05) outperformed ST model for all the traits. The highest increases in the predictive ability of MT models relative to ST models were 11.1% (Mn), 11.5 (Cd), 33.3% (Fe), 95.2% (Cu) and 126% (Zn). Accounting for the local epistatic effects using a haplotype-based model further improved the predictive ability of MT models by 4.6% (Cu), 3.8% (Zn), and 3.5% (Cd) relative to MT models with only additive effects. The predictive ability of the haplotype-based model was not improved after optimizing the marker panel by only considering the markers associated with the traits. This study first assessed the local epistatic effects and marker optimization strategies in the MT genomic prediction framework and then illustrated the power of the MT model in predicting trace element traits in rice for the effective use of genetic resources to improve the nutritional quality of rice grain.

INTRODUCTION
Over half of the world's population relies on rice as a staple crop (Bandumula, 2018). Growing and consuming rice has relative merits, as rice is the major dietary source for both toxic and essential trace elements (Yang et al., 2018). For instance, Cd is a potent environmental and human health toxicant (Arao and Ae, 2003;Uraguchi and Fujiwara, 2012;Lien et al., 2021) transported into rice grain via the same channels with other trace elements Zn, Fe, Cu, and Mn (Sasaki et al., 2012;Hao et al., 2018;Han et al., 2021) of essential nutritional and physiological functions to plants, animal and humans species (Miller, 1970;Olivares and Uauy, 1996;White and Broadley, 2009;Aschner and Erikson, 2017;Gao and Xiong, 2018).
Recent advancements in genomic research and the increasing number of germplasm resources in gene banks offer a great opportunity to develop safe and nutritious rice varieties costeffectively. The trait's heritability indicates the potential that a given trait can be genetically improved. Previously, broad sense heritability of grain Zn, Fe, Cu, Mn, and Cd was found to be low (0-0.3), moderate (0.4-0.6) to high (0.6 or higher) (Norton et al., 2010;Pinson et al., 2015;Naik et al., 2020), indicating the practical possibility to improve these traits via selective breeding methods. Furthermore, many molecular genetic studies have identified numerous quantitative trait loci (QTLs) responsible for trace element uptake, transport, and accumulation into different rice tissues through genome-wide association studies (GWAS) or QTL mapping (Lu et al., 2008;Garcia-Oliveira et al., 2009;Ueno et al., 2010;Du et al., 2013;Meng et al., 2017;Swamy et al., 2018;Yang et al., 2018;Descalsota-Empleo et al., 2019;Frouin et al., 2019;Liu J. et al., 2021). As a genomics-enabled breeding approach, marker-assisted selection (MAS) is useful to improve trace element traits when genes/QTLs with large additive genetic effects exist (Wu et al., 2020). However, prominent nonadditive gene action has also been reported for trace element traits, making MAS-based strategies ineffective (Sharma V. et al., 2021). In addition, MAS-based breeding methods are practically ineffective at simultaneously exploiting information from multiple genes (Spindel et al., 2016) or traits (Van Der Straeten et al., 2020).
In contrast, genomic selection (GS) approaches make use of total genome-wide markers with either large additive effects or minor effects to derive the genomic estimated genetic values of genotypes (Meuwissen et al., 2001), which overcomes the constraints of MAS-based methods (de los Campos et al., 2009). Also, GS models can be modified to a multi-trait (MT) form to exploit available information from multiple traits simultaneously. The MT models used in GS heavily rely on genetic correlation between traits (Henderson and Quaas, 1976). This correlation possibly results from the pleiotropic effect (multiple traits controlled by the same QTL) or linkage disequilibrium (LD) between genes (Falconer, 1996). Exploiting multi-trait information in GS has been awarded with an increase in prediction accuracy ranging from 24% to 105% relative to single trait (ST) models Sun et al., 2017;Arojju et al., 2020). Besides gains in prediction accuracy, integrating MT models with appropriate cross-validation (CV) schemes compensated for the negative effect of small population size without affecting the prediction accuracy, enabling breeders to minimize phenotyping budgets (Lado et al., 2018;Arojju et al., 2020). The benefits of MT models under various CV schemes are yet to be studied in diverse rice collections. Nevertheless, MT models have shown their potential in predicting complex traits in rice, such as grain arsenic content (Ahmadi et al., 2021), grain yield , and root index architecture (Sharma S. et al., 2021).
Most of the MT genomic prediction studies discussed above only modeled the additive genetic effects. Non-additive effects are also essential components of the genetic effect and can benefit the predictive ability of MT models if accommodated (dos Santos et al., 2016;Lyra et al., 2017). However, nonadditive effects such as dominance or global epistatic effects may not be conserved during breeding due to chromosomal recombination events (Falconer, 1996;Breseghello and Sorrells, 2006;He et al., 2017). In contrast, the local epistasis that spans short segments of chromosomes can be preserved over generations (Akdemir and Jannink, 2015), as adjacent loci normally hold a strong LD (Ardlie et al., 2002). Earlier GS studies with ST models illustrated that accounting for local epistatic effects along with the main additive model increased the prediction accuracy of agronomic traits in wheat accessions (Akdemir and Jannink, 2015;Akdemir et al., 2017;He et al., 2017;Jiang et al., 2018;He et al., 2019). However, the benefits of modeling local epistasis effects in MT models remain unknown in crop or animal species.
Genomic prediction models can be extended to incorporate markers associated with causal QTLs, such as trait-associated SNPs (TA-SNPs), bridging the gap between biology and mechanistic GS models using uninformative genome-wide markers. Also, genomic prediction with markers derived from functional QTL is less reliant on LD patterns shared by training and target populations, possibly allowing robust prediction, especially across unrelated populations where LD decays more rapidly (Snelling et al., 2013). Simulation and empirical studies have shown that accounting for known QTLs improves the performance of genomic prediction models compared to models using uninformative genomewide markers (Bernardo, 2014;Owens et al., 2014). Alemu et al. (2021); Zhou et al. (2021) reported a two-to four-fold gain in prediction accuracy using GS + de novo GWAS (Spindel et al., 2016), in which the most significant TA-SNPs from a GWAS conducted on the training population are fitted as fixed effects in the model along with the polygenic background. Other groups (Bhandari et al., 2019;Ahmadi et al., 2021) also reported gains in prediction accuracy ranging from 16% to 32% by exploiting GWAS-derived TA-SNPs using trait-specific genomic relationship matrices (Zhang et al., 2014) in which markers with stronger association signals are assigned higher weights than markers with weaker associations. However, the application of the above methods has not always been beneficial (Veerkamp et al., 2016;Rice and Lipka, 2019) and has been shown to depend on the genetic architecture of the traits of interest, trait heritability, the number of underlying causal mutations, and their effect sizes (Huang and Mackay, 2016). In addition, the use of TA-SNPs in genomic prediction has been scarcely investigated in models accounting for the non-additive effects. The potential of GWAS-derived TA-SNPs on the predictive ability of MT models accounting for the local epistatic effects in diverse rice populations is yet to be demonstrated and worth inspecting.
There is a great scope for applying MT models to evaluate trace elements in large germplasm collections such as those archived in gene banks. Earlier studies using ST models showed that GS is a robust and cost-efficient tool to predict the genetic merit of individuals in large germplasm collections for various agronomic traits, such as grain yield in rice (Tanaka et al., 2021), biomass yield in sorghum (Yu et al., 2016), oil, protein, and yield in soybean (Jarquin et al., 2016), total root length in maize (Pace et al., 2015), and days to head and days to maturity in wheat . However, the application of MT or even ST models to predict the concentrations of trace elements in food crops is still limited to a few studies involving arsenic Ahmadi et al., 2021), Mn (Leplat et al., 2016), and Zn (Guo et al., 2020) in rice, barley, and maize grain, respectively. Therefore, the overall goal of the present study is to compare the robustness of ST and MT models in predicting concentrations of four essential trace elements, Fe, Zn, Cu, and Mn, and one toxic metal, Cd, in rice grain. Different CV schemes, implying varied phenotyping patterns and costs, were examined in our study to seek the most efficient phenotyping strategy when multiple traits are planned to be measured. In addition, we investigated whether incorporating local epistatic effects and using a selective marker panel of TA-SNPs derived from GWAS into MT models could further enhance the predictive ability of MT models.

Rice Materials
Our study used 250 rice accessions, including indica and japonica ecotype accessions from Asian countries (Supplementary Table  S1). Accessions from China are mainly landrace indica varieties mostly cultivated on Cd-polluted soils in Guangdong province, China (Long et al., 2014).

Plant Cultivation and Quantification of Trace Elements in Rice Grain
The procedures followed for growing the 250 accessions and determining concentrations of trace elements in rice grain were as previously described by Liu S. et al. (2021). Briefly, seeds from the 250 accessions were first cultivated in pots filled with soil collected from the experimental station of the Agricultural Genomics Institute at Shenzhen, China. Next, germinated seeds were selected and cultivated in seedling trays for 4 weeks. Healthy seedlings were then transferred into pots containing soil amended with an initial concentration of Cd of 0.5 mg kg −1 . Finally, all the seedlings were planted in an augmented randomized complete block design with two replicates of 25 accessions from 20 July 2019 to 2 October 2019. Our study was limited to a single environment. Multi-environment data would be essential for understanding the environmental correlations and their stability and genotype effects by environment interactions (GxE). To determine grain concentrations of Zn, Fe, Cu, Mn, and Cd, grain samples were first peeled and dried at 65°C for 3 days. The dried samples were then crushed, wet-digested in concentrated nitric acid (HNO 3 ) at 120°C for 30 min, and further digested with perchloric acid (HClO 4 ) at 180°C until the samples became transparent. The samples were then diluted with ultrapure water. Finally, the grain concentration of each trace element was determined using the Inductively Coupled Plasma Mass Spectrometry (ICP-MS) machine (Houk et al., 1980).

Genotyping
The 250 accessions were genotyped following the resequencing and variants-calling procedures of the rice 3K project as reported by Wang et al. (2018). The following steps were implemented for all the genotypes to merge the variants-calling: First, raw reads were aligned to the R498 reference genome (Du et al., 2017) using the program bwamem alignment software (Li, 2013). Next, the PCR duplicates were identified with Picard software, version 2.9.0 (http:// broadinstitute.github.io/picard/), and discarded. Following that, the GATK HaplotypeCaller engine (McKenna et al., 2010), with the option "-ERCGVCF," was used to call genotypes at each site. The resulting genomic variants called format (gVCF) for each genotype were combined using the GATK Genotype GVCFs engine. Next, the GATK hard filter pipeline was used to individually call SNP and INDEL variants from the population variant file. All the variants within 5 bp of an INDEL were discarded. A variant was confirmed if at least one genotype supported it with a QUAL parameter greater than 30. After that, VCF tools indicated (Danecek et al., 2011) sites for which genotypes were not called in at least 20% of the used genotypes. The above procedures yielded 30,089,814 biallelic SNPs for the 250 genotypes. SNP quality control steps were implemented using PLINK software (Chang et al., 2015) with standards that remove SNPs with 1) minor allele frequency lower than 0.5, 2) call rate less than 0.9, and 3) pairwise LD (r 2 ) greater than 0.01. Finally, 36,171 SNPs were available for the 250 accessions.

Estimation of Genomic Heritability and Traits Genomic Correlation
The mixed linear model was used to estimate genomic heritability as follows: Where y is the vector of concentration of trace element under consideration, 1 n is the vector of ones, n is the number of genotyped cultivars, μ is the intercept, g is the vector of the genetic effects of accessions, e is the residual vector. g and e are assumed as random effects, respectively, following g~N(0, Gσ 2 g ) and e~N(0, Iσ 2 e ) where G is genomic relationship matrix estimated following (Yang et al., 2010) and I is identity matrix. The genomic heritability was estimated as Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 883853 Where σ 2 g is the additive genetic variance component and σ 2 e is the residual variance component. GCTA software (Yang et al., 2011) was used to compute the genomic relationship matrix and genomic heritability.
The genomic correlation between traits was estimated using the formula: cor covi,j σ i σ j where cov i,j is the genetic covariance between i th trait and j th trait, σ i and σ j are the square root of the genetic variance of i th trait and j th trait. The genetic covariances and variances were estimated using the R package MTM (De los Campos and Grüneberg, 2016).

Genetic Diversity
A hierarchical cluster analysis based on the Euclidean's distance matrix computed with the SNP genomic profiles was performed to inspect the genetic diversity among the 250 genotypes. In addition, a heat plot based on the cluster analysis was drawn to visualize the genetic dissimilarities.

Genomic Prediction Approaches
The genomic prediction models used were the ST and MT models. The ST model only captured additive genetic effects, while the MT models accommodated both additive and local epistatic effects. ST model was the commonly used genomic best linear unbiased prediction (GBLUP) model and same as the mixed model estimating genomic heritability: y 1 n μ + g + e, where y, 1 n , μ and ε were exactly as afore denoted, g is the vector of additive genetic effect in genotypebased model or additive plus local epistatic effects in the haplotype-based model. In the genotype-based model, we assumed g~N(0, Gσ 2 a ), where G is an n × n-dimensional additive genomic relationship matrix, and σ 2 a is the additive genetic variance component. In the haplotype-based model, we assumed g~N(0, Hσ 2 h ) , where H is an n × n-dimensional haplotypic relationship matrix derived from the haplotypic profile matrix with values 0, 1, and 2 indicating the number of copies for a specific haplotypic allele (Jiang et al., 2018;He et al., 2019). To obtain the haplotypic profile matrix, the genotypic data of SNPs were phased using software SHAPEIT (Delaneau et al., 2012) with default augment settings. The phased genotypic data was recoded to haplotypic profiles using a fixed-length haplotype of 2, 3, 4, or 5 SNPs. G and H were established using software GCTA (Yang et al., 2010;Yang et al., 2011) based on the genotypic and haplotypic data, respectively.
For MT models, we used two approaches considering no correlation between traits; the Bayesian multi-output regressor stacking (BMORS) proposed by Montesinos-López et al. (2019); the MT-GBLUP; and two methods accommodating correlation between traits; the factorial analytic (FA) model and the unstructured (UN) model. For genotype-based approaches, the MT-GBLUP model was formulated as y μ + u + ε, where , σ 2 g j and σ 2 ε j are respectively the genetic and residual variance of jth trait, and ⊗ denotes Kronecker product of matrices. BMORS was a two-stage process. The first stage is the same as the MT-GBLUP, but instead of directly using the GBLUP predicted values as the final output, BMORS implemented a second stage that integrated the GBLUP predicted values of each trait in the first stage and fitted a ridge regression model. In this way, the prediction of a single trait could be corrected by the predictions of other traits in the first stage using the second stage model (Spyromitros-Xioufis et al., 2012;Spyromitros-Xioufis et al., 2016;Montesinos-López et al., 2019). The FA model was also based on the formula of the MT-GBLUP model but assuming a covariance structure between traits, that is where k is the number of latent factors, Λ is a j × k dimensional matrix containing trait loadings, Π is a j × j diagonal matrix (Burgueño et al., 2012). Theoretically, the FA model requires at least three traits to be simultaneously included in the model. The UN model (Burgueño et al., 2012;Cuevas et al., 2017) tried to estimate all variances and covariances in Ψ u , i.e., σ 2 gj and cov gigj , i, j ∈ {1, . . ., m}, which may cause convergence problems when a large number of traits are considered. The haplotype-based approach was only implemented in the MT-UN model by replacing the relationship matrix G by H .
The ST-GBLUP model was implemented in R (R Core Team, 2016) using the BGLR package (De los Campos and Pérez-Rodríguez, 2015). The MT-GBLUP, FA, and UN approaches were realized using the R package MTM (De los Campos and Grüneberg 2016). BMORS was fitted using the R package BMTME (Montesinos-López et al., 2019). The number of iterations of all models was set to 20,000, and the first 12,000 were discarded as burn-in.

Cross-Validation Schemes and Evaluation of Genomic Prediction Accuracy
Four different CV schemes, referring to those reported by Lado et al. (2018) and Arojju et al. (2020) were used in our study ( Table 1). CV1 was applied to both ST-GBLUP and MT models, referring to a scenario where the target trait was predicted without the support of auxiliary traits (ST-CV1) or with auxiliary traits only available in the training set (MT-CV1). CV2 and CV3 were only assessed for MT models. Under CV2 scheme, genotypes in both training and test sets had phenotypic data for all the auxiliary traits. Under CV3, phenotypes of the auxiliary traits were only available in the test set.
To assess the genomic prediction accuracy across the above CV schemes, the entire population of 250 genotypes was randomly divided into five equal-sized folds. Four folds collectively constituted the training set, and the remaining fold was the test set. Stochastic partitioning of training and test sets was repeated 20 times, yielding one hundred times (5 folds × 20 replicates) calibrations and predictions. The genomic prediction accuracy of the target trait was estimated using the Pearson correlation coefficient between the genomic predicted genetic values and the observed phenotypic values of 250 accessions when incorporated in the five test sets of each repeat of CV. The Student t-test was used to test the statistical difference in genomic prediction accuracies among the prediction models.

Selective SNP Marker Panel
To investigate whether the predictive ability of MT genomic prediction on rice grain trace elements concentration could be boosted by optimizing the SNP marker panel, we applied a (GWAS) to identify the trait-associated SNPs (TA-SNPs) and establish the selective marker panel. The CV scenarios in which the MT-UN haplotype-based models disregarding the length of haplotypes (two to five SNPs) constantly showed statistically significantly (p < 0.05, t-test) higher prediction accuracies than their genotype-based counterparts and the ST-GBLUP model were used to validate the efficacy of using the TA-SNPs to train the genomic prediction models. In more detail, GWAS using the total SNP marker panel was performed in the training set of the designated CV scenarios. First, the TA-SNPs with p values less than 0.01 were recorded. Then each chromosome was divided into bins spanning 300 kb (the bin size is decided by the LD decay, with the physical distance between pairs of SNPs based on the total population). Finally, the most trait-associated SNPs with the lowest p-value in each bin was picked together with the TA-SNPs (p < 0.01) and recorded to constitute the selective SNP marker panel of each repeat of CV. The MT-UN genotype-based model was implemented using the genotypic data of the selective SNPs. Contrastingly, the adjacent selective SNPs located within 300 kb (highly possible as the position of the selective SNP from each bin is unfixed) were combined to compile the haplotypes using the phased genotypic data as the LD decay implied a non-negligible LD among them. The remaining SNPs without close neighbors within 300 kb were maintained, and their genotypic profiles were used. Therefore, the MT-UN haplotype-based model took advantage of both haplotypes and genotypes. The GWAS was implemented using GCTA software (Yang et al., 2011;Yang et al., 2014). The additive genomic relationship matrix was exclusively used in the GWAS model to account for the relatedness between accessions.

Linkage Disequilibrium Decay, Kinship, and Population Structure
The LD decay distance between all SNP markers for the 250 accessions was~250-300 kb when the cut-off value (r 2 ) was set at 0.1, assuming non-negligible SNP pairwise correlation (r = 0.3) ( Figure 1A). The kinship between accessions was determined based on pairwise Euclidean distances. Pairwise distances among accessions ranging from 0 to 0.2 accounted for less than 5% of all the pairwise distances. Pairwise distances from 0.6 to 0.8 were the most frequent and accounted for 12%-23% of all the pairwise distances ( Figure 1B). Further, no genetically structured sub-populations were observed among the 250 varieties used in this study. However, several families were detected (Supplementary Figure S1).

Distribution of Phenotypes, Genomic Heritability, and Genetic Correlation
The distributions of phenotypes (adjusted phenotypic means) based on the 250 accessions varied among the five traits studied. The distribution of Zn was almost symmetrical. The skewness was high and negative for Cu, and moderate and negative for Fe, Cd, and Mn (Supplementary Figure S2). The genomic heritability for all studied traits ranged from low (Zn: 0.14 and Cu: 0.21) to medium (Mn: 0.35) and high (Fe: 0.5 and Cd: 0.62) (     (Figure 2). Also, under the CV1 scheme, prediction accuracies of MT models were not statistically significantly (p < 0.05) superior to those of ST-GBLUP irrespective of the models and traits studied (Supplementary Tables  S2-S6).
As compared, when phenotypes of the auxiliary traits were made available in both training and test sets (MT-CV2) or merely in the test set (MT-CV3), the MT models, namely FA or UN, significantly (p < 0.05) outperformed the ST-GBLUP model (Figures 3, 4). For most of the studied traits, the highest performance of MT models was observed under the MT-CV2 scheme (Figures 3, 4). However, MT-GBLUP and BMORS MT models were not significantly (p < 0.05) superior to the ST-GBLUP model for all the CV schemes studied (Supplementary  Tables S7-S11).
We further compared scenarios where the prediction of the target traits was assisted with a single auxiliary trait or a combination of multiple auxiliary traits in MT models. Supporting the prediction of Zn with one of its correlated traits (Cu, Fe, Mn, or Cd) was sufficient to significantly (p < 0.05) increase the prediction accuracy MT-UN model relative to ST-GBLUP model in MT-CV2 ( Figure 3A). Cu was the best single auxiliary trait for predicting Zn. Incorporating observations from Cu in MT-UN model (under MT-CV2) significantly (p < 0.05) increased the prediction accuracy of Zn by 82.6% (0.23-0.42) relative to the ST-GBLUP model. However, the highest increase in prediction accuracy (126% or 0.23-0.52) of MT models was observed when observations from Mn, Fe, and Cu were combined as supporting traits for Zn under the MT-CV2 scheme ( Figure 3A). Under MT-CV3, the MT-UN model outperformed the ST-GBLUP model only after multiple auxiliary traits were used to support the prediction of Zn ( Figure 3A).
Similarly, compared to ST-GBLUP, the prediction accuracy of Cu by the MT-UN model significantly (p < 0.05) increased by  Figure 3C).
Furthermore, the prediction accuracy of Cd (the most heritable trait) with the MT-UN model was significantly (p < 0.05) improved by 3.8% (0.52-0.54) when Fe, its strongly correlated trait, was used as a single auxiliary trait under the MT-CV2 scheme ( Figure 4A). As observed for the other traits, 7.6% (0.52-0.56) and 11.5% (0.52-0.58) gains in prediction accuracy were attained after using combined information from multiple auxiliary traits (Mn, Fe, and Cu) in MT-UN and MT-FA models under MT-CV2, respectively ( Figure 4A). Similarly, under MT-CV3, MT models did not significantly outperform ST-GBLUP models in scenarios where a single auxiliary trait was used. However, when information from Fe and Zn or Fe and Cu was considered, an improvement of 5.7% (0.52-0.55) in the prediction accuracy of the MT-FA model over the ST-GBLUP model was observed ( Figure 4A).
Finally, when Mn was the target trait, MT-UN with a single auxiliary trait failed to improve its prediction accuracy under both MT-CV2 and MT-CV3 schemes ( Figure 4B). However, considering information from additional traits (Cu, Zn, and Cd), using the MT-FA model significantly (p < 0.05) improved the prediction accuracy of Mn up to 11.1% (0.36-0.40) over ST-GBLUP under MT-CV2 ( Figure 4B). On the other hand, under the MT-CV3 scheme, MT-UN or MT-FA did not significantly outperform ST-BLUP even after using multiple auxiliary traits.

Prediction Accuracy of Haplotype-Based Model Versus Genotype-Based Models
We further investigated the benefits of accommodating local epistatic effects on the prediction accuracy of MT models by using haplotypes instead of genotypes in the UN model. Comparing to the genotype-based UN model, the observed largest and significant (p < 0.05) increment of prediction accuracies using haplotype-based models was 3.8% for Zn (0.52-0.54), 4.6% for , and 3.5% (0.56-0.58) for Cd under MT-CV2. For Zn, the above improvement in prediction accuracy was achieved with a haplotype length of 3 SNPs, and when Mn, Fe, Cu, and Cd were collectively used as auxiliary traits ( Figure 5A). For Cu, the observed gains were realized with a haplotype length of 4 SNPs and when auxiliary traits Mn and Zn were used together ( Figure 5B). For Cd, the gains were from the MT-UN model with a haplotype length of 2 SNPs, and when Mn, Fe, Cu, and Zn were combined as the auxiliary traits ( Figure 5C). Under MT-CV3, the haplotype-based UN model was significantly (p < 0.05) superior to the genotype-based UN model by 12.5% (0.32-0.36) for Zn ( Figure 5A) and 6% (0.33-0.35) for Cu ( Figure 5B).

Prediction Accuracy of a Haplotype-Based Model Capitalizing on Trait-Associated SNPs
With this study, we sought to investigate whether the prediction accuracy of Zn, Cu, and Cd with the haplotypebased model can be improved by using the selective marker panel made by TA-SNPs derived from GWAS (Supplementary Table S12). We purposely selected Zn and Cu because for both traits the haplotype-based MT model performed superiorly for several scenarios irrespective of the lengths of haplotypes (2-5 SNPs) ( Figure 5). We also investigated Cd in addition to Zn and Cu since successive significant (p < 0.2) TA-SNPs were observed in GWAS based on the total population ( Supplementary Figures 4A-C). The non-negligible LD (r 2 ≥0.1) observed between the TA-SNPs, especially in Cu and Cd, underpinned the necessity of modelling local epistatic effects among TA-SNPs ( Figures 6A-C). The haplotype-based UN model accounting for TA-SNPs significantly (p < 0.05) outperformed their genotype-based counterparts; however, it was significantly (p < 0.05) inferior to the model using all genome-wide markers for all the traits and scenarios evaluated ( Figures 7A-C).

DISCUSSION
Quantifying trace element content in food crops is labor-and time-intensive. As a result, trace element traits have been the subject of few genomic prediction studies (Owens et al., 2014;Leplat et al., 2016;Frouin et al., 2019;Guo et al., 2020;Ahmadi et al., 2021) compared to agronomic or physiological traits. This study demonstrates how MT models with appropriate CV strategies can be useful in saving phenotyping resources for trace element traits in diverse rice collections without compromising the prediction accuracy. It also provides the first proof of concept in diverse rice for incorporating local epistatic effects and trait-associated SNPs into MT genomic prediction models.

Multi-Trait Models Improved the Prediction Accuracy of Trace Elements in Rice Grain
In this study, MT models did not significantly outperform ST-GBLUP under the CV1 scheme for all the scenarios evaluated (Supplementary Tables S2-S6). Earlier studies also reported insignificant differences in the prediction accuracies of MT-CV1 and ST-CV1 (Calus and Veerkamp, 2011;dos Santos et al., 2016;Bhatta et al., 2020), implying that MT models are not always robust over ST models, especially when information on auxiliary traits is only available in the training set and the unobserved accessions are predicted only based on genotypic data. In contrast, when phenotypes of the auxiliary traits were present in the training and test set (MT-CV2) or merely in the test set (MT-CV3), the prediction accuracy of MT models (MT-UN and MT-FA) for the unobserved target traits (Zn, Cu, Fe, and Cd) was significantly improved relative to ST-GBLUP (Figures 3, 4). Previous studies attributed the predictive performance of MT models to both higher heritability of the auxiliary trait and strong genetic correlation between the target and auxiliary traits (Sun et al., 2017;Fernandes et al., 2018).

Accounting for the Information From Multiple Auxiliary Traits Boosted the Predictive Ability of Multi-Trait Models
Using a single auxiliary trait in the MT-UN model significantly (p < 0.05) improved the prediction accuracy of target traits Zn, Cu, Fe, and Cd relative to the ST-GBLUP model (Figures 3, 4). When a strong genetic correlation exists between target and auxiliary trait, the prediction accuracy of MT models could still be improved under MT-CV2 or MT-CV3 regardless of trait heritability. For instance, supporting the prediction of Cu with Zn, its strongly correlated trait (cor Zn, Cu = 0.95) but with lower heritability (h 2 Zn = 0.14; h 2 Cu = 0.21) significantly improved the prediction accuracy of Cu with the MT-UN model ( Figure 3B). Also, supporting Cd with Fe, its strongly correlated trait (cor Cd, Fe = 0.95) but with lower heritability (h 2 Cd = 0.62; h 2 Fe = 0.50), improved the prediction accuracy of Cd with the MT-UN model ( Figure 4A). Arojju et al. (2020) also indicated that the genetic correlation was the main cause of the observed gain in prediction accuracy of MT models. The same study further showed that when a trait in strong genetic correlation with the target trait is used in the MT model, the predictive performance of the MT model was still superior to the ST model even after reducing the training population size by 50%.
Collectively accounting for phenotypes of multiple auxiliary traits further improved the predictive ability of the MT models compared to the MT models with a single auxiliary trait. For example, the highest increase in the prediction of Zn was 82.6% when a single auxiliary trait was used in MT-UN models. Yet, using multiple traits collectively in the same model improved the prediction of Zn by 126% compared to the ST-GBLUP model ( Figure 3A). Also, MT models with one auxiliary trait showed no benefit over the ST-GBLUP model when predicting Mn, with relative medium heritability and no strong genetic correlation with any other studied trait. However, when auxiliary traits were collectively used in the MT-FA model, significant improvements in the prediction accuracy over the ST-GBLUP model were observed ( Figure 4B). Multiple auxiliary traits would optimize MT models, though the assisting trait per se is neither strongly genetically correlated with the target trait nor highly heritable. Therefore, when no single auxiliary trait meets the criterion of heritability or genetic correlation, combining multiple auxiliary traits in the MT model could be an effective approach to enhance the predictive ability of MT models. These findings are concurrent with previous findings by Wang et al. (2017), indicating that the prediction accuracy of MT models was highest when eight different traits were used as auxiliary traits to predict grain yield in rice.

Modeling Local Epistatic Effects is Beneficial in Multi-Trait Models Irrespective of Using Total or Selective Marker Panel
Previous studies demonstrated that accounting for local epistatic effects besides the additive effect in genomic prediction could improve the prediction accuracy of ST models (Akdemir and Jannink, 2015;Akdemir et al., 2017;Jiang et al., 2018;He et al., 2019). Here, we are the first to attempt to model the local epistatic effects in the context of MT genomic prediction. Accounting for the local epistatic effects in haplotype-based MT models significantly improved the prediction accuracy of Zn, Cu, and Cd relative to genotype-based MT models, only capturing additive effects ( Figures 5A-C). Relative to ST-GLUP, the highest increase in prediction accuracy, 134.7% for Zn, was observed after incorporating the local epistatic effects into the MT-UN model ( Figure 5A). These findings imply that the potential of MT models can be maximized by accounting for local epistatic effects besides additive effects in the model. Using a selective marker panel based on approaches exploiting the trait biological and genetic background knowledge such as GWAS has been proven effective to improve the predictive ability of GS models (Owens et al., 2014;Wang et al., 2019). Our study did not show any improvements of prediction accuracy by using the TA-SNPs instead of all genome-wide SNPs for Cu, Zn, and Cd ( Figures 7A-C). These findings could be attributed to the complex genetic architectures of the trace elements we studied ( Supplementary Figures S4A-C). Our approach was slightly similar to previous methods using GWAS-derived TA-SNPs to construct the trait-associated matrix (Zhang et al., 2014;Ahmadi et al., 2021), except that we did not assign weights to haplotypeor genotype-based genomic relationship matrices. Though numerous studies reported improved gains from using the above strategy (Bhandari et al., 2019;Ahmadi et al., 2021), Veerkamp et al. (2016) showed that the proportion of total variance explained by the TA-SNPs combined in a GRM was considerably smaller than that explained by all variants in Holstein-Friesian cattle population. A potentially more promising way to use TA-SNPs would be to fit them as fixed effects in the GP model along with all other SNPs as random effects (Spindel et al., 2016). However, the latter approach is best suited for features with a few large-effect QTLs in a polygenic context Bian and Holland, 2017;Rice and Lipka, 2019). Therefore, the genetic architecture of the target traits must be studied before applying this strategy to a breeding program. The marked advantage of the haplotype-based UN model over their genotype-based counterparts using the TA-SNPs ( Figures 7A-C) substantiates the existence of local epistasis in trace element traits (Sharma V. et al., 2021) and the merit of modelling local epistatic effects in MT-GP program.

Factors Affecting the Observed Prediction Accuracies: Trait Heritability, Genetic Correlation, and Population Relatedness
Various factors affect the predictive ability of GP models used in GS Xu et al., 2021). In this study, the genomic heritabilities spanned a wide range from 0.14 to 0.62 (Table 2), which enabled the evaluation of the performance of MT models under contrasting levels of genomic heritability. The genomic heritability of Zn was the lowest (h 2 = 0.14, Table 2), which contradicts several previous studies that reported moderate to high heritability of Zn (Norton et al., 2010;Pinson et al., 2015;Naik et al., 2020). The poor heritability estimate of Zn in this study could be due to potential environmental effects. Unfortunately, our study does not include multi-environment trials and therefore does not provide insight into environmental factors and GxE interactions on genomic prediction of trace element traits.
Strong genetic correlations (Cor > 0.75) were observed between pairs of the studied trace element traits ( Table 2). This was expected due to their overlapping genetic and physiological mechanisms (Sasaki et al., 2012;Cu et al., 2020). For example, transporter gene families like zinc-iron permease (ZIP), natural resistance-associated macrophage proteins (NRAMPs), and heavy metal transporting, ATPases (HMAs) have been associated with uptake and translocation of several trace elements in plants (Fernández-Paz et al., 2021;Vanderschueren et al., 2021;Zhang et al., 2022). As a result, borrowing information from correlated traits overall improved the prediction accuracy of the MT models.
The success of GS is also highly dependent on the LD between markers and unknown causal variants. The genetically distant training and test sets would have different LD decay patterns and consequently impede the prediction (Snelling et al., 2013;Desta and Ortiz, 2014;Thistlethwaite et al., 2020). Such a problem is typically prominent in germplasm accessions, limiting the power of GP . This is also the case for the diverse rice population used in this study, as most accession pairs were distantly unrelated ( Figure 1B). As a result, the predictive ability for ST-GBLUP was poor (Figure 2), particularly for traits with the lowest estimated heritabilities, namely Zn and Cu ( Table 2). Adding related materials to the training population has been suggested to overcome the problem of low relatedness between training and test populations and improve the accuracy of genomic prediction (Arenas et al., 2021). Nevertheless, increasing relatedness will damage genetic gain in the long term because genetic variation will be limited or exhausted if related populations are overused (Jannink et al., 2010;Moeinizade et al., 2019). We show that MT models are powerful tools for predicting trace element traits in populations with diverse backgrounds. However, further studies with larger datasets are needed to elucidate the utility of different populations and marker optimization strategies in the context of MT genomic prediction. Expected prediction accuracy and phenotyping cost for different CV schemes and MT models. Green dots represent models which account only for additive effects. Black dots represent models considering both additive and local epistatic effects. The GBLUP model under ST-CV1 is economically advantageous because the main effort is just devoted to phenotyping one target trait in the training set. However, in terms of prediction accuracy, it is less robust than the UN and FA MT models under MT-CV2 or MT-CV3. Compared to ST-GBLUP, using MT models under MT-CV1 has no advantage in phenotypic resource-saving or prediction accuracy. In contrast, MT models implemented under MT-CV2 and MT-CV3 can improve prediction accuracy. However, high phenotyping efforts can be expected with MT-CV2, mainly when multiple auxiliary traits need to be phenotyped. MT-CV3 saves resources by only phenotyping the test set population (20% of the total population in our case). Accounting for local epistatic effects may further improve the predictive ability of MT models under MT-CV2 or MT-CV3.
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 883853 The Prospect of Multi-Trait Models for the

Evaluation of Genetic Resources in Gene Banks
Expediting genomic selection in gene banks to predict the genetic merit of the unobserved accessions would enable accurate identification of promising donor accessions without a comprehensive phenotypic test of all the accessions in the field (Pace et al., 2015;Yu et al., 2016;Crossa et al., 2017;Tanaka et al., 2021;He et al., 2022). In fact, as the traits of breeders' interest are extensive, the genetic resources archived in gene banks would be evaluated for several traits. MT genomic prediction is an effective method to realize this comprehensive evaluation. Our study tested several MT models under three different CV schemes, implying different phenotyping layouts and costs. We found that high prediction accuracy of MT models can be achieved under prediction schemes MT-CV2 and MT-CV3. The MT-CV2 scheme requires more budget for phenotyping auxiliary traits in both training and test sets. Therefore, breeders may kindly MT-CV2 if phenotypes for the auxiliary traits can be inexpensively obtained. Otherwise, MT-CV3 is more cost-effective as phenotypes for the auxiliary traits are only required for the test set (i.e., 20% of the entire population). Since using multiple auxiliary traits collectively in the MT model can improve prediction even if the individual auxiliary traits do not fully meet the heritability and genetic correlation conditions, an ideal situation would be to phenotype less expensive and more manageable traits (e.g., root system architecture, 100-grain weight, data to heading, etc.) to support the prediction of expensive target traits with the MT models. Besides, accounting for local epistatic effects in MT models would help to improve the predictive ability. The different scenarios studied here and their respective potentials in terms of prediction accuracy and phenotyping cost are illustrated in Figures 8A,B.
To breed safe and nutritious crop varieties, further studies using the genomic selection index (Habyarimana et al., 2020), for example, are desired to provide a comprehensive understanding of the strategies to optimize essential nutrients and toxic metals such as Cd in food crops.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
GY and SH conceived the project. WZ, JZ, and MB prepared the samples and collected the data. SH, JZ, and MB analyzed the data. WZ and MB drafted the manuscript. SH and GY revised the manuscript. All the authors have read, edited, and approved the current version of the manuscript.