The use of mathematical models in the mineral nutrition of plants allows to simulate the dynamics of the absorption of water and dissolved ions in response to different internal and external factors (Juárez-Maldonado et al., 2017). The models contribute to the quantitative understanding of the factors involved in the absorption, transport, and accumulation of mineral elements; additionally, they allow to explore different environmental or endogenous situations that modify the nutrition of the plant (Mankin and Fynn, 1996). Regarding Si modeling, (Sakurai et al., 2017) presented a dynamic model of Si absorption and transport for rice. The model was integrated considering the activity of different transporters and the distribution of Si through different nodes in the entire plant; the model was able to predict the dynamic behavior of silicon in the plant successfully. However, in the case of vegetables, there are no models that consider silicon, although there are published models that effectively simulate nutrition with other mineral elements (Juárez-Maldonado et al., 2014b, 2017; Ramírez-Pérez et al., 2018).

The aim of this manuscript was the integration of an estimated dynamic model that approximates the availability, absorption, transport, and accumulation of silicon in a tomato crop in soil and soilless conditions.

### BENEFITS OF SILICON IN HORTICULTURAL PLANTS

The ferns, horsetails, and grasses such as corn, wheat and rice, and sugar cane, are the plants that naturally accumulate more silicon (Liang et al., 2015). However, in the presence of adequate amounts of silicon in the form of Si(OH)4, all plants, including horticultural species such as tomato and cucumber, absorb it. Plants use silicon in a manner not yet well understood to stimulate the antioxidant metabolism, the processes of plant's hardening, defense, and adaptation to environmental factors. Depending on whether they are species that carry out silicification, Si(OH)4 is concentrated in polymeric form (amorphous hydrated silica) in different cellular and extracellular compartments, finally transforming it using a deposition process dependent on transpiration into insoluble biogenic silica (SiO2.nH2O) which forms structures called phytoliths or opal (Sangster et al., 2001; Katz, 2014; Exley, 2015). The biogenic silica is subsequently incorporated into the soil contributing with 1–3% of the total Si in the soil (Desplanques et al., 2006).

The silicon absorbed by the plants seems to be maintained under constant exchange between the soluble forms (Si(OH)4) and the insoluble fraction (polymeric silicic acid and biogenic silica) (Exley, 2015). Most of the Si deposited as biogenic silica remains as such throughout the life of the plant (Sangster et al., 2001). The soluble part is directly available to humans when they consume plant foods, while the insoluble fraction could perhaps be considered as an integral part of the fiber. The environmental factors and differences between plant species that modify the ratio Si soluble/Si insoluble, which ultimately determines the dietary utility of the product, have not been studied.

Considering that, (i) all plants seem to have the capacity to absorb silicon (Exley, 2015), (ii) and in view of the rise that has taken in recent years the production of vegetables using soilless production systems (Pignata et al., 2017), (iii) in addition to the fact that irrigation water and horticultural substrates provide little bioavailable silicon (Liang et al., 2015), then, it would be advisable to include silicon on a daily basis in fertilizer formulations used in those soils with low Si bioavailability as well as in soilless crops for the production of vegetables under protected conditions (Epstein, 1994, 2001).

The different groups of plants have different capacity to mobilize Si toward their various organs, but practically all absorb the silicon from the soil when it is available in the soil solution or the nutrient solution. The species with low mobilization capacity accumulate it in the roots and stems, while the species with high mobilization capacity accumulate it in stems, leaves, fruits, and seeds. Si appears to be absorbed in the form of Si(OH)4 by channels belonging to the aquaporins' group. Thus the rate of absorption and transport depends on the flow of water linked to transpiration (Exley, 2015; Sakurai et al., 2017).

The cereals are plants with a high capacity of silicification and therefore represent a significant amount of Si in the diet. However, the silicon contained in cereals will be encountered almost all in the form of insoluble biogenic silica (Sangster et al., 2001) which would be partially dissolved by the acids of the digestive system; on the other hand, in fruits and vegetables, due to their lower silicification capacity, it is expected that there will be more soluble silicon, which theoretically would be more available to be assimilated during intake. Considering the above, it is possible that the fruits of horticultural species such as tomato can be excellent sources of Si for the diet.

Additionally, it is known that in comparison with the dicotyledons, cereals contribute less Ca and Mg (White and Broadley, 2005). Therefore, a diet high in cereals that provides a significant amount of Si on average will contain less Ca and Mg than a mixed diet with base in cereals and dicots. On the other hand, the consumption of vegetables and fruits has grown considerably in recent decades among the human population and it is desirable that species such as tomatoes, eggplants, strawberries, cucumbers, avocados, melons, watermelons, carrots, onions, chilies, pumpkins, among others, contain a higher amount of silicon, considering the double benefit already mentioned of the crop higher tolerance to stress and the contribution of Si to human consumers.

In soilless crops, it is necessary to consider the contribution of Si in fertilizers since irrigation water does not provide enough, only from 5 to 24 mg L−1 Si (Liang et al., 2015). The lowest recommended concentration of Si in the nutrient solution for plants growing on substrates other than soil is 28 mg L−1 Si (Epstein, 1994), which can be achieved with 123 mg L−1 of Na2SiO3.

### EDAPHIC FACTORS THAT DETERMINE THE AVAILABILITY OF SILICON

Si is found in soil as an inert mineral in the form of quartz or aluminosilicates such as micas and feldspars. The weathering of these materials by rainwater, irrigation water, or by the acid metabolites of microorganisms and plant roots produce Si(OH)4 that under a balanced condition reaches a concentration of up to 1.8 mM (173 mg L−1 , equivalent to 50.4 mg L−1 Si). Above this level, reaching 2 mM (192.18 mg L−1 ), Si(OH)4 forms hydrated amorphous silica polymers containing Si unavailable for plants (Epstein, 2001; Liang et al., 2015).

The actual value of the concentration of Si(OH)4 in the soil solution is much lower than 1.8 mM, commonly found between 0.1 and 0.6 mM (9.61–57.66 mg L−1 ), but with such low values as 0.02 mM (1.92 mg L−1 ) in very eroded soils (Epstein, 2001; Liang et al., 2015).

The concentration of bioavailable Si in soil solution results from the release rate of Si(OH)4. Bioavailability is dependent on the silicon content of the soil minerals, organic matter, the temperature, the amount of precipitation and the acidity of the soil or soil pore water. The incorporation of Si in plants occurs at a rate dependent on the intensity of the transpiration (Exley, 2015), so that conditions of rapid growth can rapidly decrease the availability of Si in the soil solution (Epstein, 2001; Liang et al., 2015).

Soils of tropical areas where high precipitation occurs as well as calcareous and sandy soils of semi-arid and arid regions with low vegetation provide low quantities of Si to the soil solution, so it is recommended to use fertilizers with Si (Epstein, 1999). An affordable source of Si is siliceous sand that is offered in different granulometries, and that is used in quantities of between 500 and 4,000 kg ha−1 . On the other hand, Mollisol and Vertisol soils of the temperate and subtropical regions are soils that can provide an adequate amount of silicon (Epstein, 2001; Gérard et al., 2008). However, this has not been corroborated in regards of the actual availability of Si(OH)4 in soil pore water, since there is little-published information about concentrations, dynamic behavior, and association with other edaphic components of Si in the solution of the soil.

The temperature exerts a substantial impact on the soil solubilization rate of Si. However, the seasonal changes in temperature are significant as a determinant of the Si concentration in the soil solution only in the cold seasons of temperate zones, because the range of temperatures suitable for the growth of a crop is also adequate for the solubilization reactions of silicon in the soil (Sommer et al., 2006). Therefore, the temperature is not considered as a factor subject to management regarding the bioavailability of Si for crop plants. Possible exceptions would be crops in soil mulching and crops grown in greenhouse soil or tunnels. In both cases, soil or substrate temperatures are more stable, and on average higher than those of uncovered soil, in addition to water management more precise in time and quantity, so the bioavailability of Si is expected to be higher.

Another factor regulating the availability of Si(OH)4 is the pH of the soil pore water, that depends on the pH of the rainwater or irrigation water and is also modified by the respiratory activity and extrusion of organic acids by microorganisms and plant roots (Pérez-Labrada et al., 2016). In fact, the presence of Si induces the synthesis of citric acid in plants (Hernandez-Apaolaza, 2014). The pK1 of Si(OH)4 is 9.6, which indicates that its bioavailability in a nutrient solution is practically unaltered with pH values lower than 9. In the study of Gérard et al. (2008), there was little impact of pH on the bioavailability of Si in the soil solution, but the study conditions were developed under a very narrow range of pH variation. It will be necessary to collect data in different types of soil, or in soils subjected to treatments that modify its reaction or the pH of the soil pore water, to determine the effect of pH on the concentration of Si(OH)4.

Both a nutrient solution and the soil solution contain components that modify pH and interact with Si. With a pH> 7 that promotes the formation of Fe hydroxides, an adsorption process occurs that causes the polymerization of Si(OH)4. With pH< 6 Si(OH)4 begins to polymerize on surfaces with minerals containing Fe, while Al3+ would promote the stabilization of Si(OH)4 polymers, which would make Si unavailable for plants (Sommer et al., 2006). Considering this, it is possible that the availability of Si in the soil solution is higher with pH values between 6.0 and possibly 7.5 (maybe showing some resemblance to the pattern of bioavailability of P), as long as the soil parent material provides Si in sufficient quantity. Calcareous soils, which naturally have pH values> 8 in the soil solution (Pérez-Labrada et al., 2016), do not provide enough Si (Liang et al., 1994). Thus the fertilizer contributions with Si in crops in calcareous soils are beneficial (Zhang et al., 2017).

Another factor to consider regarding the availability of Si(OH)4 in the soil pore water is soil organic matter (SOM) and its dissolved forms. SOM have a profound impact on the availability of mineral elements (Diacono and Montemurro, 2010), either directly by chemical processes or indirectly by the promotion of bacteria and fungi that solubilize Si and other elements of soil minerals (Landeweert et al., 2001). An expected effect of SOM would be the adsorption of Al3+ through organic acids (Rustad and Cronan, 1995), which would decrease the Al-Si association and increase the concentration of Si(OH)4 available in the soil solution. The organic acids derived from SOM are also agents that promote dissolution in mineral surfaces (Drever and Stillings, 1997) so that in soils with silicon-rich parent materials or agricultural soils with the application of Si fertilizers would be very helpful. In nutrient solutions for soilless crops, the use of organic acids can also be useful to improve the solubility of fertilizers with silicon. On the other hand, indirect evidence is available that indicates that SOM is a factor that increases the bioavailability of Si for crop plants (Ding et al., 2008; Sun et al., 2017), thereby SOM management should be considered to increase the availability of Si for crops.

### THE MODEL

The information in the previous section highlighted the factors that can be subjected to management in a crop, both in soil and soilless, with the purpose of increasing the availability of Si(OH)4 for plants. In the crops grown in soil, a primary factor is the silicon content of soil's parent material. In the fertile soils of temperate and subtropical zones, Si inputs are rarely required in the fertilizers since the soil will undoubtedly provide the necessary amount. On the other hand, in the calcareous soils of arid and semi-arid regions, and in the soils of tropical regions subject to regimes of intense precipitation, the application of Si with fertilizers will be necessary, but also the consideration of the pH and organic matter management of the soil to ensure adequate availability of Si(OH)4.

In soilless crops, the critical factor to consider will be the concentration of Si in the irrigation water. Values below 28 mg L −1 Si point out the need to provide Si up to a maximum of 50 mg L −1 . The management of pH is the next factor to be considered. However, the data presented indicate that pH management aimed at ensuring the bioavailability of P in nutritive solution (that is, maintaining it between 5.5 and 7) will be adequate.

A Matlab-Simulink model (the archives are included in Supplementary Material) is presented below which allows verifying the impact of different environmental scenarios, both in a soil crop and in a soilless crop, using as a model tomato plants. There is also an example of the use of the software to obtain the estimated impact of the environmental variables on the absorption of Si by the tomato plants. The data presented in the previous parts of the manuscript can be tested in this model by verifying the result regarding the concentration of Si in the plants. The purposes of the use of the model are educational or technical, and from our perspective, the model can be useful in the agronomic management of Si in a tomato crop and, possibly applicable to other horticultural crops.

### Description of the Model

Tomato (Solanum lycopersicum L.) was used as a biological model to describe the distribution of silicon accumulation in the different organs. To describe the effects of the various environmental factors mentioned, the deterministic mathematical model proposed by Tap (2000) and modified by Juárez-Maldonado et al. (2014b) will be used as a basis.

The model consists of six state equations, using as inputs the radiation (PAR, µmol m−2 s −1 ), temperature (◦C), and CO2 concentration (µL L−1 ). The model allows to directly considering the effect of these three variables on the accumulation of silicon in the different organs of the tomato plant.

The scope of the model refers to environmental conditions where intense stress does not prevail since it is assumed that the growth rate of the plants will be a direct function of the irradiance and temperature.

Within the plant, silicon accumulates in different organs depending on the corresponding transpiration rates. Thus, it is necessary to calculate the transpiration by a tomato plant dynamically. For this, the equation 1 is used, based on the fact that a linear correlation can be considered between the biomass accumulated by the tomato plant and its transpiration (Juárez-Maldonado et al., 2014a).

$$Transipation = Biomass \ast plm/t\text{cg} \tag{1}$$

Biomass is the mass of the tomato plant in g m−2 ; plm is a parameter of the linear model (8.5714); and tcg is the time of crop growth (10279801 s).

Assuming that the tomato plant does not present a substantial accumulation of Si (Liang et al., 2015), the maximum absorption limit was set for the model at 1% (as SiO2) of the dry biomass (Miyake and Takahashi, 1978). In this condition, and while there is an unlimited availability of silicon in the soil solution, the accumulation of maximum total silicon (MSiT) in the plant would be as follows:

$$\text{MSiT} = \frac{Biomass}{100} \ast PD \tag{2}$$

Where PD is the planting density expressed in plants m−2 , which for this model was established in 3 plants m−2 . This plant density was used by Juárez-Maldonado et al. (2014b) and provide the best financial margin, high yield, and fruit quality (Peet and Welles, 2005).

The distribution of accumulated silicon in the tomato plant will then follow the different transpiration rates of its organs, that is, leaves> stem> fruits ≥ root. In the particular case of tomato, organ transpiration can be approximated to the following percentages of total transpiration: leaves = 90%, stem = 5%, fruits = 2.5%, and root = 2.5%.

Even though the potential availability of Si(OH)4 in the soil solution is 192.18 mg L−1 (Epstein, 1999; Liang et al., 2015), disponibility is affected by temperature, pH, and organic matter content of the soil. In addition to the factors that are modified with agricultural management such as soil moisture and soil profile.

According to the literature, the availability of silicon in soils is directly affected by soil temperature (Epstein, 1999; Liang et al., 2015). Although there is no clear explanation of how this behavior occurs, it is possible to approach it with a third-order model (Equation 3). This is due to the disponibility of silicon is between 8 and 35◦C, being its highest availability at 25◦ C.

$$T\mathfrak{B}\*Temp^3 + T\mathfrak{Z}\*Temp^2 + T1\*Temp + T0\tag{3}$$

Where T3, T2, T1, and T0 are the parameters of the thirdorder model (equivalent to −0.0003; 0.0127; −0.1093; and 0.1674 respectively), and Temp is the 0–30 cm soil temperature (◦C).

Concerning the SOM, it is known that there is a positive correlation with the availability of silicon (Ding et al., 2008; Sun et al., 2017). View from an agricultural perspective, soil is rich in organic matter when it has a concentration of 5%. An adjustment with a Michaelis-Menten function was used to describe the higher availability of Si, due to the effect of SOM. For this, the following Equation (4) was used.

$$Vma\overline{\times} \times \frac{OM}{(Km+OM)}\tag{4}$$

Where Vmax is the parameter of maximum availability of silicon due to organic matter normalized to 1 (Vmax = 1). OM is the amount of organic matter contained in the soil (%, w/w). And Km is the Michaelis-Menten parameter (Km = 2.5).

The pH is also a determining factor in the availability of silicon (Liang et al., 2015). This factor, as well as temperature, is related to the availability of silicon that can be approached to a thirdorder model. The availability of silicon in soil occurs in the pH range from 2 to 9, with a pH of 7 being the highest availability. Therefore, its effect can be described as follows:

$$pH\text{3} \* pH^{3} + pH\text{2} \* pH^{2} + pH\text{1} \* pH + pH\text{0}\tag{5}$$

Where pH3, pH2, pH1, and pH0 are the parameters of the thirdorder model (with values −0.0235, 0.325,−1.1563, and 1.2262, respectively). And pH represents the pH of the soil studied.

Therefore, the Si(OH)4 available (SiAv) to be absorbed by the tomato plant is described by the following equation:

$$\text{SiAv} = (\text{SiP} - \text{SiWater}) \ast ETem \ast EOM \ast EpH + \text{SiWater} \quad \text{(6)}$$

Where SiP is the maximum amount of silicon in a soil without polymerization [192.18 mg L−1 Si(OH)4]; ETem represents the effect of temperature on the availability of silicon (Equation 3); EOM describes the impact of organic matter on the availability of silicon (Equation 4); EpH represents the effect of pH on the availability of silicon (Equation 5); and SiWater is the amount of Si(OH)4 available in the irrigation water. The model supposes that under no condition will be the available Si(OH)4 be higher than the SiP value.

The accumulated Si (as SiO2) in the tomato plant (SiT) was determined with the silicon [Si(OH)4] available and the transpiration (Equation 1) in the following way:

$$\text{SiT} = \text{SiAv} \ast \text{Transjunction} \ast \text{SiSi} \tag{7}$$

Where SiSi is the fraction of cumulative Si in the plant of the total available Si(OH)4 (g).

In soilless cultivation conditions, the only source of Si(OH)4 for the crop will be the content of the irrigation water since there is no such source of soil mineral replacement as in the soil. Therefore, the accumulation of Si in the tomato plant in soilless culture (ASiTSC) will depend entirely on the transpiration of the plant (Equation 1) and the availability of Si(OH)4 in the irrigation water (SiWater). This relationship is expressed as:

$$\text{ASiTSC} = \text{Transjunction} \ast \text{SiWater} \ast \text{SiSi} \tag{8}$$

As previously described, the availability of silicon in soil depends on three primary conditions: pH, organic matter, and soil temperature. Of these conditions, it is feasible to modify the amount of organic matter or the pH. In the case of temperature, the easiest way would be to use covers as plastic mulches, which could increase the soil temperature by 3–4◦C (Ruíz-Machuca et al., 2015). Therefore, these factors can be considered as crucial factors to the agronomic management of silicon availability (Liang et al., 2015).

### Silicon Accumulation in Tomato

According to the simulations carried out using the proposed model, a soil with pH 7.0 and organic matter content of 6% can obtain the maximum availability of Si(OH)4, which can be > 4,500 mg m−2 at 15◦C; or > 8,300 mg m−2 at 25◦C (**Figures 1A,D**). Considering the availability, and two average temperature conditions of soil (15 and 25◦C), the highest availability of silicon is obtained with an average soil temperature of 25◦C (**Figure 1D**). On the contrary, when the organic matter content is low (<1%) along with a pH >8 (like a soil of a semi-arid region), the availability of Si(OH)4 in the soil can drastically decrease to <20 mg m−2 at 15◦C (**Figure 1A**), or <34 mg m−2 at 25◦ C (**Figure 1D**). These results describe the effect of pH, organic matter, and temperature factors on the availability of Si(OH)4 in soil. In addition to demonstrating the potential sensitivity of the availability of SiOH4 in the soil to the modifications on any of the conditions mentioned.

The availability of silicon in the soil will directly impact the accumulation of plants grown on it (Epstein, 2001). The higher availability of silicon derived from the factors evaluated (**Figure 1D**), results in a more significant accumulation of Si (as SiO2) based on the dry weight of the tomato plant (up to 1 %), as can be seen in **Figure 1E**. The **Figure 1E** fits well the reported Si concentration in tomato plants grown under pH 5.5 (Miyake and Takahashi, 1978) to pH 8.48 (Gunes et al., 2007).

The same behavior is observed concerning accumulated silicon per plant, reaching up to 13.7 g per plant, that represent the maximum accumulation of silicon for tomato plants under this conditions (**Figure 1F**). However, low availability of silicon in the soil can lead to a small accumulation of silicon in the plant. In the example the conditions of a soil corresponding to the situation of a semi-arid region, there would be an accumulation <1.3 g per plant for both soil temperature conditions (**Figures 1C,F**). This equals to a silicon concentration based on plant dry weight <0.1% (**Figures 1B,E**).

Since silicon accumulates in the different organs of the tomato plant as a function of the rate of transpiration, then the highest accumulation will be observed in the leaves, since they represent 90% of total transpiration. In fruits, a lower accumulation of silicon will be seen, since the rate of transpiration is little compared to that of the leaves (∼5%) (Leonardi et al., 2000). However, the availability of silicon in the soil will ultimately define the accumulation of silicon both in the entire plant and in its various organs.

Contrary to a crop established in soil, in soilless cultivation, e.g., hydroponics, the primary factor that will modify the availability of Si for the plant will be its concentration in the irrigation water used. It has been reported that irrigation water can have a Si content [as Si(OH)4] of 5–20 mg L−1 , while it is considered that an adequate concentration of Si would be 28 mg L −1 (Epstein, 1994, 1999; Liang et al., 2015). However, the transpiration rate of the crop will finally define the amount of SiOH4 absorbed and accumulated in the different organs. Since the growth of the plant and the proper distribution of biomass in the various organs will affect the rate of transpiration, then crop growth should be considered as an additional factor that will change the accumulation of silicon in a soilless crop system. Therefore, the environmental factors (PAR, CO2, air temperature) that affect the growth of the crop will, in turn, affect the accumulation of silicon in the different organs.

According to the simulations carried out, a low concentration of CO2 and a low incidence of PAR generate little accumulation of biomass in the tomato plant (**Figures 2A,D**). This same result is observed when the temperature of the environment changes, 30◦C generates biomass of up to 3,000 g per plant (**Figure 2A**); while 20◦C produces up to 2,240 g per plant (**Figure 2D**), at the highest conditions of PAR and CO2 concentration. As a consequence, the total transpiration of the plant is modified when

environmental factors changes, and therefore the accumulation of silicon. With an average temperature of 30◦C, the biomass accumulated in the fruits can represent around 60% of the total of the tomato plant, and the leaves less than 10%. However, when the temperature drops to 20◦C, the biomass distribution in the tomato plant changes. Under this condition, the biomass accumulated in the fruits is ∼33%, whereas in the leaves it increases up to 33%.

Since the leaves constitute the largest area of transpiration, modifying the temperature of the environment in which the tomato grows substantially alters the final accumulation of silicon. Therefore, according to the model, a temperature of 30◦C will lead to a lower accumulation of silicon in the tomato plant (**Figure 2B**), while at 20◦C there will be more significant accumulation (**Figure 2E**). The **Figure 2C** fits the reported Si concentration in tomato plants grown under temperatures of 28–32◦C (Cao et al., 2015). The result will be a higher level of silicon in dry weight of the tomato plant at low temperatures (**Figure 2F**), while high temperatures will decrease the concentration considerably (**Figure 2C**).

When considering the availability of silicon in the irrigation water, it can be observed that a condition of low Si content (5 mg L−1 as Si[OH]4) will result in less concentration and accumulation of silicon in the plant (**Figures 3A,B**). On the contrary, adequate availability of Si (28 mg L−1 as Si[OH]4) in the irrigation water will result in increased accumulation and concentration of silicon in the plants (**Figures 3C,D**). High accumulation appears to occur regardless of the environmental conditions in which the crop develops, when there is adequate availability of Si in the irrigation water. However, when conditions are favorable for the development of leaves in tomato plants, the maximum concentration of silicon for this species can be reached (**Figure 3D**).

The model presented focuses on the impact of external factors on the growth of tomato plants, under the assumption that as long as exists the availability of Si the plants absorb it and transport it at a rate proportional to the growth rate and transpiration rate. The point highlighted with the model is that the biofortification of the fruits with Si depends on the availability of the element in the form of Si(OH)4 both in soil and in soilless culture.

Sakurai et al. (2017) developed a successful dynamic model of Si absorption and transport for rice. The model is based on endogenous variables, as the activity of different transporters and the distribution of Si through different nodes in the entire plant. The model presented in this manuscript is focused on exogenous variables, susceptible to agronomic management both in field cultivation as in the greenhouse, and it has been used successfully

(C,D). For the simulation, 12 h of PAR and an average air temperature of 20◦C were considered.

to simulate the absorption of other elements in tomato and other crops (Juárez-Maldonado et al., 2014b; Ramírez-Pérez et al., 2018).

However, it must be taken into account that the presented model is used to describe the accumulation of silicon in the plant under relatively favorable environmental situations. The presence of stresses such as water deficit, salinity, deficiency of mineral nutrients and pathogens, results in loss of precision. With a certain amount of PAR and with a particular temperature regime, the stressed plants would have real biomass lower than the estimated by the model, which means an overestimation of the absorbed silicon.

As far as we know, except those published by (Sakurai et al., 2015, 2017) for monocotyledons, there are no similar dynamic models published about the absorption of silicon in dicots. The model described for the tomato crop is a first preliminar advance that we believe substantially can improve the understanding of some factors that regulate the bioavailability of silicon.

### CONCLUSIONS

After the results obtained from the presented model, the following is proposed:

1) When crops are grown in soil, the bioavailability of silicon can be increased by adding organic matter from organic amendments or humic substances, or by modifying the pH of the soil solution, using organic or inorganic acids, to be the closest to 7.0.

### REFERENCES

## AUTHOR CONTRIBUTIONS

All authors were responsible for processing information and manuscript writing. AB-M, AJ-M, FP-L, and AM-D: Conceptualization; AJ-M, ML-P, LR-P, and AM-D: Model design and implementation; ML-P, FP-L, SG-M, LG-D, JG-M, and AB-M: Manuscript drafting. All authors read and approved the final manuscript.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018. 00647/full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 López-Pérez, Pérez-Labrada, Ramírez-Pérez, Juárez-Maldonado, Morales-Díaz, González-Morales, García-Dávila, García-Mata and Benavides-Mendoza. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Genetic Basis of Variation in Rice Seed Storage Protein (Albumin, Globulin, Prolamin, and Glutelin) Content Revealed by Genome-Wide Association Analysis

Pingli Chen1 , Zhikang Shen1 , Luchang Ming1 , Yibo Li1 , Wenhan Dan1 , Guangming Lou1 , Bo Peng1 , Bian Wu1 , Yanhua Li2 , Da Zhao1 , Guanjun Gao1 , Qinglu Zhang1 , Jinghua Xiao1 , Xianghua Li1 , Gongwei Wang1 and Yuqing He1 \*

1 National Key Laboratory of Crop Genetic Improvement and National Center of Plant Gene Research, Huazhong Agricultural University, Wuhan, China, 2 Life Science and Technology Center, China National Seed Group Co., Ltd., Wuhan, China

### Edited by:

Huixia Shou, Zhejiang University, China

### Reviewed by:

Umesh K. Reddy, West Virginia State University, United States Shahidul Islam, Murdoch University, Australia

> \*Correspondence: Yuqing He<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Plant Nutrition, a section of the journal Frontiers in Plant Science

Received: 21 January 2018 Accepted: 18 April 2018 Published: 09 May 2018

#### Citation:

Chen P, Shen Z, Ming L, Li Y, Dan W, Lou G, Peng B, Wu B, Li Y, Zhao D, Gao G, Zhang Q, Xiao J, Li X, Wang G and He Y (2018) Genetic Basis of Variation in Rice Seed Storage Protein (Albumin, Globulin, Prolamin, and Glutelin) Content Revealed by Genome-Wide Association Analysis. Front. Plant Sci. 9:612. doi: 10.3389/fpls.2018.00612 Rice seed storage protein (SSP) is an important source of nutrition and energy. Understanding the genetic basis of SSP content and mining favorable alleles that control it will be helpful for breeding new improved cultivars. An association analysis for SSP content was performed to identify underlying genes using 527 diverse Oryza sativa accessions grown in two environments. We identified more than 107 associations for five different traits, including the contents of albumin (Alb), globulin (Glo), prolamin (Pro), glutelin (Glu), and total SSP (Total). A total of 28 associations were located at previously reported QTLs or intervals. A lead SNP sf0709447538, associated for Glu content in the indica subpopulation in 2015, was further validated in near isogenic lines NIL(Zhenshan97) and NIL(Delong208), and the Glu phenotype had significantly difference between two NILs. The association region could be target for map-based cloning of the candidate genes. There were 13 associations in regions close to grainquality-related genes; five lead single nucleotide polymorphisms (SNPs) were located less than 20 kb upstream from grain-quality-related genes (PG5a, Wx, AGPS2a, RP6, and, RM1). Several starch-metabolism-related genes (AGPS2a, OsACS6, PUL, GBSSII, and ISA2) were also associated with SSP content. We identified favorable alleles of functional candidate genes, such as RP6, RM1, Wx, and other four candidate genes by haplotype analysis and expression pattern. Genotypes of RP6 and RM1 with higher Pro were not identified in japonica and exhibited much higher expression levels in indica group. The lead SNP sf0601764762, repeatedly detected for Alb content in 2 years in the whole association population, was located in the Wx locus that controls the synthesis of amylose. And Alb content was significantly and negatively correlated with amylose content and the level of 2.3 kb Wx pre-mRNA examined in this study. The associations or candidate genes identified would provide new insights into the genetic basis of SSP content that will help in developing rice cultivars with improved grain nutritional quality through marker-assisted breeding.

Keywords: GWAS, storage protein, grain quality, endosperm, nutrition, Oryza sativa L.

### INTRODUCTION

fpls-09-00612 May 7, 2018 Time: 23:5 # 2

Rice is one of the major staple cereal foods and is an important source of total protein in human food. SSP account for approximately 8% of the dry grain weight and are the second most abundant ingredient after starch in rice. Rice has the lowest protein content among cereal grains, but net protein utilization is highest (Juliano, 1992). High-protein rice is likely to increase human nutrition in poor families, especially where rice is a staple food. Therefore, increasing the SSP content has become one of the main breeding objectives in improving nutritional quality in rice.

The SSP in rice can be classified into four fractions: albumin, globulin, prolamin, and glutelin, according to differences in solubility. Glutelin encoded by 15 genes accounts for as much as 80% of the total SSPs and is concentrated in the milled fraction, whereas prolamin, the most evenly distributed protein, accounts for less than 5% (Yamagata et al., 1982). Based on amino acid sequence similarity, glutelins are classified into four subfamilies: GluA, GluB, GluC, and GluD (Kawakatsu et al., 2008). In rice, SSP genes have been cloned and characterized mostly by mutant screening (Ren et al., 2014). Glutelins are synthesized in the rough endoplasmic reticulum as a 57 kDa precursor. Previous studies have identified rice 57H mutants that accumulate relatively high levels of 57 kDa pro-glutelin and have floury/opaque endosperm phenotype (Wang et al., 2009; Ren et al., 2014). Of 57H mutants, only gpa3, Osvpe1, and OsRab5a have been successfully cloned (Wang et al., 2009; Wang Y. et al., 2010; Ren et al., 2014). Prolamins are encoded by a multigene family of 34 gene copies. Based on the molecular mass, prolamins classified into three groups: 10 kDa prolamin (RP10), 13 kDa prolamin (RM1, RM2, RM4, and RM9), and 16 kDa prolamin (RP16) (Yamagata et al., 1982; Kawakatsu et al., 2008). Both albumin and globulin are concentrated in the bran and polishing during milling removes a major portion of these proteins (Shewry, 2007). Globulin is also easily digested (Yamagata et al., 1982; Zhang et al., 2008), and only a limited number of genes has been cloned and characterized (Bhullar and Gruissem, 2013). RA16 and RA17 that are associated with seed allergenic protein have been reported as albumin genes in previous studies (Adachi et al., 1993). The nutritional value of glutelin is higher than prolamin because it has a greater digestive capacity by humans and higher lysine content in rice (Ogawa et al., 1987). However, patients with kidney disease and diabetes need low glutelin diet (Mochizuki and Hara, 2000; Nishimura et al., 2005; Morita et al., 2009). Some proteins that belong to the albumin and globulin are considered to be allergenic. Collectively, emphasis in rice breeding should not only be on the concentration, but also on the quality of rice protein.

Seed storage protein is quantitatively inherited, and is affected by growing environment (Shewry, 2007). QTL mapping based on molecular markers and linkage maps has always been a common method of genetic studies (Li et al., 2014; Peng et al., 2014). Many QTL for crude protein content in rice have been reported (Aluko et al., 2004; Wang et al., 2008; Lou et al., 2009; Yu et al., 2009; Peng et al., 2014), but fewer studies have investigated the individual protein fractions in milled rice. qPC1 (OsAAP6), controlling the natural variation in SSP content, has been cloned using a mapbased cloning strategy (Peng et al., 2014). Zhang et al. (2008) identified 16 QTL for contents of crude protein and the four protein fractions.

Genome-wide association study (GWAS) by means of single nucleotide polymorphism (SNP) has become the method of choice for investigation of the genetics of important traits in Arabidopsis thaliana (Chan et al., 2011), rice (Huang et al., 2010), maize (Xiao et al., 2016), sorghum (Morris et al., 2013), and others (Ogura and Busch, 2015). Although GWAS is widely used in genetic analysis of grain quality traits in rice, such as gelatinization temperature, amylose content, grain appearance, and milling quality (Borba et al., 2010; Huang et al., 2010), few studies have used this approach to investigate total SSP in rice (Huang et al., 2012; Bryant et al., 2013). The analysis of genetic basis of nutritional quality has been reported in maize (Deng et al., 2017). Maize opaque2 (o2) mutation could increase free lysine levels and create the foundation for quality protein maize (QPM) breeding. Combining GWAS and linkage mapping, a gene duplication at the 27-kDa γ-zein locus qγ27 was identified (Liu et al., 2016). qγ27 increases the level of 27-kDa γ-zein gene expression in QPM, which is essential for endosperm modification. GWAS on amino acids has been carried, and 247 and 281 significant loci were identified in two different environments (Deng et al., 2017). However, no studies have been reported about the genetic basis of the four SSP fractions by GWAS in rice.

In this study, we performed GWAS of Alb, Glo, Pro, Glu, and total SSP in milled rice using 527 Oryza sativa accessions grown in two environments with an aim to identify loci involved in the genetics. Haplotype analysis and expression pattern of candidate genes then provided valuable in better understanding the genetic basis of variation in SSP content. The results shown here promote our understanding of the genetic basis of the storage protein groups, should be of use for breeders attempting to improve nutritional quality by means of marker assisted selection.

### MATERIALS AND METHODS

### Plant Materials, Field Experiments, and Trait Measurements

A diverse worldwide collection of 527 O. sativa landraces and elite accessions (Supplementary Table S1) was used in this study. Structural analysis indicated that the entire collection belonged to nine subpopulations: indI, indII, indica intermediate, Tej, Trj, japonica intermediate, Aus, VI, and intermediate (Chen et al., 2014) and are available at the RiceVarMap1 . The indica subpopulation (indI, indII, and indica intermediate) included 294 accessions, whereas the japonica subpopulation (Tej, Trj, and japonica intermediate) included 155 accessions.

Lines were planted in two environments in Hubei province: Ezhou in 2014 (Env. 1) and Wuhan in 2015 (Env. 2). The sowing dates were 25 May in both years. Seedlings about 25 days old were transplanted to the field. There were three

1http://ricevarmap.ncpgr.cn

rows with 10 plants each in each plot. The planting density was 16.5 cm between plants within a row, and 26.4 cm between rows. Field management basically followed recommended practice of agriculture, with fertilizer applied (per hectare) as follows: 48.75 kg nitrogen, 58.5 kg phosphorous, and 93.75 kg potassium as the basal fertilizer; 86.25 kg nitrogen at the tilling stage; and 27.6 kg nitrogen at the booting stage.

At maturity, three plants in the middle of the second row of each accession were harvested and bulked. Dry seeds were threshed in bulk and the rough rice was air-dried, and stored at room temperature for 3 months, and then stored at 4◦C. 50 g of rough rice were dehulled into brown rice using a TR 200 dehuller (Kett, Tokyo, Japan). The embryo and aleurone layer of brown rice were removed into milled rice through a Pearlest mill (Kett, Tokyo, Japan). The rice was ground into flour with a CT 410 Cyclotec mill (FOSS, Hillerod, Denmark), passed through an 80-mesh sieve and stored at −20◦C until the Alb, Glo, Pro, and Glu contents were determined based on previously published previous methods (Kumamaru et al., 1988). Briefly, 0.1 g sample of milled rice flour was placed in a centrifugation tube with 1.0 ml solvent containing 10 mM Tris–HCl buffer (pH7.5) for Alb extraction; 1.0 ml solvent containing 1 M NaCl, for Glo extraction; 60% n-propanol containing 1 mM EDTA-2Na, for Pro extraction; and 0.05 M NaOH for Glu extraction. The mixture was stirred for 2 h at room temperature, and extracts were separated from residues by centrifugation at 12,000 rpm for 15 min at 4◦C The procedure was repeated three times. The extracts were stored at −20◦C until further analysis. The contents of each fraction were determined by the Coomassie brilliant blue G-250 dye-binding method (Bradford, 1976) using bovine serum albumin as a standard, and quantitative analysis was carried out using Infinite M200 (Tecan Group, Männedorf, Switzerland) (Peng et al., 2014). Total SSP was the sum of the Alb, Glo, Pro, and Glu contents. The 2-year field experiment was designed with three replicates per year. The average SSP contents across three replicates within 1 year were used for GWAS. The SSP contents of the 527 O. sativa accessions are listed in Supplementary Table S1. Amylose content was measured as previously described method (Tan et al., 1999).

### Genome-Wide Association Study

All 527 accessions were genotyped via sequencing (Chen et al., 2014). SNP information was available on RiceVarMap (see foot note text 1), a comprehensive database for rice genomics. The physical locations of the SNPs were identified based on the rice annotation version 6.1 of variety Nipponbare from Michigan State University. A total of 3,916,415 SNPs in the whole population; 2,767,159 SNPs in the indica subpopulation; and 1,857,845 SNPs in the japonica subpopulation (minor allele frequency ≥0.05; number of accessions with minor alleles ≥6) was used for GWAS (Chen et al., 2014). A linear mixed model (LMM) was used for detecting associations using Fast-LMM (Lippert et al., 2011). Population structure was controlled using a kinship matrix constructed with all SNPs (Chen et al., 2014). Effective independent SNPs were detected (Li et al., 2012), and were 757,578 in the whole population; 571,843 in the indica subpopulation; and 245,348 in the japonica subpopulation. The thresholds were set at a P-value of 5.0 × 10−6 to identify significant association signals. To obtain independent association signals, multiple SNPs, exceeding the threshold in a 5 Mb region, were clustered based on an r 2 of LD ≥ 0.25, and SNPs with the minimum P-value in a cluster were deemed to be lead SNPs.

### Statistical Analysis

Based on the standardized disequilibrium coefficients (D'), linkage disequilibrium (LD) was investigated. LD heatmaps were constructed using the TASSEL5.02 program and R package "LDheatmap"3 . Statistical analysis, including a correlation analysis, was conducted using IBM SPSS Statistics 22.0. Differences in SSP values were examined by Student's t-tests. Broad-sense heritability (H2 ) for each phenotype was estimated using repeatability between 2 years of phenotypic data, calculated as the variance among variety grand means divided by their total phenotypic variance.

### Candidate Genes and Haplotype Analysis

Candidate genes within a 200 kb genomic region ( ± 100 kb from the lead SNP) in the associated loci were selected based on (i) biochemically related proteins or protein clusters; (ii) homologous genes with known function, and (iii) expression profiles. The genotypes of RP6, RM1, Wx, PROLM1, and other three candidate genes in the 527 rice accessions were obtained from the RiceVarMap database (see foot note text 1). The haplotypes were classified according to all SNPs (except sites in intron) including their intragenic region and 2 kb upstream with an MAF > 0.05 in a candidate gene. There were at least five rice accessions in the haplotypes for comparative analysis. One-way ANOVA and Student'st-tests were applied to compare differences in SSP content among all possible haplotype pairs.

### RNA Extraction and Quantitative RT-PCR Analysis

According to the manufacturer's instructions, the total RNA was extracted from rice different tissues using TRIzol reagent (Invitrogen). About 3 µg of RNA sample was processed by RNase-free DNaseI (Invitrogen) and reverse transcribed using M-MLV reverse transcriptase (Invitrogen) with Oligo(dT)15. Quantitative RT-PCR was carried out using Fast Start Universal SYBR Green Master (Rox) superMIX (Roche, Mannheim, Germany) in a ViiA 7 Real-Time PCR system (Applied Biosystems), according to the manufacturer's introductions. Measurements were obtained using the relative quantification method. Actin was used as a reference gene in the qRT-PCR experiments. The experiment was designed with three biological replicates and three technical replicates per material. Error bars indicate standard error. The measurements were obtained using the relative quantification method. The significant difference was analyzed statistically by One-way ANOVA and Student's t-tests. All primers for qRT-PCR analysis are listed in Supplementary Table S2.

2https://tassel.bitbucket.io/

3https://www.r-project.org/

### RESULTS

fpls-09-00612 May 7, 2018 Time: 23:5 # 4

### Phenotypic Variation and Heritability of SSP Content

The result analysis revealed a large variation in all phenotypes evaluated and the traits appeared to be normally distributed (**Figures 1A–E** and Supplementary Table S3). Average SSP contents in Env. 1 and 2 were 70.9 and 53.7 mg/g; Alb was 2.9 and 3.3 mg/g; Glo was 5.8 and 5.3 mg/g; Pro was 2.8 and 2.1 mg/g; and Glu was 59.4 and 43.5 mg/g, respectively (**Figures 1A–E**). Glu accounted for approximately 80% of total SSP; Alb and Pro each accounted for about 5%; and Glo accounted for about 10% (**Figure 1F**). Compared with other three storage protein contents, the average content of Pro was the lowest, but the variation of Pro was largest in the whole population and each subpopulation in both environments. In four storage protein, Total SSP showed the lowest heritability (29.5%), whereas Pro showed the highest (76.8%) (Supplementary Table S3).

Correlation coefficients between each pair of components, and between components and total SSP were significant and positive in both environments, except for those between Glo and Pro in any environment and between Alb and Glo and between between Alb and Pro in Env. 2 (Supplementary Table S4). High correlations were found only between Glu and total SSP in both environments with coefficients of 0.99 in Env. 1 and 0.98 in Env. 2.

### Genome-Wide Association Study for SSP Contents

We performed GWAS on the entire population and on the indica and japonica subpopulations for each year. The FaST-LMM program reduced the effect of population structure (Yang et al., 2014). Quantile-quantile plots of all five traits for the whole population, and indica and japonica subpopulations are illustrated in **Figure 2** and Supplementary Figures S1, S2. Some associations were detected in different subpopulations, and some of the associations for different traits were in the same chromosomal regions. Any two lead SNPs within a 100 kb region were considered to be a single association locus.

The association analysis for the whole population identified 34 loci (phenotypic variance >10%) associated with three traits with a suggestive threshold value at 5.0E-06 (**Table 1**). Most of them were detected for Alb (21 associations) and Pro (12 associations). Lead SNPs for Alb were widely distributed in the rice genome: chromosome 1, 2, 3, 4, 5, 6, 7, 8, and 9, with chromosomes 6 and 7 having more associations. Associations explained phenotype variation of 10.4–20.9%, with the association on chromosome 6 (sf0601764762) making the largest effect. For Pro, associations accounting for 10.1–34.6% of the phenotypic variance were identified on chromosomes 2, 5, 7, 10, and 11, with chromosome 5 and 7 exhibiting more associations. Only one lead SNP sf0317000156 on chromosome 3 with a phenotype variation of 14.4% was detected for Glo. No lead SNP with phenotype variation of more than 10% was detected for Glu or total SSP.

A large number of peaks (phenotypic variance >10%) were also detected by GWAS in the indica (33 associations) and japonica (40 associations) subpopulations (Supplementary Table S5). For Alb, 23 associations were distributed on all 12 rice chromosomes, except chromosome 3 and 10, but only three associations were identified in the indica subpopulation. Alb associations identified in the indica and japonica subpopulations explained phenotype variation of 10.5– 12.3, and 11.3–41.6%, respectively. Four SNPs, sf0142207782, sf0209990680, sf0605251091, and sf0804866973 individually explaining more than 30% of the Alb variation, were detected in the japonica subpopulation in Env. 1. For Glo, 14 associations were identified on chromosomes 1, 3, 4, 5, 6, 9, 11, and 12, and equal numbers of associations were detected in the indica and japonica subpopulations, explaining 10.3–16.8 and 12.9–17.0% of the variation, respectively. For Pro, there were 23 associations, involving chromosomes 1, 3, 4, 5, 7, 9, 11, and 12. Among them 14 (with phenotype variation of 10.4–38.3%) and 9 (with phenotype variation of 15.5–22.6%) associations were identified in the indica and japonica subpopulations, respectively. Ten associations for Glu were detected on chromosomes 1, 4, 5, 6, 7, 8, 9, and 10, in both subpopulations. For total SSP, only three associations on chromosomes 10 and 11 with phenotype variation of 10.4–14.8% were detected in the indica subpopulation.

Among the 107 associations detected in the whole population and in the indica and japonica subpopulations, 16 were detected in different populations and nine involved two or three different traits (**Table 1** and Supplementary Table S5). Examples include lead SNP sf0519612378 with phenotypic variance of 14.9– 36.6% that was detected in the whole population and indica subpopulation in Env. 1 and 2; and lead SNP sf1022972496 detected in indica subpopulation in Env. 1 was associated with phenotypic variances of 12.0 and 14.8% for traits Glu and total SSP, respectively. Seven associations in the whole population and four associations in the indica subpopulation were detected both in both environments. For Pro, two lead SNPs (sf0514987630 and sf0515211855) in the whole population and one lead SNP (sf0515706446) in the indica subpopulation were detected in both environments. The significance levels of the associations ranged from P = 5.0E-06 to P = 8.6E-16, P = 4.9E-06 to P = 9.0E-08, P = 5.0E-06 to P = 1.3E-17, P = 2.1E-06 to P = 1.4E-08, and P = 4.5E-06 to P = 3.1E-08 in LMM for Alb, Glo, Pro, Glu, and total SSP, respectively, and the most significant association for sf0706363663 located in chromosome 7 (**Table 1**).

### Co-localization of Associated Sites With QTLs Previously Reported and Grain Quality-Related Genes

There were many overlaps between the present associations detected by GWAS and reported QTLs or intervals related to SSP content in rice. A total of 28 associations from this study were located at previously reported QTLs or intervals (shown in **Table 1** and Supplementary Table S5 with corresponding references) of which 10 and eight associations were for Alb and Pro, respectively.

In the indica subpopulation, lead SNP sf0709447538 in a notable hotspot region at the interval of 9.1–9.5 Mb, explaining 12.7% of the Glu variation, was detected on chromosomes 7

in 2015 (**Figures 3A,B** and Supplementary Table S5). Interestingly, we found that the lead SNP sf0709447538 was overlapped with the amino acid content QTLs (7–4, 7–5, and 7–6), identified in a previous study using an F9 recombinant inbred line population, which derived from a cross between Zhenshan97 (ZS97) and Delong208 (DL208) (Wang et al., 2008; Zhong et al., 2011). To validate the QTL, we developed near-isogenic lines (NILs) (**Figure 3C**). NILs of QTL were developed by successive crossing and backcrossing ZS97 (high protein content) and DL208 (low protein content), three times (BC3) to ZS97. The QTL was selected by two molecular markers MRG186 and MRG4499 (Supplementary Table S2). Self-pollinating the BC3F1 plants heterozygous for this fragment produced NIL(ZS97) and NIL(DL208). Analysis of NIL(ZS97) and NIL(DL208) showed that NIL(ZS97) was significantly higher Glu than NIL(DL208), which was the same as the phenotype in two parents (**Figure 3D**). The result indicated that the QTL was reliable, which helps further identify the underlying genes and their genetic basis.

On the other hand, 13 associations were detected in regions close to previously identified grain-quality-related genes. Five genes (PG5a, Wx, RM1, RP6, and AGPS2a) were less than 20 kb from lead SNPs. PG5a, RM1, and RP6 were associated for Pro; Wx

(P-value). Lead SNPs in significant peaks are red. The horizontal dotted line indicated the genome-wide significance threshold (P = 5.0E-06). Total, total SSP content.

and AGPS2a was associated with Alb. Glutelin genes GluA1 and OsAAT2 were associated with total SSP in Env. 2. Additionally, several lead SNPs were located close to starch-metabolism-related genes, such as PUL, ISA2, Wx, GBSSII, OsACS6, and AGPS2a. SNPs located close to other reported grain quality-related genes (Prol14, RA17, RA16, and PGL) were also examined using a LMM (**Table 1** and Supplementary Table S5).

### Haplotype Analyses for the Reported Genes RP6 and RM1

The association SNP sf0705739605 for Pro was each <10 kb away from two prolamin genes RP6 and RM1 reported to encode prolamins in rice (Wen et al., 1993; Kawakatsu et al., 2009). In particular, the lead SNP with highly significant P-values (P = 6.5E-11 in Env. 1 and P = 5.4E-16 in Env. 2) was identified in the whole population and indica subpopulation (**Figures 4A–C**, **Table 1**, and Supplementary Table S5). The association explained 24.7 and 26.3% of the phenotypic variances in the whole population in Env. 1 and 2, respectively. Lead SNPs sf0705735351 and sf0705739605 and all polymorphic sites in RP6 and RM1 were in high linkage disequilibrium (in high LD with each other; r 2 = 0.94–0.99) with most polymorphic sites (**Figure 4H**). We performed haplotype analyses for RP6 and RM1 and identified three main haplotypes (Hap1-3) at each locus (Supplementary Table S6). Hap1 and hap2 of RP6 and hap2 and hap3 of RM1 were not identified in japonica. Hap2 of RP6 and hap3 of RM1 had significantly higher Pro than the alternative haplotypes in both

TABLE 1 | Associated single nucleotide polymorphisms (SNPs) identified by linear mixed model (LMM) method in the whole population.

Pop, population; Chr., chromosome; P, P-value estimated in LMM; P.V (%), proportion of phenotypic variance explained; Total, total seed storage protein (SSP) content. aAssociated SNPs within 100 kb for the same trait are considered the same locus; for Alb, Glo, and Pro, only lead SNPs with P.V (%) more than 10% are shown; bold scripts indicate detections in different populations and italic indicates detection in different traits.bNegative value means the gene is upstream of the SNP site.cRefers to interval or QTL reference.

environments (**Figures 4D,F**). In the region contained coding region and 2 kb upstream of RP6 and RM1, 17 and 26 SNPs were found, respectively (Supplementary Table S6). At the RP6 locus, there were two synonymous SNPs, two non-synonymous SNPs in the exon, three SNPs in the 30 untranslated regions, and 10 substitutions in the 2 kb cis-regulatory region. At the RM1 locus, there were five synonymous SNP, no non-synonymous SNPs in the exon, two SNPs in the 50 and 30 untranslated regions, and 19

plots of LMM for Pro in the all accessions in 2015 (A). (B) Local Manhattan plot surrounding the peak in 2015 on chromosome 7. Arrow indicates the position of the lead peak. The corresponding colors of r 2 represent linkage disequilibrium levels. (C) Plant architectures of near-isogenic lines. (D) Phenotypes of Glu of two parents ZS97, DL208, NIL(ZS97), and NIL(DL208). ∗∗Indicates the differences of Glu between two materials are significant at P < 0.01.

substitutions in the 2 kb cis-regulatory region. To get an overview of the expression profiles of RP6 and RM1, the CREP database4 , a website that contains the dynamic gene expression profile of indica rice, was searched (Wang L. et al., 2010). The results showed that RP6 and RM1 displayed high-level expressions in endosperm but low-level expressions in other tissues in ZS97 (**Figures 4E,G**). Considering the complexity of population structure and genetic background, we checked the expressions of RP6 and RM1 with different haplotypes in 20 and 34 accessions randomly chosen from indica group, respectively. Using qRT-PCR analysis, The results indicated that in the endosperm at 7 days after pollination (d.a.p.), expression levels of RP6 in hap2 accessions were much higher than those in hap3 accessions (**Figure 4I** and Supplementary Table S7), and expression levels of RM1 in hap3 accessions were much higher than those in hap1 accessions (P < 0.01) (**Figure 4J** and Supplementary Table S7). These results show that two genes might be good candidates for the GWAS locus. In conclusion, two genes had high-level and specific expressions in endosperm of rice, and genotypes with

### Analysis of One Candidate Gene Wx

We also found a highly significant association signal for Alb involving sf0601764762 (P = 2.0E-08 in Env. 1 and P = 8.6E-16 in Env. 2) on chromosome 6 (**Table 1** and **Figures 5A–C**). The lead SNP, explaining 12.9 and 20.9% of the phenotypic variances in the whole population in Env. 1 and 2, respectively, was located in the first intron of Wx (**Table 1** and **Figure 5D**). Other SNPs in Wx showed different LD associations with sf0601764762 (**Figure 5D**). Wx is the most important genetic determinant of amylose content (Tian et al., 2009). We identified eight major Wx haplotypes with 23 SNPs (**Figures 5E,F** and Supplementary Table S6). Hap1, 6, and 7 showed lower Alb values than hap4 and 5 in both environments. The lead SNP sf061764762 is located to the first intron of Wx, the major gene determining starch content. This SNP generates alleles Wxa with a normal GT sequence at the 5 0 splice junction of intron 1, and Wxb with a G to T mutation in intron 1. Wxa and Wxb produce a mature 2.3 kb Wx mRNA and a 3.3 kb Wx pre-mRNA, respectively (Wang et al., 1995).

higher Pro exhibited much higher expression levels in indica group in this study.

4http://crep.ncpgr.cn/crep-cgi/home.pl

signals of RP6 (E) and RM1 (G) in various tissues of ZS97 based on the microarray data. The y-axis represents the expression signals. (H) Representation of pairwise r 2 values among polymorphic sites in RP6 and RM1. The lines in red represent lead SNP. Expression levels of RP6 (I) and RM1 (J) in the endosperm at 7 days after pollination in indica group. Error bars, SE of 3 replicates. ∗∗Indicates the differences of expression levels between two haplotypes are significant at P < 0.01. Hap, haplotype; HD, heading date; DAP, day after pollination.

We checked the quantity of 2.3 and 3.3 kb Wx RNA with four haplotypes in 35 and 37 accessions randomly chosen from indica and japonica, respectively. Using qRT-PCR analysis, it showed that quantity of 2.3 kb Wx RNA in the endosperm at 7 d.a.p. in hap7 indica accessions were much higher than those in hap2 indica accessions, hap3 japonica accessions and hap5 japonica accessions (P < 0.01). In contrast, quantity of 3.3 kb Wx RNA in hap7 indica accessions were much lower than those in hap2 indica accessions and in hap5 japonica accessions (**Figures 5G,H** and Supplementary Table S7). With highest Alb content among four haplotypes, hap5 had relatively high 3.3 kb Wx RNA quantity and low 2.3 kb Wx RNA quantity. However, with lowest Alb content among four haplotypes, hap7 had relatively low 3.3 kb Wx RNA quantity and highest 2.3 kb Wx RNA quantity. It was found that Wx exhibited higher level quantity of 2.3 kb Wx mRNA in endosperm in ZS97 (hap7) than Minghui 63 (hap2), two indica cultivars (**Figure 5I**). We further compared the quantity of 2.3 and 3.3 kb Wx RNA in various tissues between two varieties Zhonghua 11 (hap2) and ZS97 by quantitative RT-PCR (**Figure 5J**). The results showed that ZS97 had lower 3.3 kb Wx RNA quantity and higher 2.3 kb Wx RNA quantity than Zhonghua 11 in endosperm of 7 and 14 d.a.p. However, Wx displayed very low-level expressions in stem, sheath, and flag leaf in both ZS97 and Zhonghua 11. Then we measured amylose content of the corresponding accessions in 2015 (Supplementary Table S7). The correlation analysis among Alb, amylose content, quantity of 2.3 and 3.3 kb Wx RNA was performed, and the results are presented in **Table 2**. Amylose content had significant correlations with both quantity of 2.3 and 3.3 kb Wx RNA, but amylose content was positively correlated with quantity of 2.3 kb Wx RNA and was negatively correlated with quantity of 3.3 kb Wx RNA, which was consistent with previous study (Wang et al., 1995). Significant and negative correlations were observed between Alb and amylose content or quantity of 2.3 kb Wx RNA (**Table 2**), suggesting that Wx may negatively regulate Alb. However, Alb and quantity of 3.3 kb Wx RNA had no

significant correlations. We speculated that Wx might influence Alb content.

### Analyses of Four Candidate Genes

Association analysis at the haplotype level greatly increases mapping power (Han et al., 2016). We re-detected these grain-quality-related genes mentioned above at the haplotype level (Supplementary Table S8). Four SSP loci (OsAAT2,

Alb of 46 accessions in 2015 was measured; AC, amylose content of 44 accessions in 2015; 2.3 kb, relative expression level of 2.3 kb Wx; 3.3 kb, relative expression level of 3.3 kb Wx. ∗P < 0.05 and ∗∗P < 0.01.

RA17, RM1, and RP6) and four starch-metabolism-related genes (AGPS2a, ISA2, PUL, and Wx) were detected in both environments; however, GBSSII, GluA1, and RA16 were only detected in one environment. The order of magnitudes of the P-values for these genes was greatly reduced compared with that at the SNP level. To further verify the association possibility, we validated some of candidate genes detected via GWAS by haplotype analyses and expression profiles in public databases.

To determine whether novel functional loci were implicated by GWAS, we investigated SNP sf0515211855 (P = 6.5E-07 in Env. 1 and P = 2.7E-11 in Env. 2), in Chr.5 (14.7- 15.8 Mb) (**Figures 6A–C**). This lead SNP was significantly associated with Pro and explained 8.5 and 13.7% of the phenotypic variance in the whole population in Env. 1 and 2, respectively (**Table 1**). SNPs in the region were grouped into two LD block, and sf0515211855 was in the first block. Within this block, there were a gene cluster of 18 prolamin precursors and a candidate gene encoding an expressed protein (LOC\_Os05g25500) (**Figure 6B**). Within the

represents the expression signals. HD, heading date; DAP, day after pollination.

cluster we focused on PROLM1 (LOC\_Os05g26240), located close to the lead SNP, and identified twelve major haplotypes (Hap1-12) (**Figures 6D,E** and Supplementary Table S6). Hap4, 5, and 6 had significantly higher Pro than other haplotypes in both environments. Microarray data (see foot note text 4) indicated that PROLM1 had extremely high expression enrichment in endosperm of ZS97 at the ripening stage (**Figure 6G**). We obtained five major haplotypes (Hap1-5) for LOC\_Os05g25500 (Supplementary Table S6). Hap3 and 5 had significantly higher Pro content than Hap1 and 4 in both environments (**Figure 6F**). LOC\_Os05g25500 exhibited higher expression levels in endosperm than other tissues of ZS97 (**Figure 6H**). We focused on other two candidate genes LOC\_Os03g29750 (encoding an expressed protein) and LOC\_Os02g13130 (encoding a KH domain-containing protein) for Glo in indica group in 2015 (**Figure 7A**) and Glu in all group in 2014 (**Figure 7B**), respectively. Hap5 of LOC\_Os03g29750 had significantly higher Glo than other haplotypes, and LOC\_Os03g29750 showed very high-level expression in endosperm of ZS97 (**Figure 7C**). Hap4 of LOC\_Os02g13130 had significantly higher phenotypes than other haplotypes. For LOC\_Os02g13130, homologous with maize gene encoding high molecular weight glutenin subunit x, it had high-level expressions in most tissues and organs (**Figure 7D**). These results indicated that these genes might be good candidates for the GWAS locus.

HD, heading date; DAP, day after pollination.

### DISCUSSION

### Phenotypic Variation and Trait Correlation

The contents of the four components of SSP and the total SSP were normally distributed in two environments and influenced by environment in this study (**Figure 1**). The indica subpopulation showed wider variation in Pro than the japonica subpopulation, but narrower variation compared to japonica for Alb. These results indicated that there were probably different genetic bases underlying Alb and Pro in the two subpopulations. The contents and proportions of different components as well as total SSP in the present study were in agreement with those reported by Kumamaru et al. (1988). Glu is highest composition of SSP in rice endosperm and has more essential amino acid required for human (Tanaka et al., 1995), suggested that it is more effectively to improve protein content and nutritional value of rice by increasing Glu than other three SSPs.

In previous studies, different conclusions were reached on the heritability of SSP content in rice. Hillerislambers et al. (1973) reported that the heritability of SSP content was 13.3–37.2% in rice. However, Shenoy et al. (1991) obtained values as high as 71% and considered that it is effective to select protein content phenotypes in early generation in the protein content breeding. Here, different SSP showed different heritabilities (29.5–76.8%) (Supplementary Table S3). It suggested that four SSP may have differences in generation selection. To some extent, Alb with higher heritability was effectively selected in early generation, but it may not be good for Glu and Total with lower heritability. The relationship between different SSP components in rice is complex. In the present study, there were significant positive correlations between Glu and the other three SSPs in both environments (Supplementary Table S4), that was consistent with

reported study (Zhang et al., 2008). These results indicated that Glu and other three SSPs might partially share the common genetic mechanism.

### GWAS for Seed Storage Protein in Rice

The use of high-density, genome-wide SNPs in GWAS not only can detect candidate genes, but also have a comprehensive understanding of the regulatory mechanism of related traits. Among the 107 associations detected in the study, few sites were repeatedly detected in both environments. Compared with Glo and Glu, Alb and Pro have larger variations in the populations (Supplementary Table S5). Similarly, association mapping was more efficient in identifying associations for Alb and Pro with relatively higher heritability than total SSP and Glu with lower heritability (Supplementary Tables S3–S5). On the other hand, we noticed that all associations detected for Alb and Pro explained more than 100% of the phenotypic variance. This suggests that phenotypic variance explained by the interaction between some of these associations might be very large.

Although genome-wide association studies are becoming more sophisticated, it should be noted that association mapping may lead to false positive associations largely caused by population structure. In order to reduce false-positives resulting from genetic structure, we also analyzed the indica and japonica panels separately by LMM. Thirteen associations were adjacent to previously known grain-quality-related genes (**Table 1** and Supplementary Table S5). These results suggested that association mapping was an effective way to find candidate genes for SSPs in rice. And 28 associations were detected in previously reported intervals or QTLs (Aluko et al., 2004; Wang et al., 2008). QTL pro6 was repeatedly detected for Alb (sf0605251091 and sf0601764762) and Glu (sf0602321094). Lead SNPs (sf0709447538 and sf0712842943 for Glu, sf0705739605 and sf0706196757 for Pro) were co-located in intervals 7–4 and 7–5. Although large numbers of QTLs for grain protein content were detected in the past; only one major QTL has been cloned (Peng et al., 2014). For the lead SNP sf0709447538, explaining 12.7% of the Glu variation and located to the same chromosome region as QTL (7–4 and 7–5), we developed NILs to confirm the QTL. The results showed that it was significantly difference in Glu between two NILs and suggested that the QTL was reliable. To further purify the genetic backgrounds and fine map the QTL, it is needed to backcross the NILs for recombinant screening. It has been reported that a quantitative trait locus (qγ27) affecting expression of 27-kDa γ-zein has been successfully cloned by GWAS and linkage mapping analysis (Liu et al., 2016).

We detected a lead SNP responsible for Alb in the Wx gene region (**Figure 5** and **Table 1**), This region is considered a hotspot of major QTLs (6-2; qPC-6; and pro6) for protein content in rice (Aluko et al., 2004; Wang et al., 2008; Lou et al., 2009; Yu et al., 2009). Analysis of the correlation among Alb, amylose content and quantity of 2.3 kb Wx RNA showed that Wx might influence Alb content, that needs transgenic experiment to verify the function of Wx. Additionally, Starchmetabolism-related genes, such as PUL (Fujita et al., 2009), ISA2 (Utsumi et al., 2011), GBSSII (Hirose and Terao, 2004), Flo4 (Kang et al., 2005), and AGPS2a (Akihiro et al., 2005) were associated with SSPs in different populations. Overexpression of albumin gene RAG2 significantly increased total protein content, prolamin, glutelin, and amylose content, but decreased total starch (Zhou et al., 2016). Carbon and nitrogen metabolisms show a cooperative modification, and consequently, validation the function of these associations might help to better understand the genetic relationship between SSP and starch contents in rice grain.

### Application in the Improvement of Grain Quality in Rice

Grain with good quality could be developed by regulating the SSP content. SSP is a typical quantitative trait typically affected by environment (Shewry, 2007). Combination of conventional breeding and molecular techniques, e.g., markerassisted selection (MAS), may provide a more efficient approach for improving the SSP content of the grain than classical breeding alone (Zhang et al., 2008). The Lgc1 mutation has been used to development new low easy-to-digest protein rice varieties with low glutelin content and high prolamin content, which is useful for patients with chronic renal failure (Mochizuki and Hara, 2000; Nishimura et al., 2005; Morita et al., 2009). O2 mutation could increase free lysine and tryptophan levels by reducing the synthesis of zeins in maize, which is useful for QPM breeding (Mertz et al., 1964; Liu et al., 2016). Genes and possible causative SNPs identified in the present study could be used as potential targets for rice grain nutritional quality improvement. Glutelin is the most easily to digest and contains high lysine (Tanaka et al., 1995). The lead SNP sf0709447538 co-located in QTLs (7–4, 7–5, and 7–6) was validate to have an effect on glutelin, which could be targets for map-based cloning of the candidate genes to illuminate the molecular mechanism of glutelin and improve grain nutritional quality by MAS. RM1 and RP6 are adjacent genes on chromosome 7, and with higher levels expression and higher Pro, hap2 of RP6 and hap3 of RM1 are distributed mainly in the indica subpopulations (**Figure 4** and Supplementary Table S6). Therefore, the region covering RM1 and RP6 from japonica subpopulation can be a promising target for reducing Pro content and achieving better nutritional quality in indica cultivars. Wx may be a key gene for regulating the content of Alb and amylose content. The availability of the Wx gene sequence provides the possibility of improving the protein content via Wx gene modification. Different haplotype combinations of candidate genes for SSP would produce grains with different eating and nutritional quality. Genes and possible causative SNPs identified in the present study could be useful for breeding rice cultivars with favorable eating and nutritional quality.

### CONCLUSION

In the present study many associations were identified for five SSP traits using GWAS in two environments. We detected novel associations, known SSP genes, and SNPs adjacent to known starch-metabolism-related genes. We also analyzed haplotypes of known grain-quality-related genes. Our results suggested that

GWAS was an effective way to identify genes for rice SSP traits and the level of 3.3 kb Wx pre-mRNA is positively correlated with Alb content, providing new insights into the genetic basis of rice quality. Overall, we provided useful information that could be used in future gene functional studies and rice quality improvement.

### AUTHOR CONTRIBUTIONS

YH designed and supervised the experiments. PC, ZS, GL, and BW performed all the phenotypic evaluations. LM, YiL, PC, WD, and BP performed analysis and interpretation of the data. PC wrote the paper. GW, YaL, and DZ provided rice germplasm samples. GG, QZ, JX, and XL participated in the field management and logistic work.

### FUNDING

This work was supported by grants from the from the Ministry of Science and Technology of China (Grants

### REFERENCES

2016YFD0100501, the National Program on R&D of Transgenic Plants (2016ZX08009-003-004 and 2016ZX08001002-002), the National 863 Project (2014AA10A604), and the earmarked fund for the China Agriculture Research System (CARS-01-03) of China.

### ACKNOWLEDGMENTS

We thank Xufeng Bai, Peng Yun, Qiuxiang Luo, Haijiao Dong, Hao Zhou, Pingbo Li, Quanxiu Wang, Dujun Wang, Yuanyuan Zheng, Zhongmin Han, Xiaokai Li, Hu Zhao, and Wei Chen for editing, suggestions, and assistance.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00612/ full#supplementary-material

(Oryza sativa L.) and the function of PUL on starch biosynthesis in the developing rice endosperm. J. Exp. Bot. 60, 1009–1023. doi: 10.1093/jxb/ern349

amylopectin biosynthesis in rice endosperm. Plant Physiol. 156, 61–77. doi: 10.1104/pp.111.173435

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Chen, Shen, Ming, Li, Dan, Lou, Peng, Wu, Li, Zhao, Gao, Zhang, Xiao, Li, Wang and He. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Banana21: From Gene Discovery to Deregulated Golden Bananas

Jean-Yves Paul1 \*, Robert Harding1 , Wilberforce Tushemereirwe2 and James Dale1

1 Centre for Tropical Crops and Biocommodities, Queensland University of Technology, Brisbane, QLD, Australia, 2 National Agricultural Research Organisation, Kampala, Uganda

Uganda is a tropical country with a population in excess of 30 million, >80% of whom live in rural areas. Bananas (Musa spp.) are the staple food of Uganda with the East African Highland banana, a cooking banana, the primary starch source. Unfortunately, these bananas are low in pro-vitamin A (PVA) and iron and, as a result, banana-based diets are low in these micronutrients which results in very high levels of inadequate nutrition. This inadequate nutrition manifests as high levels of vitamin A deficiency, iron deficiency anemia, and stunting in children. A project known as Banana21 commenced in 2005 to alleviate micronutrient deficiencies in Uganda and surrounding countries through the generation of farmer- and consumer-acceptable edible bananas with significantly increased fruit levels of PVA and iron. A genetic modification approach was adopted since bananas are recalcitrant to conventional breeding. In this review, we focus on the PVA-biofortification component of the Banana21 project and describe the proof-of-concept studies conducted in Australia, the transfer of the technology to our Ugandan collaborators, and the successful implementation of the strategy into the field in Uganda. The many challenges encountered and the potential future obstacles to the practical exploitation of PVA-enhanced bananas in Uganda are discussed.

Keywords: East African highland banana, staple crop, Uganda, micronutrient deficiency, vitamin A deficiency, pro-vitamin A, carotenoids, biofortification

### INTRODUCTION

Vitamin A (VA) or retinol is an important nutrient which supports vital physiological and developmental functions. Since it cannot be synthesized de novo, VA has to be acquired through a diversified diet (Fraser and Bramley, 2004; Kimura et al., 2007; Beyer, 2010; Fitzpatrick et al., 2012). Whereas retinol is derived directly from animal sources, plant-derived pro-vitamin A carotenoids (PVACs) such as α- and β-carotene must first be converted into retinol by the body (van den Berg et al., 2000). The majority of populations living in developing countries depend on starchy food staples such as cassava (Manihot esculenta), maize (Zea mays), potato (Solanum tuberosum), rice (Oryza spp.), plantain, and banana (Musa spp.) which unfortunately are largely deficient in essential micronutrients such as PVACs.

Vitamin A deficiency (VAD) causes a number of VAD disorders including night and total blindness, premature death (Sommer and Vyas, 2012), and reduced immunity leading to increased risk of childhood infections and high infant mortality (Herbers, 2003; Bai et al., 2011). VAD affects an estimated 190 million pre-school children worldwide, most of whom live in developing countries, with a reported 5.17 million registered cases of clinical or severe levels of VAD (WHO, 2009). It is also estimated that 250,000–500,000 children become blind due to VAD each year,

### Edited by:

Felipe Klein Ricachenevsky, Universidade Federal de Santa Maria, Brazil

### Reviewed by:

Michael A. Grusak, Children's Nutrition Research Center, United States Hamid Khazaei, University of Saskatchewan, Canada

> \*Correspondence: Jean-Yves Paul<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Plant Nutrition, a section of the journal Frontiers in Plant Science

Received: 13 February 2018 Accepted: 09 April 2018 Published: 26 April 2018

### Citation:

Paul J-Y, Harding R, Tushemereirwe W and Dale J (2018) Banana21: From Gene Discovery to Deregulated Golden Bananas. Front.  In Uganda, 20% of children aged 6 months to 5 years and 19% of women aged 15–49 years suffered from VAD in 2006, predominantly in the low-income demographics within central Ugandan communities, which heavily rely on banana as a staple food (UDHS, 2006).

In the Great Lakes region of East Africa highland bananas (EAHBs) are an important food security crop and the main food staple (Adeniji et al., 2010). In Uganda, it is estimated that 75% of farmers grow bananas contributing to about 7% of global banana and plantain production (Ssebuliba et al., 2006). In addition, Uganda is the largest banana consumer in the world with an estimated per capita consumption of 220–250 kg/year (Tushemereirwe et al., 2006). In rural populations of Uganda where EAHBs form a major and sometime unique part of the diet, 20 µg/g dry weight (dw) β-carotene equivalent (β-CE) is the minimum amount of PVA required in banana fruit to provide 50% of the estimated average requirement (EAR) of VA (Paul et al., 2017). Despite their considerable biodiversity, most EAHB varieties grown in East Africa have low levels of essential micronutrients such as iron (Fe), zinc (Zn), and PVACs (Davey et al., 2009; Fungo et al., 2010). Overreliance on banana in this region, and in Uganda in particular, has contributed to the exacerbation of micronutrient deficiency-related illnesses.

Strategies such as diet diversification, supplementation, and food fortification have been used to help alleviate some of these ailments with varying levels of success (WHO, 2009; Gómez-Galera et al., 2010; Fitzpatrick et al., 2012; Bruins and Kraemer, 2013). However, whereas these interventions are very successful in an urban context where the target population is in the vicinity of service providers, they fail to reach rural communities that are most in need (Dalmiya and Palmer, 2007; Victora et al., 2008; WHO, 2009). As such, in Uganda, VADs in children are relatively low in urban areas but remain elevated in rural communities (UDHS, 2006).

More recently, biofortification has emerged as a complementary, cost-effective, and sustainable approach to deliver micronutrient-dense crops to the poorest and hardestto-reach communities. This can be achieved through either conventional breeding, where the necessary traits are available within the accessible "breeder's gene pool," or through genetic modification. Examples of successful conventional breeding approaches include the biofortification of maize (Egesel et al., 2003; Harjes et al., 2008), sweet potato (Mwanga et al., 2009), and cassava (Ceballos et al., 2012). These biofortified products are already being disseminated in various parts of the world including Uganda (Anderson et al., 2007; Hotz et al., 2012; Talsma et al., 2016). Arguably the most successful example of biofortification through genetic engineering is the development of Golden Rice (GR) (Ye et al., 2000; Paine et al., 2005). Considering the popularity of EAHBs in East Africa, biofortification of this food staple with enhanced levels of PVACs (or other micronutrients) is now believed to be the best long-term, sustainable, and cost-effective strategy to ease the burden of VAD in high risk populations of Uganda. The use of conventional breeding to develop PVA-biofortified EAHBs is constrained by their low fertility and the lack of high PVA EAHB varieties in the known gene pool. Further, any new varieties developed are unlikely to possess the attributes of locally preferred landraces of EAHBs. Genetic modification is therefore the fastest and most reliable approach to improve the existing preferred varieties.

Here, we review the foundations, goals, achievements, future prospects, and challenges of Banana211 , a project that undertook the challenges of alleviating micronutrient deficiencies in Uganda by enhancing the nutritional content of its staple food, EAHBs, through genetic modification.

### THE PROMISE OF BANANA21 – MORE NUTRITIOUS BANANAS MADE IN UGANDA, BY UGANDAN SCIENTISTS, FOR THE UGANDAN PEOPLE

Banana21 is one of four original projects funded by the Grand Challenges in Global Health (GCGH) program of the Bill and Melinda Gates Foundation (BMGF) to "Create a Full Range of Optimal, Bioavailable Nutrients in a Single Staple Plant Species." Banana21 is a collaborative research project between the Centre for Tropical Crops and Biocommodities at Queensland University of Technology (QUT) in Australia and the National Banana Research Program of the National Agricultural Research Organisation (NARO) of Uganda. The project aims to help "Alleviate VAD and iron deficiency anemia through the micronutrient enhancement of the staple food of Uganda, bananas."

From the outset, it was recognized that the success of the project depended on its objectives being widely understood and accepted in Uganda. As such, we considered it extremely important to train and empower young Ugandan scientists to deliver the project milestones thus ensuring the biofortified genetically modified (GM) bananas were generated in Uganda, by Ugandan scientists, from Ugandan varieties, and for the benefit of the Ugandan people. As a consequence, in addition to scientific discoveries and their applications, Banana21 has been a capacity building project where technology transfer continues to be at the core of every phase, milestone, and decision-making activity.

An overview of the Banana21 project strategy is shown in **Figure 1**. Initially, the proof-of-concept research was developed at QUT using the locally grown dessert banana cultivars "Cavendish." The QUT components included gene and promoter discovery, tissue culture and transformation, field trials of GM bananas, and downstream fruit sampling and analysis. The development and implementation of a comprehensive stewardship plan was also a key project component. Following proof-of-concept in Australia, the technology and know-how was to be continuously transferred to NARO in Uganda to generate PVA-biofortified local EAHB varieties. Further, it was necessary for the infrastructure to be implemented at NARO to ensure that the GM plants produced would be tested in the laboratory and the field with the consistency and rigor necessary to generate the data required for a GM product deregulation dossier.

1http://www.banana21.org/

### THE CHALLENGES AND ACHIEVEMENTS OF BANANA21 – TRANSGENIC BIOFORTIFIED EAHBs WITH SIGNIFICANTLY HIGHER FRUIT PVA CONTENT

### Phase 1 – Early Discovery and Technology Transfer

HDR, higher degree research.

At the commencement of the project in 2005, the only practical demonstration of PVA biofortification in a staple food crop was GR. The initial GR strategy involved the re-engineering of a carotenoid biosynthesis pathway in normally carotenoid-free rice endosperm by the endosperm-specific [glutelin 1 (gt1) promoter] expression of a daffodil-derived phytoene synthase (psy) transgene in combination with the constitutive (CaMV 35S) expression of a bacterial phytoene desaturase (crtI) gene. Although successful, the level of carotenoid accumulation in rice endosperm was still considered to be too low for practical exploitation (Ye et al., 2000). Subsequent research revealed that the origin of the psy transgene and choice of promoter were important factors affecting PVAC accumulation levels (Paine et al., 2005). This led to the development of Golden Rice 2 (GR2) whereby the use of a maize-derived psy transgene and bacterial crtI, both under the control of the gt1 promoter, resulted in high levels of PVAC accumulation (Paine et al., 2005). Based on this success, the GR2 strategy was deemed to be the most logical approach to develop VA-biofortified bananas. One of the major initial challenges was the identification of suitable transgenes and promoters for expression in banana. As such, the early phase of the project focused on designing and testing large numbers of expression constructs containing a suite of different promoters and transgene combinations.

Due to the paucity of information regarding transgene expression in banana fruit, it was necessary to assess a range of different promoters for their spatiotemporal activity in banana, particularly fruit. As such, several constitutive promoters were isolated from different sources in addition to promoters that controlled the expression of genes involved in banana fruit development. These were fused to the β-glucuronidase reporter gene (uidA), transformed into banana embryogenic cell suspensions (ECS) (Khanna et al., 2004) and transgenic plantlets regenerated. A list of all promoters tested, their origin, and specificity is presented in **Table 1**.

In combination/parallel with the promoter characterization study, a range of different transgenes were also assessed. Initially, the strategy used to develop GR2 was adapted to banana and the first generation of expression constructs was made to express the maize phytoene synthase 1 (ZmPsy1B73) transgene alone or in combination with the bacterial (Pantoea ananatis) carotene desaturase transgene (PaCrtI) and controlled by various promoters (**Table 1**). Our efforts also focused on isolating and using PVA-associated cisgenes from banana in the hope of minimizing gene silencing and thus achieving more stable gene expression over several generation. The use of banana-derived cisgenes/s was also considered to be more advantageous from a future deregulation perspective. The biodiversity of banana and plantain is huge especially in Southeast Asia and Papua New Guinea where bananas are believed to have originated (Perrier et al., 2011). As such, the ability of various cultivar to accumulate PVACs in their fruit is also very diverse (Davey et al., 2009; Ekesa et al., 2015). Of particularly interest to us was a small group of bananas originating from the Pacific called the Fe'i bananas (Musa troglodytarum) which accumulate extremely high levels of fruit PVACs (Englberger et al., 2003). One such Fe'i banana variety, "Asupina," became our model cultivar for not only understanding carotenoid metabolism and accumulation in banana fruit but also as the source of a cisgene encoding a phytoene synthase, MtPsy2a (Mlalazi et al., 2012; Buah et al., 2016). Following the cloning and molecular characterization of MtPsy2a, expression cassettes containing this gene with and without PaCrtI and controlled by various promoters were also constructed and used to generate transgenic "Cavendish" banana lines in order to assess the levels of PVA accumulation in fruit.

During Phase 1, an important training program was also initiated whereby several students from Uganda and Kenya commenced their higher degree research (HDR) doctoral studies at QUT. The research projects developed for these students focused primarily around banana tissue culture and transformation technologies so that the knowledge and experience gained in the QUT laboratories would be transferred back to Africa for the benefit of the project. Simultaneously,

TABLE 1 | List of promoters and genes tested in AFT-1 of the Banana21 project.

training of technical staff and capacity building around infrastructure and laboratory equipment began at NARO, Uganda. A major component of this training program was the development of protocols and standard operating procedure (SOP) for generating and transforming ECS of local banana varieties. From the beginning, the Ugandan component of the project focused on establishing ECSs from three popular Ugandan banana cultivars, namely "Nakitembe," "M9," and "Sukali Ndiizi" since these were considered the most appropriate target cultivars for biofortification. True EAHB varieties such as "Nakitembe" are the preferred cultural choice among Ugandans and are usually grown in the highlands where disease pressure is minimal. The disease-resistant EAHB hybrid "M9" is less popular but was chosen as it is more productive in the lowlands where diseases such as black Sigatoka are a serious limitation to banana production. "Sukali Ndiizi" was chosen because it is a small sweet banana which is popular among children, the most vulnerable target population.

### Phase 2 – Proof of Concept, Field Trials in Australia and Uganda

The banana fruit PVAC target levels necessary to deliver 50% of the EAR of VA in vulnerable populations was estimated at 20 µg/g dw β-CE. The PVA-biofortification proof-of-concept research was done in Australia using transgenic "Cavendish" and "Lady Finger" bananas as the models with the aim of obtaining the target fruit PVA levels with no changes in agronomical characteristics. One of the major limitations of the project was the lengthy timeframe from transformation of banana ECS through to fruit harvest which is approximately 30 months. This limitation precluded the serial testing of our expression constructs. Consequently, our strategy was to transform "Cavendish" and "Lady Finger" ECS in Australia with a large number of different expression constructs, regenerate transformants, and field trial all the plants in parallel. At the end of Phase 1, between 10 and 30 independents transgenic "Cavendish" banana lines for each of the expression constructs had been generated. Considering the very large number of lines to be tested, and the fact that bananas are a very large crop, it was not practical to conduct a trial in the glasshouse. With these limitations in mind, we assessed these lines directly in the field without prior glasshouse characterization. The first Australian field trial (AFT-1) of GM "Cavendish" banana lines was subsequently established in 2009. This trial contained all available independent lines for each construct but only a single plant per line. A total of 28 constructs were tested, 14 to test promoter activity (as promoter/uidA reporter gene fusions) and 14 to test the same promoters in combinations with three transgenes (MtPsy2a, ZmPsy1B73, and PaCrtI) encoding proteins involved in the biosynthesis of PVACs (**Table 1**). This was the first GM banana field trial in Australia and was conducted under a license issued by the Australian Office of the Gene Technology Regulator (OGTR)2 .

Of the 14 promoter/uidA fusion combinations tested in the field, three promoters (Ubi, Exp1, and Aco) consistently conferred the strongest levels of GUS expression. The constitutive maize polyubiquitin (Ubi) promoter showed consistently strong activity from the earliest stages of fruit development through to maturity. In contrast, the banana fruitspecific expansin 1 (Exp1) and ACC oxidase (Aco) promoters were only active in the later stages of fruit development (Paul et al., 2017). When we subsequently investigated the accumulation of PVACs during fruit development, a similar trend was seen whereby constitutive expression of either

2http://www.ogtr.gov.au/internet/ogtr/publishing.nsf/Content/dir076rarmp-toc~ dir076rarmp-ch1~dir076rarmp-ch1s5

MtPsy2a or ZmPsy1B73 using the Ubi promoter increased PVA accumulation from the earliest stages of fruit development while PVAC accumulation was restricted to the late stages of fruit development when the same transgenes were tested under the control of either the Exp1 or Aco promoters. Analysis of fruit samples from plants transformed with MtPsy2a, ZmPsy1B73, and/or PaCrtI revealed that constitutive (Ubi) expression of MtPsy2a resulted in the highest fruit β-CE levels which reached almost 19 µg/g dw in the plant crop (**Figure 2A**). Although this was just below the target level, it nonetheless represented a 11-fold increase over wild-type banana PVA baseline levels and was highly encouraging (Paul et al., 2017).

A second Australian field trial (AFT-2) was planted in September 2012 to test the efficacy of new PVA biofortification constructs as well as to monitor transgene stability in multiple plants of selected lines from AFT-1 through three successive generations. The additional PVA genes and promoters tested in this second field trial (**Table 2**) were tested in combination with the MtPsy2a gene under the control of various promoters. Unfortunately, none of the new combinations tested were found to elevate fruit PVA levels as high as those achieved by the constitutive expression of MtPsy2a alone in AFT-1. However, an extremely interesting and important outcome from AFT-2 was the observation that fruit PVA levels increased over successive generations in some lines to levels exceeding the target. For example, the fruit PVA levels of a line containing Ubi-MtPsy2a was found to increase from 11.7 µg/g dw β-CE in the first generation to 75.1 µg/g dw β-CE in the fourth generation. Qualitative and quantitative phenotypic and agronomical data have been recorded for all the lines tested and has revealed some interesting trends (Paul et al., 2017). The origin of the transgene had the biggest impact with, for example, banana lines expressing the banana MtPsy2a appearing normal with no variation in critical agronomical features such as yield and cycle time. In contrast, although expression of the maize Psy1B73 gene resulted in increased fruit PVA levels, many lines with undesirable phenotypes such as stunting and photo bleaching were observed.

During Phase 2, scientists at NARO had established regenerable ECSs, as well as efficient transformation protocols, for two varieties, "M9" and "Sukali Ndiizi." These breakthroughs allowed the generation of transgenic banana plants from both varieties containing either the ZmPsy1B73 or the "Asupina" derived phytoene synthase 2a (MtPsy2a) transgenes. Since the results of AFT-1 were not available when this activity commenced, these gene were placed under the control of the banana-derived Exp1 promoter (Mbabazi, 2015). The transgenic lines generated from both cultivars were subsequently grown in a field trial at NARL-NARO, Kawanda Uganda in 2010 to assess the fruit PVAC levels. This first Ugandan field trial (UFT-1), conducted under authorization from the Ugandan National Biosafety Committee (NBC), was a very important milestone for Banana21 as it represented the first GM crop field trial in sub-Saharan Africa where the events had been created by local scientists from a national laboratory.

At the completion of Phase 2, proof-of-concept for PVA biofortification had been demonstrated in Australia using "Cavendish" and "Lady Finger" bananas and revealed that target fruit PVA levels could be achieved using the constitutive (Ubi) or fruit-specific (Aco) expression of MtPsy2a alone (Paul et al., 2017). Further, very few off-type traits, such as reduced yield and increased cycle time, were observed in the transgenic lines and the trait appeared stable over successive generations. Importantly, the field trial in Uganda (UFT-1) produced multiple lines with fruit PVA levels higher than their respective controls including one Exp1-MtPsy2a line of M9 with 33.1 µg/g dw β-CE (Mbabazi, 2015). With this knowledge, a new generation of plant expression vectors were made at QUT and two constructs containing Ubi-MtPsy2a and Aco-MtPsy2a

FIGURE 2 | (A) PVA-biofortified "Cavendish" banana in Australia and (B) PVA-biofortified (left) and wild-type control (right) "Nakitembe" EAHB in Uganda.

TABLE 2 | List of new promoters and genes tested in AFT-2 of the Banana21 project.

were transferred to NARO for transformation into EAHB varieties.

### Phase 3 – Product Development: Early Events Selection Field Trial in Uganda

This phase of the project (2012–2017) initially involved generating a total of 200 independent transgenic lines each of EAHB cultivars "M9" and "Nakitembe" containing MtPsy2a under the control of either the Ubi or the Aco promoters. Following molecular characterization by the now fully trained technical staff and highly qualified scientists at NARO, the transgenic lines were field planted (UFT-2) at NARL-NARO, Kawanda in August 2014 to identify suitable elite lines for further testing in multi-location field trials (MLTs). During those 5 years, the NARO team showed exceptional professionalism in learning and implementing good practices around generating, handling, and tracking GM products. The results from UFT-2, which will soon be published, clearly demonstrate that fruit PVACs can consistently accumulate at levels above the required target of 20 µg/g dw β-CE in GM-biofortified "M9" and "Nakitembe" without phenotypic alteration of the plants (**Figure 2B**).

From the current transgenic line selection trial at Kawanda, 10 elite lines each of "M9" and "Nakitembe" will be selected to progress through to future MLTs. The initial selection process began in 2017 by selecting lines with fruit PVA levels equal or greater than 20 µg/g dw β-CE at the full green developmental stage (harvesting stage in Uganda) and yield within 20% of non-GM controls plants. From this initial selection, molecular analysis was used to identify lines containing fewer than three copies of the integrated expression cassettes with a preference for single integrations. Although selection of "single copy" events is preferred in seed crops to produce homozygous lines that do not segregate for the transgenic trait in future generations, in a vegetatively propagated crop such as banana, it is only rationalized to increase the likelihood of events with "clean insert."

During this phase, NARO scientists also conducted a preliminary blind sensory panel test of traditionally prepared banana meal (matooke) using fruit from both non-GM and PVA-biofortified "M9" and "Nakitembe" to compare appearance and texture but not taste. Interestingly, 80% of the panelists (n = 15) rated fruit from one of the "M9" PVA-biofortified lines as the most preferred whereas fruit from its non-GM counterpart had the highest dislike proportion (29%).

On the basis of these encouraging results, Banana21 entered its fourth and final phase of funding in October 2017. This final phase is focused on generating all the data required for the deregulation of a GM PVA-biofortified EAHB in Uganda. The prospect of the release of the world's first deregulated GM banana developed from an African laboratory is as much exciting as it is daunting.

### THE FUTURE CHALLENGES OF BANANA21 – THE NOT SO LONG ROAD AHEAD TO DEREGULATION

From a technology perspective, the groundwork in Uganda has been completed and the exceptional results from UFT-2 have provided the ideal platform to further select elite lines to progress through to MLTs and ultimately to farmer release.

Prior to MLTs, critical data necessary for the compilation of a deregulation dossier must be obtained. Whole-genome sequencing of the pre-selected lines is necessary to allow the site(s) of transgene insertion into the host genome to be identified. This is necessary to ensure that (i) the transgene is intact, (ii) there are no new open-reading frames created, (iii) there is no disruption of endogenous open-reading frames, and (iv) no plasmid sequence has been integrated into the host genome. Finally, it is essential that the composition of the fruit from the deregulated lines is very similar to that derived from the wild type. Therefore, compositional analysis for food characteristics such as calories, calories from fat, carbohydrates, protein, ash, and moisture will be done, and lines with values >15% different to controls will be discarded.

The next phase of the project involves the identification of two lines (one lead and one reserve) each of "M9" and "Nakitembe" from the MLTs, and obtaining the necessary agronomic, biochemical, and molecular data to ultimately prepare the dossier for submission to regulators in Uganda for deregulation and general release. At the end of this phase of the project (December 2021), these four lines are expected to meet all the agronomic, biochemical, and molecular analyses and biosafety assessment required for deregulation in Uganda under the proposed Biosafety Bill.

One of the major hurdles for Banana21 is the lack of a regulatory framework for biotechnology in Uganda. Without a regulatory body controlling the safe application

of biotechnologies, Banana21 will not be able to release PVA-biofortified "Golden bananas" to farmers and consumers in need. After ratifying the Cartagena Protocol on Biosafety in 2002, it took 6 years for Uganda to approve a policy on Biotechnology and Biosafety in 2008. Under the current policy, the NBC supervises all GMO activities up to the stage of Confined Field Trial (CFT) under the supervision of the National Council of Science and Technology (UNCST) Act 1990. Therefore in 2012, the National Biotechnology and Biosafety Bill, 2012, was introduced into parliament to provide a regulatory framework and guide the implementation of modern biotechnology in Uganda to minimize any potential risks to the environment, human, and animal health. After 5 years, Uganda's Parliament passed the National Biotechnology and Biosafety Bill into law becoming the Biosafety Act, 2017 on October 4, 2017. The new law will not only benefit Banana21 but also a multitude of other biotechnology products developed in Uganda such as bananas with bacterial wilt resistance, drought tolerant maize, bollworm resistance and herbicide tolerant cotton, and new cassava varieties with resistance to Cassava mosaic and brown streak viruses.

Our initial target of 20 µg/g dw β-CE was calculated using a bioconversion factor of PVACs to retinol of 6:1 based on the results of a study using Mongolian gerbils (Bresnahan et al., 2012). Underestimating this bioconversion factor would raise the target above its current value. For this reason, Banana21 with financial help form HarvestPlus and the BMGF, commissioned a nutrition study at Iowa State University to determine a more accurate bioconversion ratio in humans. Although the results from this study are not yet available, some of the PVA-biofortified lines that have been developed under Banana21 have over four-times the initial target value with no yield penalty and could potentially be substituted if all other deregulation criteria are met.

The development of robust diagnostics for banana viruses was another important component of this project. These diagnostics form the basis of a banana virus indexing protocol that has been rolled out to ensure that the plantlets derived from the Banana21 project are virus tested, thus reducing the potential distribution of infected planting material to farmers.

Since the majority of banana growers in Uganda are subsistence farmers, the distribution strategy adopted in the future will need to minimize the cost of planting material while maximizing the rate of distribution. Initial propagation of the lines to be released will be done at NARO and the plantlets will be sent to small banana micropropagation laboratories and also used to establish small "mother gardens" for the initial production of suckers. A cost-effective and self-sustaining strategy for dissemination will then involve identifying "innovative farmers" that will be given suckers. For every sucker, they will be asked to give away two suckers to neighbors who in turn will be asked to give away two suckers for each one received under the scheme.

### REFERENCES

Adeniji, T. A., Tenkouano, A., Ezurike, J. N., Ariyo, C. O., and Vroh-Bi, I. (2010). Value-adding post harvest processing of cooking bananas

A key component of the next phase of the project will be the implementation of a comprehensive stewardship and communication plan. This includes (i) forming a Technical Advisory Committee (TAC) that meets regularly and provides scientific, strategic, and biosafety expertise, (ii) implementing and regularly updating SOPs, (iii) keeping accurate and safe records of the data with tools such as the BananaTracker software developed by QUT, and (iv) meetings with the Australian OGTR and Food Standards Australia and New Zealand (FSANZ) to seek advice on the requirements for deregulation if these lines were to be deregulated in Australia. The NARO team has also been involved in various communication activities in an attempt to educate the public and de-mystify the use of GMOs. Important stakeholders are targeted through workshops and information sessions, as well as various paper-based and audio-visual communication materials.

### FINAL REMARKS AND CONCLUSION

From the outset in 2005, Banana21 has been on a trajectory to develop lines of EAHBs with levels of fruit PVACs that would provide 50% of the EAR of VA with consumption of only 300 g per person per day. Based on the significant progress thus far, it is highly likely that the transgenic lines developed under Banana21 will be released by 2021 and have a significant impact in alleviating VAD in a sustainable way, especially in rural Uganda where bananas are a fundamental part of the culture. The PVA-enhanced, disease-resistant "M9" line will have the greatest impact in lower elevations of Uganda where the disease pressure is high, while the PVA-enhanced "Nakitembe" line will have greatest impact in the highlands where there is much lower disease pressure. The importance of banana as a food security crop (perennial nature, year-round production, and ability to cope with long periods of drought) associated with a low cost, farmer-driven distribution strategy should ultimately see "Golden bananas" adopted as a widespread and efficient VAD alleviating strategy in the next decade.

### AUTHOR CONTRIBUTIONS

J-YP and RH drafted the initial manuscript while WT and JD reviewed and provided the constructive criticisms.

### FUNDING

The Banana21 team project was, and still is, supported by a Grant from the Bill & Melinda Gates Foundation and the Department for International Development (United Kingdom).

Anderson, P., Kapinga, R., Zhang, D., and Hermann, M. (2007). "Vitamin A for Africa (VITAA): an entry point for promoting orange-fleshed sweetpotato to

(Musa spp. AAB and ABB genome groups). Afr. J. Biotechnol. 9, 9135–9141.

combat vitamin A-deficiency in sub-Saharan Africa," in Proceedings of the 13th ISTRC Symposium, Apartado, 711–720.

maize for plant breeding trials. Food Chem. 100, 1734–1746. doi: 10.1016/j. foodchem.2005.10.020

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Paul, Harding, Tushemereirwe and Dale. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

fpls-09-00558 April 24, 2018 Time: 17:15 # 8

# Improving the Yield and Nutritional Quality of Forage Crops

Nicola M. Capstaff\* and Anthony J. Miller\*

John Innes Centre, Norwich, United Kingdom

Despite being some of the most important crops globally, there has been limited research on forages when compared with cereals, fruits, and vegetables. This review summarizes the literature highlighting the significance of forage crops, the current improvements and some of future directions for improving yield and nutritional quality. We make the point that the knowledge obtained from model plant and grain crops can be applied to forage crops. The timely development of genomics and bioinformatics together with genome editing techniques offer great scope to improve forage crops. Given the social, environmental and economic importance of forage across the globe and especially in poorer countries, this opportunity has enormous potential to improve food security and political stability.

#### Edited by:

Felipe Klein Ricachenevsky, Universidade Federal de Santa Maria, Brazil

### Reviewed by:

Zhipeng Liu, Lanzhou University, China Giovanna Attene, University of Sassari, Italy

#### \*Correspondence:

Nicola M. Capstaff<EMAIL_ADDRESS>Anthony J. Miller<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Plant Nutrition, a section of the journal Frontiers in Plant Science

Received: 06 December 2017 Accepted: 06 April 2018 Published: 24 April 2018

#### Citation:

Capstaff NM and Miller AJ (2018) Improving the Yield and Nutritional Quality of Forage Crops. Front. Plant Sci. 9:535. doi: 10.3389/fpls.2018.00535 Keywords: forage, nutritional enhancement, grass production, legumes, breeding, management

## INTRODUCTION

Forage grasslands are used to feed livestock and globally it has been estimated that they represent 26% of the land area, and 70% of agricultural area (FAO, 2010). Such crops are significant economically, as the European example shows (see **Figure 1**). Forage crops are usually grasses (Poaceae) or herbaceous legumes (Fabaceae). Some tree legumes such as mulga (Acacia aneura) and leadtree (Leucaena leucocephala) are also grown in desert and tropical grasslands (Muir et al., 2011). In the tropics, popular grasses include Napier grass (Pennisetum purpureum), Brachiaria, and Panicum species. In the poorest parts of the world livestock production is critically important for smallholders' livelihoods. Sub-Saharan Africa is an example and frequently women maintain the livestock production systems (Njuki and Sanginga, 2013). In temperate climates, the main grasses include bentgrass (Agrostis spp.), fescue (Festuca spp.), ryegrass (Lolium spp.) and orchard grass (Dactylis spp.) or hybrids of these. For example, Festuca and Lolium hybrids has been developed from 1970s (Ghesquière et al., 2010) giving rise to crops such as Festulolium pabulare which combines the superior forage quality of Lolium multiflorum with the persistence and stress tolerance of Festuca arundinacea. Some maize (Zea mays) cultivars have been specifically bred for forage. The commonly cultivated herbaceous legumes are trefoil (Lotus corniculatus), medics (Medicago spp.), clover (Trifolium spp.) and vetches (Vicia spp.). Brassica forage species include cultivars of oilseed rape (Brassica napus) and kale (Brassica oleracea). Fodder beet (Beta vulgaris) is another temperate forage. The combination of forage crops grown in any country varies depending on climate and livestock needs, however, the perennial legume lucerne or alfalfa (Medicago sativa) is the most widely cultivated as it can be grown with both temperate and tropical grasses, or as a standalone crop. This is a huge topic to review as there are so many species grown across the world, therefore we have chosen to focus on a few examples, the tropical grasses Pennisetum and Brachiaria, and more prominently the temperate crops Lolium and alfalfa.

In an ideal world, we would all eat pulses rather than the animal products generated from them, as grain legumes are the food that offers the most sustainable future (Foyer et al., 2016). There is continued pressure from many groups to lower human consumption of animal products due to livestock efficiency issues and for human health (Cramer et al., 2017). There is a lack of reliable statistics for the proportion of adults adopting a plant-based diet, but it is estimated to be between 1 and 10% of the population in developed western countries such as within the European and United States (Mcevoy and Woodside, 2010) and studies support these diets as healthy and nutritionally adequate (American Dietetic Association and Dietitians of Canada, 2003). However, the consumption of livestock products can be regarded as important to a healthy diet due to their high nutrient density (CAST, 2013) regardless of the numerous efficiency and environment concerns (Di Paola et al., 2017), particularly true in developing countries where undernourishment incidences are estimated as ∼4–22% of the population (Alexandratos et al., 2006). Livestock production can convert non-edible crops such as the forages into human food, with sustainable intensification possible when inputs and outputs of the system are balanced (Derner et al., 2017).

Moreover, the cultural and social significance of livestock cannot be underestimated and the trend of increased global production is set to continue (Thornton, 2010). Livestock feature prominently across all cultures both in cuisine, but also music and literature. Additionally, in many developing countries the rearing of livestock such as cattle and goats are vital in times of hardship; many view animals as living 'piggy-banks,' that can for example pay the family school fees (Herrero et al., 2013). Therefore, in practice livestock production is set to continue throughout the world and forage crops will be grown for coming decades. Plant research has chiefly focussed on grain crops, but here we argue that there is enormous potential for improving forages. Improving the yield and nutritional quality of forage crops can help mitigate the unsustainable negative impacts of livestock production.

### FORAGE CROPS IN LIVESTOCK DIETS

Forage crops can be feed directly to livestock or can be processed by partial drying or pre-digestion. Because of this processing, animal feeds can be categorized as either bulky feeds or concentrates. Bulky feeds are also termed forage and are produced from grass, cereal and legume cropping as described above, such as alfalfa, Lolium or a mixture of the two. This forage can be provided to animals directly through grazing pasture land or in a processed form, such as hay (where water content is >15%) or dried (pelleted) biomass. Concentrates are generally cereal, oilseed and legumes seeds, or bi-products of their preparation for human food, biofuel and textile. They can also include high energy feedstuffs such as sugar-rich crop molasses and fats of animal origin, for example fish by-catch discards. In industrialized countries, production of both these categories of feed can surpass the amount produced for plant-based food for human consumption; in United States over double dry matter per-capita per year (DM cap/yr.) is produced for animal feed than for foodstuffs (Krausmann et al., 2008).

Livestock diet can therefore be exclusively forage or largely forage with concentrate supplementation. Concentrate supplementation is used to compensate nutritional deficiencies in the forage supply, increase animal performance such as milk production or at particularly challenging periods of development, for example calving. Due to most livestock diet being of forage this review focuses on the main crops used worldwide and will not discuss concentrates. The amount plant science has contributed to improvements in concentrates has been underappreciated and undervalued in literature, however, the role these crops have on livestock production has been reviewed previously (Erb et al., 2012).

Forage crops can be grown in mixed species cultivation to provide nutritional and environmental benefits. By offering livestock mixed grazing pastures or blending feeds, nutritional quality can be enhanced. For example, alfalfa is the highest-yielding perennial forage legume and produces more protein per unit area than other forage legumes and so can be grown alone or in combination with a range of different grass species. Well-managed alfalfa is normally grown successively for 3 or more years, but if harvested too late in the season the crop cannot survive the winter (Bélanger et al., 2006).

### FORAGE NUTRITIONAL CONTENT

### Digestibility

The nutritional status of a forage crop depends upon the concentration (and ratios) of carbohydrates, proteins, and lipids. The composition of these organic nutrients determines the digestibility (D-value) of each crop which along with mineral and vitamins provides the amount of energy which can be derived by the animal (ME measured in MJ/kg DM) (Osbourn, 1980). Such calculations are becoming increasingly prevalent when growers are deciding which crop to grow based and particularly dependent on if the animal is non-ruminant or ruminant.

In forage crops 50–80% of DM is carbohydrate; if this percentage is too low then supplements of grains can be added. The primary types of carbohydrate are the insoluble structural saccharides cellulose and hemicellulose, or the storage forms such as starch and water-soluble polymers (e.g., fructans). These are degraded into simple sugars through cleavage of glycosidic bonds, either by the animal itself (non-ruminant and ruminants) or via microbial digestion and subsequent animal absorption (ruminants only). Different ratios of carbohydrates within the forage crop will have altered downstream digestibility for the animal, especially if the cell-wall structure constrains digestion by the microbial population or limits plant cell wall penetration (Weimer, 1996). Although lignin, a polyphenolic compound within forage, is not a carbohydrate, it has a dramatic impact on the digestibility of cellulose hemicellulose; lignin binds with structural carbohydrates and cell wall proteins and reduces nutrient availability. For forages increased lignin concentration in the growing crop will increase the percentage of indigestible DM. Of the major forage crops grown globally grasses, particularly Lolium perenne, have high digestibility due to high soluble sugar content alongside low lignin content (Ruckle et al., 2017).

Animal digestion of simple carbohydrates produces monosaccharides which can be readily metabolized. In ruminants, only microbial digestion of structural carbohydrates produces simple sugars which are subsequently metabolized to pyruvate. Pyruvate is absorbed by the animal and is metabolized further into volatile fatty acids (VFAs) which are a major energy source, (Bergman, 1990). Ruminants absorb VFAs in their rumen, and the rate of this is dependent on the concentration of individual VFAs, rumen pH and the absorptive area in the ruminal lining.

### Protein

Nitrogen (N) availability to animals is predominantly from forage proteins and are estimated using crude total protein Kjeldahl measurements. Protein is usually abundant in the major form of Ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO), although relative amounts vary between species (Wallace et al., 1997). This is especially true when comparing content in grasses with herbaceous legumes, with red clover (Trifolium pratense), white clover (Trifolium repens) and lucerne (Medicago sativa) grown widely due to their high protein value (Ruckle et al., 2017). Again, lignin will severely affect the digestibility of protein. Some micronutrients like proanthocyanidins or condensed tannins also change the digestibility of protein, but they inhibit protein degradation through binding. This can be advantageous as rapid protein degradation is causative of bloat, however, too high a tannin content will mean protein passing through the digestive track is unabsorbed and therefore a loss in nutrition value (Lees, 1992; Piluzza et al., 2014). This means there is a balance between reduced bloat and animal productivity (Mueller-Harvey, 2006). All grasses contain little or no proanthocyanidins, whereas many legumes especially big trefoil (Lotus pedunculatus) and Sericea lespedeza (Lespedeza cuneata) can have levels as high as 18% DM (Barry and Manley, 1984; Mueller-Harvey, 2006). Other N-containing compounds can be found in forage such as nucleic acids, nitrate and ammonia (Wallace et al., 1997).

### Lipids