To apply this filter it is necessary to indicate the wavelength segment for analysis, the cutting criteria (Group or Individual), the maximum percentage of variations accepted (%), and the number of neighbors (N size). The selection of each wavelength segment and the criteria for each one (% and N Size) will depend on the user experience and the environmental conditions where the measurements were taken; for example, noise at 1800– 1950 nm and 2350–2500 nm is usually wider and stronger than at 1300–1400 nm, so the criteria should consider higher values of % and N Size for the first two segments. In this operational example, the filter was applied to the whole spectral range (350–2500 nm) considering a group filter, with five wavelengths as N size and a maximum accepted variation among them of 20%. **Figure 2** shows the results prior to (A) and after (B) the filter was applied. In this case, the filter was able to detect two main noise zones from 1833 to 1935 nm and from 2422 to 2500 nm.

### Scan Analysis

The Scan analysis module allows detection of abnormal variations among samples within the same scan. In this operational example, the Scan analysis was applied using the

function Run from the current data, that is to say, considering the results obtained using previous filter (without spectral noise). The Maximum variation coefficient was set at 0.5%. The software was able to select 383 scans or plots without problems and 17 where the threshold was exceeded (5, 26, 36, 112–113, 119, 144, 181, 223, 233, 274, 348, 356–358, 395, and 399). In the **Figure 3**, scan or plot 399 is graphed, and the first sample (1195 on red) was selected for deletion. This result could be an indicator of a measurement problem associated with the operator (modification of the measurement angle) or external conditions (e.g., wind speed) during the first sample integrations. In case of all samples from a specific scan need to be deleted, the software will maintain this scan as empty rows, avoiding problems in further analyses.

### Outlier Analysis

This module is a simple and exploratory analysis to identify outlier scans, allowing the user to detect field measurement problems (e.g., calibration). Four scatterplot graphs will show the relationship between any SRI available on the software database and the loaded response variables. If a problem is detected, it is possible to use the Edit option to manually remove the samples. In this operational example, the relationship between NDVI and Yield was used to inspect the possible outlier samples. **Figure 4A** shows how different SRIs (NDVI, SR, PRI, and WI) can generate different data distributions, helping the user in cases where problems are not so evident; on the top left graph (NDVI vs. Yield) two clouds of data points can be identified, divided at the NDVI value of 0.31.

Once the information from the smaller data cloud was analyzed (NDVI < 0.31), it was evident that the data set corresponded to 96 contiguous scans or plots (104–200), suggesting that there were problems associated with the measurement. When information from the spectrometer was checked, it was concluded that the operator had skipped one calibration. It is always important to check the pertinence of negative SRI values because they are probably related to measurement errors.

After identification of the origin of a particular problem, any graph can be selected for editing. In this example (NDVI vs. Yield), the scan with the problem can be selected (**Figure 4B**) and deleted (**Figure 4C**).

### Collinearity Analysis

Using linear regression or ANN, the collinearity analysis module identifies wavelengths that delivering the same predictive information for a given response variable, keeping only those that best explain it. In this operational example, collinearity analysis was applied by considering the results obtained from the scan analysis (Run from current data), with Yield being the response variable in the linear regression (R square cutoff = 0.95). Results of this analysis found 131 wavelengths without collinearity (**Figure 5**).

### Individual Wavelength Analysis

In this module it is possible to assess the relationship of individual wavelengths and a given response variable. Three regression models were selected (Polynomial 1 and 2, and Exponential) to search for wavelengths with determination coefficients above 0.3 in relation to Yield (**Figure 6A**). If the user selects Plot all results, a graph will show the wavelengths below and above the determination coefficient cutoff (**Figure 6B**). These results can be exported to a spreadsheet; for each selected regression model, only wavelengths where the chosen statistic surpassed the cutoff will be shown (**Figure 6C**). In this operational example, there were three groups of wavelengths with determination coefficients above the threshold: 733–1139, 1409–1815, and 1936–2421 nm.

### Spectral Reflectance Index (SRI) Analysis

As in the previous module, when different regression models were considered, SRI analysis has the option to evaluate the relationship between all loaded response variables and all SRIs available in the software database. For this example, three regression models were selected (Polynomial 1 and 2, and Exponential) to search the SRIs and response variables with an adjusted determination coefficient higher than 0.25 (**Figure 7A**). When Export data is selected, all relationships with adjusted determination coefficients higher than 0.25 will be reported (**Figure 7B**); the results, which are organized according to the loaded variables (column A) and SRIs (column B), show which regression model had the highest determination coefficient (Best) for each SRI, as well as its statistics [adjusted and non-adjusted determination coefficient, root mean squared error (RMSE), sum of squares due to errors (SSE) and degree of freedom (DFE)] (columns C–H). The results for each evaluated regression model are also described (Polynomial 1: columns I–M; Polynomial 2: columns N–R, and so on). In this screen example, the adjusted R 2 varied between 0.257 (Datt 850;710;680) and 0.406 (DLAI 1725;970), with these SRIs having the highest and lowest RMSEs, respectively (**Figure 7B**).

The selection of Open selection dialog (Detailed index report) enables the user to select specific SRIs and response variables for figure elaboration (Index report, **Figure 7C**). In this operational example, three SRIs (AI, BI, and CI) and one response variable (Yield) were selected. The SRI value for each scan or plot is given (**Figure 7D**) so the user can generate XY scatter plots for each tested SRI (X) and response variable (Y).

### CONCLUSIONS

Spectral Knowledge (SK-UTALCA) is a software package that allows an easy and fast exploratory analysis of highresolution spectral reflectance data, providing the user with tools to detect measurement problems and the generation of key information for later modeling. SK-UTALCA is especially useful for plant breeding or any other research area where the number of measurements (big data files) involves long working hours that increase the risk of making involuntarily mistakes. This freely-available software is the result of several years of measurements and analysis of spectral data oriented toward the prediction of traits in plant breeding.

### AUTHOR CONTRIBUTIONS

GL, CP-E contributed equally toward the intellectual input into the final version of this paper, including development of the software, data analysis, interpreting, and discussion of the results, and writing and editing the manuscript.

### ACKNOWLEDGMENTS

This work was supported and financed by the research program "Adaptation of Agriculture to Climate Change (A2C2)," "Nucleo Científico Multidisciplinario" and the "Vicerrectoría de Innovación, Desarrollo y Transferencia Tecnológica (VIDTT)" from Universidad de Talca. We also received significant funds from the National Commission for Scientific and Technological Research CONICYT (FONDEF IDEA 14I10106 and FONDECYT N◦ 11130601). We would like to express our gratitude to Sebastian Romero, Félix Estrada, Miguel Garriga, and especially to Alejandro Escobar for continued technical assistance in field experiments and laboratory analysis. We thank also Rodrigo Aguilar and Felipe Ojeda for their outstanding programming work, Genberries Ltda. for equipment support, and Ivan Matus (INIA- Chile) for genetic material and trial maintenance. Finally, our special gratitude goes to Greg Matyjewicz (ASD Inc., Boulder, CO, USA) for technical definitions and valuable discussion.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 01996/full#supplementary-material

### Software Availability

The compiled version of the software will be available for free downloading at http://www.fenomica.utalca.cl/ and source code is available as Supplementary Material.

### REFERENCES

the chlorophyll content in plants. Remote Sens. Environ. 69, 296–302. doi: 10.1016/S0034-4257(99)00023-1

canopy scale derived from hyperspectral and CO2 flux measurements in rice. Remote Sens. Environ. 112, 156–172. doi: 10.1016/j.rse.2007.04.011

spectrometry and radiometry. Turk. J. Agric. For. 30, 421–428. Retrieved from: http://dergipark.gov.tr/tbtkagriculture/issue/11618/138426

**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 © 2017 Lobos and Poblete-Echeverría. 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) or licensor 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.

# A Comprehensive Phenotypic Investigation of the "Pod-Shattering Syndrome" in Common Bean

Maria L. Murgia1 , Giovanna Attene1 , Monica Rodriguez 1 , Elena Bitocchi 2 , Elisa Bellucci 2 , Davide Fois 1 , Laura Nanni 2 , Tania Gioia3 , Diego M. Albani 4 , Roberto Papa2 \* and Domenico Rau1 \*

1 Dipartimento di Agraria, Sezione di Agronomia, Colture Erbacee e Genetica, Università degli Studi di Sassari, Sassari, Italy, 2 Dipartimento di Scienze Agrarie, Alimentari ed Ambientali, Università Politecnica delle Marche, Ancona, Italy, 3 Scuola di Scienze Agrarie, Alimentari, Forestali ed Ambientali, Università degli Studi della Basilicata, Potenza, Italy, 4 Dipartimento di Agraria, Sezione di Economia e Sistemi Arborei e Forestali, Università degli Studi di Sassari, Sassari, Italy

#### Edited by:

John Doonan, Aberystwyth University, UK

#### Reviewed by:

M. Teresa Sanchez-Ballesta, Institute of Food Science, Technology and Nutrition (CSIC), Spain Yuhui Chen, Samuel Roberts Noble Foundation, USA

\*Correspondence:

Roberto Papa<EMAIL_ADDRESS>Domenico Rau<EMAIL_ADDRESS>
#### Specialty section:

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

> Received: 30 September 2016 Accepted: 09 February 2017 Published: 03 March 2017

#### Citation:

Murgia ML, Attene G, Rodriguez M, Bitocchi E, Bellucci E, Fois D, Nanni L, Gioia T, Albani DM, Papa R and Rau D (2017) A Comprehensive Phenotypic Investigation of the "Pod-Shattering Syndrome" in Common Bean. Front. Plant Sci. 8:251. doi: 10.3389/fpls.2017.00251 Seed shattering in crops is a key domestication trait due to its relevance for seed dispersal, yield, and fundamental questions in evolution (e.g., convergent evolution). Here, we focused on pod shattering in common bean (Phaseolus vulgaris L.), the most important legume crop for human consuption in the world. With this main aim, we developed a methodological pipeline that comprises a thorough characterization under field conditions, including also the chemical composition and histological analysis of the pod valves. The pipeline was developed based on the assumption that the shattering trait itself can be treated in principle as a "syndrome" (i.e., a set of correlated different traits) at the pod level. We characterized a population of 267 introgression lines that were developed ad-hoc to study shattering in common bean. Three main objectives were sought: (1) to dissect the shattering trait into its "components," of level (percentage of shattering pods per plant) and mode (percentage of pods with twisting or non-twisting valves); (2) to test whether shattering is associated to the chemical composition and/or the histological characteristics of the pod valves; and (3) to test the associations between shattering and other plant traits. We can conclude the following: Very high shattering levels can be achieved in different modes; shattering resistance is mainly a qualitative trait; and high shattering levels is correlated with high carbon and lignin contents of the pod valves and with specific histological charaterstics of the ventral sheath and the inner fibrous layer of the pod wall. Our data also suggest that shattering comes with a "cost," as it is associated with low pod size, low seed weight per pod, high pod weight, and low seed to pod-valves ratio; indeed, it can be more exaustively described as a syndrome at the pod level. Our work suggests that the valve chemical composition (i.e., carbon and lignin content) can be used for a high troughput phenotyping procedures for shattering phenotyping. Finally, we believe that the application of our pipeline will greatly facilitate comparative studies among legume crops, and gene tagging.

Keywords: domestication, domestication syndrome, shattering, common bean, phenotypic analysis, element composition analysis, cell wall analysis

## INTRODUCTION

The loss of seed shattering occurred independently in several crops and in different areas of the world during the domestication of many food crops, as this loss was crucial for adaptation of the plants to the agro-ecosystem, to provide ancient farmers with easier and more abundant harvests (Tang et al., 2013). Non-shattering/indehiscent types emerged in maize, barley, and rice (see Li and Olsen, 2016, for a review). Maize was domesticated in the New World, in Mexico, while barley and rice were domesticated in the Fertile Crescent of the Old World and in south-east Asia, respectively. Similarly, among the legume crops, indehiscent phenotypes emerged in soybean and common bean, which were domesticated in the Old World and the New World, respectively (Hymowitz, 1970; Harlan, 1992; Bitocchi et al., 2012; Schmutz et al., 2014). However, fully indehiscent phenotype emerged in common bean only after domestication with the development of snap varieties that are used for the production of green beans due to the absence of fiber strings along the pod valves. In other domesticated commercial classes (e.g., dry beans) shattering traits it is only reduced from that observable in wild populations.

Thus, deciphering the genetic basis of pod shattering is important for evolutionary studies, particularly to unravel the mechanisms of parallel evolution (Lin et al., 2012; Dong and Wang, 2015), and also because this will provide breeders with key information to manipulate this trait to reduce yield loss (Singh, 2001; Santalla et al., 2004)**.** For the same reason, the genomic information would be a great tool to facilitate the exploitation of exotic germplasm in common bean breeding. The potential of these studies is well-represented by those that have been conducted in cereals (Lin et al., 2012). However, in legumes, "studies of the identification of pod-shattering genes lag far behind those of the cereal crops" (Li and Olsen, 2016).

The shattering system of legume crops is distinct from that of cereals (Li and Olsen, 2016). In legumes, dehiscence is subsequent to the "hygroscopic movement" of the pod valves following dehydration. The release of the accumulated elastic tension during dehydration results in the splitting of the valves along their suture lines (Elbaum and Abraham, 2014). The ability to undergo this movement has often been attributed to specific patterns of lignification of the pod-valve tissues.

Among legumes, The most relevant studies on pod dehiscence have been conducted in soybean. Histological analysis has shown that shattering wild genotypes differ from non-shattering varieties in terms of the degree of lignification of the cells along the suture lines of the pod valves (Dong et al., 2014). Among the cultivated germplasm, differential lignification of the lignin-rich inner sclerenchyma of the pod walls also influences the level of shattering (Funatsuki et al., 2014). Single major genes underlying these histological differences have also been cloned (Dong et al., 2014; Funatsuki et al., 2014). The loss of pod dehiscence has been studied to some extent in lupin, chickpea, pigeonpea, pea yardlong bean, and wild cowpea (Ladizinsky, 1979; Muehlbauer et al., 1998; Boersma et al., 2007, 2009; Weeden, 2007; Abbo et al., 2009; Suanum et al., 2016).

In common bean, there are few such data available. The pioneering studies date to almost a century ago (Lamprecht, 1932; Prakken, 1934). These attempted to index pod-shattering resistance not only based on the occurrence of valve splitting (presence/absence), but also depending on the mode of shattering; i.e., based on the degree of torsion (twisting/spiral coiling) of the pod valves after dehiscence (Lamprecht, 1932), and on suggested histological differences between shattering and non-shattering types, mainly in the lignification patterns of the valves tissues (Prakken, 1934). Oligogenic (Lamprecht, 1932) and monogenic (Prakken, 1934) bases for the genetic control of this trait were also proposed. Several decades later, in the pioneer study of Koinange et al. (1996) the pod strings locus (St) was mapped on chromosome 2, and it was proposed to control the differences in shattering between the wax snapbean Midas, an Andean commercial cultivar, and the wild Mesoamerican accession G12873 (Koinange et al., 1996). This locus did not cosegregate with two candidate genes PvSHP1 and PvIND, even if PvIND is linked to the St locus (Nanni et al., 2011; Gioia et al., 2012).

The aim of this study was to conduct a comprehensive phenotypic investigation of pod shattering in common bean. With this aim, we also set up a phenotyping pipeline that comprises characterization under field conditions, including the chemical composition, and histological analysis of the pod valves. Following this pipeline, we characterized a population of 267 introgression lines (ILs) that were developed ad-hoc to study pod shattering in common bean. In more detail, we pursued the following three goals:

### MATERIALS AND METHODS

### Plant Materials

A population of 267 introgression lines was phenotyped, which was representative of a larger set of about 1200 introgression lines developed in the Papa laboratory (Università Politecnica delle Marche, Ancona, Italy) in collaboration with the Attene laboratory (Università degli Studi di Sassari, Sassari, Italy). The population was developed starting from a backcross between the line MG38 belonging to the recombinant inbred line (RIL) population used by Koinange et al. (1996) and the recurrent parent MIDAS (**Figure 1**). The MG38 line is a RIL obtained from a cross between the wild Mesoamerican genotype, G12873, and the Andean snap bean variety MIDAS. The MG38 genotype was selected for some wild pod traits (small size, curved shape, pigmented valves, pod shattering), and seed characteristics (very small size). However, for other traits (e.g., determinacy, seed dormancy, photoperiod sensitivity), MG38 was selected for domesticated phenotypes to facilitate the population development and increase. Based on amplified fragment-length polymorphism analysis, MG38 has 55% of the genome attributable to the wild Mesoamerica parent G12873 (Papa, unpublished data). MIDAS is characterized by large and relatively straight, yellow and snap bean non-shattering pods, with relatively large seeds.

To obtain the introgression lines, MG38 was back-crossed with MIDAS as a recurrent parent, and different cycles of backcrossing and selfing were carried out together with selection for the wild characteristics of the pods and seeds. Among the 267 lines analyzed in this study, 70 belonged to BC3/F4:F5 families, and 217 to BC3/F6:F7 families. Overall, in the field 130 families were represented. Among these families, 101 families were represented by at least two ILs. In some case, there were three ILs per family (i.e., 29 families were represented by one individual). Precisely, there were 19 BC3/F4:F5 families and 82 BC3/F6:F7 families represented with at least two ILs summing up to 232 ILs.

### Phenotypic Characterization

The phenotyping data presented here were obtained in 2014, between May and October (sowing date, May 19). The experimental layout comprised eight rows, each with 35–38 holes; the distance between rows was 1.5 m; the distance between holes (within the rows) was 0.8 m. For each line, a single plant was grown in each single hole. The two parents, MIDAS and MG38, were replicated three times. The positions of the lines were completely randomized. A plastic sheet was used along each row to facilitate weed control (Supplementary Figure 1A). Standard agronomic practices were adopted, in terms of irrigation, fertilization, and pest control. The meteorological conditions were hot and dry with many days with maximum temperature over 30◦ (Supplementary Figure 1B). Under these conditions, ILs had the opportunity to fully express their shattering phenotype.

### Measuring Pod Shattering in the Field

We evaluated shattering after each plant reached full maturity. For each plant, we first distinguished between fertile and sterile pods. The numbers of "naturally" shattering and non-shattering pods were then counted. Fertile pods were further classified into different types, as exemplified in **Figure 2A**. Four pod categories were recognized: Indehiscent; "fissured," with valves that were not perfectly closed along the ventral suture; dehiscent

with non-twisting valves; and dehiscent with twisting valves. It was sometimes difficult to distinguish between these last two categories because of the presence of intermediate cases. Nonetheless, on the basis of this classification, the number of pods falling into each category was counted for each plant (**Figure 2B**). For the statistical analysis, the number of pods was expressed as the percentage of the fertile pods produced by each plant.

Furthermore, for each line separately, the shattering of indehiscent pods was promoted by hand, for the evaluation of the "resistance to manual shattering," based on a scale from 1 (i.e., very low resistance to shattering, where valves abruptly shattered under very light pressure on the distal part of the pod) to 9 (i.e., very high resistance to shattering, where valves did not separate and it was necessary to "break" them) (**Figure 2C**; see also Supplementary Information and videos). To avoid bias, the determination of the resistance to manual shattering was

"manual shattering" of the pods, into discrete scores from 1 (very easy) to 9 (very difficult) (see also Supplementary videos).

performed independently (i.e., at a different time) from the pod classification.

### Chemical Characterization of the Pod Valves

The chemical composition of the pod valves was investigated to determine whether the pod shattering was correlated to these characteristics. This element composition analysis looked at carbon, hydrogen, and nitrogen. Here, for each introgression line, 2 g dried pod valves was pulverized in a grinder (18,000 rpm, 1 min). The pulverized tissue was transferred into plastic 50-mL tubes and stored for a few days at room temperature in a cool, dry place. The analyses were performed using 0.080 g pulverized tissue from each line. The samples were combusted at 1,000◦C in an excess of oxygen using an element analyser (LECO CHN 628; Leco Corporation, St. Joseph, MI, USA), to determine the carbon, hydrogen, and nitrogen contents. The instrument was calibrated using the "oat meal 502276" forage standard with 46.43% carbon and 2.64% nitrogen. For each run, three independent samples of the standard were included.

The analysis was first performed considering the two parental lines, MIDAS and MG38, with each as three biological replicates (i.e., three plants were grown for each parent). For each biological replicate, there were three technical replicates (i.e., three independent analyses). As there were highly significant differences between the two parents (see Results), the analysis was extended to all of the introgression lines. Three technical replicates were performed for each introgression line.

### Cell-wall Analysis

For each individual plant, 6 g dried valves were pulversed in a mill (Retsch SM 100) for 10 min. The procedure of Van Soest and Wine (1967) was then followed. First, the neutral detergent fiber was quantified, which represents the total content in the cell wall of the analyzed sample. Thus, the acid detergent fiber was determined, which mainly represents an intermediate step that is necessary to extract the acid detergent lignin, which correlates with the lignin content of the sample analyzed. We also calculated the differences for the neutral detergent fiber minus the acid detergent fiber, and the acid detergent fiber minus the acid detergent lignin, which provided rough estimations of the hemicellulose and cellulose contents, respectively (Van Soest and Wine, 1967). All of these chemical fractions are expressed as percentages of the dried organic matter after subtracting the weight of the ashes (see Supplementary Information for further details).

This analysis was initially performed for MIDAS and MG38, for which three biological replicates were available. For each biological replicate, three technical replicates were included. As the analysis of variance (ANOVA) showed clear-cut differences between the parent lines, the cell-wall analysis was extended to 12 indehiscent introgression lines, and 12 high-shattering introgression lines (>65% shattering, as seen for MG38, the wild-like parent).

### Anatomical and Histological Study of the Pod Valves

This study had the specific aim to look for differential patterns between the shattering and non-shattering lines, particularly for lignin deposition. The analysis was conducted considering 5- and 20-days-old pods, and pods at the maturation stage. The pods were kept in a solution of 95% ethanol and glacial acetic acid (5:2, v/v) for 3 days, and then stored at 4◦C in 70% ethanol. Sections of the ventral and dorsal suture sheath were obtained manually. The sections were treated with Javelle water (an aqueous solution containing sodium hypochlorite and some sodium chloride, used as a bleach and disinfectant) for ∼10 min. After this washing, the sections were immersed in 50% acetic acid for a few minutes.

The pod valves were also embedded in paraffin, and 10-µm sections were obtained using a sliding microtome (Reicher) (see Supplementary Information for further details). The manually obtained sections were stained according to two different methods: Toluidine blue O (TBO), and carmine-iodine green; whole microtome sections were stained only with toluidine blue O. The toluidine blue O was used to differentially stain polysaccharides and lignin, whereby cells with thick lignified walls are sky blue, and cellulose and hemicellulose are dark blue (Mitra and Loqué, 2014). With carmine-iodine, lignin is green, and cellulose is pink (Deysson, 1954).

### Phenotyping of the Other Plant Characteristics

To allow the study of the relationships between shattering and the other plant traits, a total of 27 traits were recorded (7 qualitative, 20 quantitative). These were: Number of cotyledonary leaves (two, three); angle of the cotyledonary leaves (60◦ , 120◦ , 180◦ ); lobature of the cotyledonary leaves; stem color (green, red); growth type (non-climbing, intermediate, climbing); flower color (white, light purple, dark purple); pod color (yellow, striped, with from 1 to 3 stripes); plant height (cm); plant vigor (height per width; cm2 ); flowering time and pod setting (days from May 19); pod weight per plant (g); valve weight per plant (g); seed weight per plant (g); number of pods per plant; number of seeds per plant; mean pod weight (g); mean valve weight (g); 100-seed weight (g); weight of seeds per pod (g); number of seeds per pod; Harvest Index at pod level. To avoid loss of seeds at the maturation stage in the shattering plants, mature pods were (almost) enveloped in plastic nets (Supplementary Figure 2). Moreover, at the end of the ripening stage (i.e., before shattering occurred), 10 pods per introgression line were randomly sampled. These were scanned, and the acquired images (600 dpi) were processed with the Tomato Analyzer software (Rodríguez et al., 2010), to determine the following pod traits: Perimeters; area; curved height; maximum height; maximum width (Supplementary Figure 3). All of these measures were in pixels. We also calculated the ratio of the curved length to maximum height, where a ratio of 1.0 indicates a perfectly straight pod, while ratios <1 indicate a more or less marked "C" shapes of the pods. All of these variables must be referred to the projection area of the pod on the scanner glass. The procedure was first set up for the parental lines, MG38 and MIDAS. The analysis was then extended to all of the other introgression lines. For statistical analysis, 10 pods per introgression line were considered, and the means were calculated.

### Statistical Analysis

For each variable used to describe shattering, the frequency distribution was first determined. Associations between variables were quantified using Pearson "r" coefficient (quantitative traits) or contingency analysis (qualitative traits). Differences among groups of lines for the various phenotypic and chemical traits were tested using one-way analysis of variance (ANOVA), considering each line as a "replicate" of the group.

Resistance to manual shattering was modeled based on the other six indicators of pod shattering: Indehiscent (%); valves separated to some degree (%); fissured (%); shattering (%); nontwisting (%); twisting (%) (see **Figure 2**). With this aim, the method of recursive partitioning was adopted, which is also known as decision-trees analysis. This is particularly indicated to investigate relationships among variables without having an apriori model, and it is particular powerful as it considers a very high number of possible partitions, and takes into consideration only the best one (see JMP version 7, User Manual; SAS Institute Inc., Cary, NC, USA). In this case, the categorical X variable was the degree of "resistance to manual shattering" (scored as 1 to 9), while all of the other six indicators of shattering were considered as possible explanatory Y variables. Thus, it was possible to obtain a hierarchal system of (dichotomic) criteria that allowed the prediction of the manual shattering resistance from the observed level and mode of shattering. Statistical analysis of the phenotypic data was all performed using JMP version 7 (SAS Institute Inc., Cary, NC, USA).

## RESULTS

### Shattering Level and Mode

As expected, the two parental lines showed highly contrasting phenotypes for pod shattering: MIDAS was completely indehiscent, while MG38 was highly dehiscent, with a mean of 65% shattering pods per plant. Moreover, 98% of the variance for shattering occurrence was located among-families indicating a very limited role of environmental factors influencing this trait in the population of ILs grown under our field conditions (see Supplementary Information and Supplementary Table 1).

**Table 1** gives the descriptive statistics for the six variables measured to characterize the introgression lines for the podshattering trait.

The introgression lines were highly variable for the podshattering trait, as the indehiscent pods per plant ranged from 3.9 to 100%, with a mean of 50.4%. Twenty-nine introgression lines (∼10% of the total) were completely indehiscent. The pods per plant with valves separated to some degree ranged from 0 to 96.1%, with a mean of 49.6%. The distribution of these two variables tended to bimodality (Supplementary Figure 4). The fissured pods per plant ranged from 0 to 71.7%, with a mean of 18.0%.

The levels of shattering were highly variable, as the shattering pods per plant ranged from 0 to 82.6%, with a mean of 31.6%. The modes of shattering were also highly variable, as the non-twisting and twisting pods per plant both ranged from 0 to ∼60%. Nontwisting pods were more frequent than twisting pods, with means of 11.1 and 20.1%, respectively.

The distribution of the trait "resistance to manual shattering" is illustrated in **Figure 3**. MIDAS had a score of 8 (i.e., high resistance), while MG38 had a score of 2 (i.e., low resistance). The mean for this trait was 4.12 (i.e., medium-low resistance; σ = 1.96; S.E. = 0.12), and the distribution appeared to be bimodal. About 15% of the introgression lines showed scores of 1 and 2 (i.e., ≤MG38), while about 10% showed scores of 8 and 9 (i.e., ≥MIDAS).

### Relationships among the Measures of Shattering

**Figure 4A** shows the relationships between the levels and modes of shattering. Introgression lines with the same or very similar levels of shattering (percentage shattering pods per plant) showed a very different ratio between the twisting and non-twisting types. For example, among the introgression lines with very high levels of shattering (>65%), the ratio of non-twisting to twisting pods per plant varied from ∼25:50 (1:2) to about 60:15 (4:1). Moreover, these data also suggested that transgressive variation probably occurred for ∼10% of the lines, which showed higher shattering than MG38 (>65%), the highly shattering parental line (**Figure 4A**).

The level of shattering was more strongly correlated with the frequency of non-twisting pods (R 2 = 0.71, P < 10−4 ) than with the frequency of twisting pods (R 2 = 0.57; P < 10−4 ) (**Figure 4B**). In particular, while a low number of twisting pods corresponded to different levels of shattering, a low number of non-twisting pods was more indicative of low levels of shattering.

The resistance to manual shattering was modeled considering the six variables measured to dissect out the shattering trait (**Figure 2**, **Table 1**), with recursive partition analysis applied (**Figure 5**). The variable that best predicted resistance to manual shattering was the shattering pods per plant; i.e., the level of shattering. Indeed, the threshold of 10% shattering pods defined two groups of plants with different mean shattering resistance scores of 3.3 and 7.1; this partition captured 65% of the total variance for shattering resistance (P < 0.0001; **Figure 5**). A second partition suggested a role for the mode of shattering. Indeed, within the group of introgression lines showing <10% shattering pods per plant, a threshold of 9% twisting pods defined two subgroups of plants that had mean shattering resistance scores of 2.9 and 4.1. This partition captured a small portion, 6%, of the total variance for resistance to manual shattering (P < 0.0001; **Figure 5**). A third partition was found within the group of introgression lines with shattering pods ≥10%. Here, the threshold of 4.2% of non-sigaroid pods defined two subgroups of plants with mean shattering resistance scores of 5.3 and 7.5, and these explained an additional 4% of the total variance for shattering resistance (**Figure 5**). Thus, cumulatively these three partitions explained 75% of the total variance. The

fourth partition (not shown) explained 0.8% of the total variance for shattering resistance, which indicated that dealing with a more complex model was not necessary.

### Chemical Analysis

The pod valves of the two parental lines, MG38 and MIDAS, had significantly different carbon contents (ANOVA, P < 0.0001; **Table 2**). The highly dehiscent MG38 had a carbon content of 43.8% dry weight, which gave a 6.8% increase in the carbon content of the indehiscent MIDAS, from 41.0% dry weight (**Table 3**). ANOVA also revealed a marginally significant difference for the hydrogen contents (P < 0.047; **Table 2**), again in favor of MG38 (6.7% dry weight) compared to MIDAS (6.5% dry weight; **Table 3**). The difference in the nitrogen contents was not significant (P = 0.502) (**Tables 2, 3**).

The comparison of the indehiscent vs. dehiscent introgression lines was highly significant for the carbon contents (**Table 2**). The dehiscent introgression lines showed a 6.9% increase in the carbon content of the indehiscent introgression lines, according to dry weight (**Table 3**). The indehiscent introgression lines had the same carbon content as MIDAS, while the dehiscent introgression lines had the same carbon content as MG38. The difference between these dehiscent and indehiscent introgression lines was small, but significant for the hydrogen content although not for the nitrogen content (**Tables 2, 3**). The frequency distribution for the carbon contents tended to be bimodal, while this was less evident for the hydrogen and nitrogen contents (Supplementary Figure 5).

The relationship between the carbon contents and the shattering pods per plant (**Figure 6**) showed an abrupt transition in the carbon content that occurred between 5 and 10% shattering pods per plant (**Figure 6A**). Partition analysis showed that this transition occurred at 7.14% shattering pods per plant (not shown). The definition of the introgression lines into two groups based on this transition, with the first with <7.14% and the second with ≥7.14% shattering pods per plant, captured 47% of the total variance for the carbon contents (**Table 2**). When the introgression lines with <7.14% shattering pods were excluded from the analysis, there was a weak, but significant, negative correlation (r = −0.296; P < 0.0001) between shattering level and carbon content (**Figure 6B**).

non-twisting). Deep blue, no shattering; red, high shattering. (B) Left: Association between levels of shattering and frequencies of twisting pods. Right: Association between levels of shattering and frequencies of non-twisting pods. In both cases, linear and smoothing spline (λ = 10,000) fits are shown. The associations were tested excluding the completely indehiscent plants.

### Cell-wall Analysis

The pod valves of MIDAS and MG38 had significantly different total fiber contents (ANOVA, P < 0.0001; **Table 4**). The fiber content of the highly dehiscent MG38 (62.0%) showed a ∼48% increase of that of the completely indehiscent MIDAS (42.0%; **Figure 7A**). The contents of the three cellwall components of lignin, hemicellulose, and cellulose were always higher for MG38 than MIDAS, with the greatest difference seen for lignin, followed by hemicellulose and cellulose (**Figures 7B–D**). Statistical analysis was carried out to compare the two groups of introgression lines, the first comprising the completely indehiscent plants (i.e., non-shattering), and the second including the plants with >65% shattering pods per plant (i.e., high shattering) (**Table 4**, **Figures 7A–D**). These two groups strongly differed in their lignin contents, with the difference for the latter representing a 180% increase of the former, which was much greater than the increase from the MIDAS to MG38 parental lines (80%) (**Figure 7B**). The high shattering group also showed increases of the non-shattering group for hemicellulose (33.9%) and cellulose (7.6%) contents, which were here less than for the parental lines (79.5%, 30.1%, respectively) (**Figures 7C,D**). These highly dehiscent (i.e., high shattering) introgression lines showed lower total fiber content than MG38 (**Figure 7A**), which was mainly due to reduction in the cellulose content, and to a slight, although not statistically significant, reduction in the hemicellulose content (**Figure 7D**). In contrast, these lines had higher lignin content compared to MG38 (**Figure 7B**). This suggested that the achievement of the very high pod shattering ability (here even higher than the wildlike MG38) is associated with an increase in the proportion of lignin in the cell wall.

### Correlations between the Element Compositions and the Cell-wall Analysis

The carbon content was strongly correlated with the total fiber content of the pod valves (r = 0.685); moreover, among the three cell-wall components (i.e., lignin, hemicellulose, cellulose), the carbon content showed the best correlation with the lignin content (r = 0.672; **Table 5**). Stepwise multiple regression analysis was performed with the carbon content as the dependent variable and lignin, hemicellulose, and cellulose as the independent variables (Supplementary Table 2). Here, the only variable that entered into the model was the lignin content (P = 1.85 × 10−5 ). This thus indicates that the hemicellulose and cellulose correlation to the carbon content was mainly due to their correlation with the lignin content.

### Anatomical and Histological Analysis of the Pod Valves

The analysis conducted with 5-day-old pods showed no obvious differences between the shattering and non-shattering genotypes (Supplementary Figure 6). The ventral sheath showed only a few cells with very low levels of lignification. A similar situation was observed for the dorsal sheath (not shown). There was no lignin deposition in correspondence with the inner parenchyma cells of the pod walls.

Analyses of 20-day-old pods showed evident lignin deposition in the ventral sheath of the pod valves, and a clear-cut difference between the shattering (i.e., MG38) and nonshattering (i.e., MIDAS) genotypes (**Figures 8A,B**). Indeed, the proportion of cells with thick secondary cell-wall formation (i.e., sclerenchymatic cells) was clearly greater for MG38 (highly dehiscent), compared to MIDAS (indehiscent). The absence of cells with thick secondary cell-wall formation for MG38 was limited only to the external layer of the cells of the sheath and to the dehiscence zone, while for MIDAS this involved all of the sheath (**Figure 8B**). Moreover, for MG38, the cell-wall thickness tended to reduce when moving from the sheath to the dehiscence zone (**Figure 8A**), where there was the tendency to easily "fracture" (**Figure 8A**). A similar pattern was observed for the dorsal sheath (Supplementary Figure 7). A clear-cut difference was also seen between the parental MG38 and MIDAS for the degree of lignification in the inner cells of the pod walls, with very strong lignification (i.e., sclerenchyma) for MG38, and complete absence of lignin deposition for MIDAS (**Figure 9**).

As not all anatomical or histological differences between the wild-like parent (MG38) and the cultivated varieties

TABLE 2 | Results for the ANOVA performed for the chemical element analysis of the pod valves.

ANOVA was applied to the following comparisons: Parental lines MG38 (dehiscent) vs. MIDAS (indehiscent); indehiscent vs. dehiscent introgression lines; and introgression lines with <7.14 vs. ≥7.14% shattering pods. M.S., Mean Square; Fx/y, F ratio with x and y degrees of freedom for the numerator and denminator, respectively. R2 adj, adjusted R2 ; P, significance level.

(MIDAS) are necessarily correlated with the shattering traits, two introgression lines (one with shattering >MG38, the other without shattering) were also compared (Supplementary Figures 8A–C). Encouragingly, here the patterns were similar to those observed between MG38 and MIDAS, which suggests that the histological differences seen do indeed underlie the shattering/non-shattering phenotypes.

It was difficult to obtain good sections of the pod valves at the maturation stage because of the fragility of the tissue. However, it can be noted that at this stage, the ventral sheath of MIDAS had

TABLE 3 | Mean contents of carbon, hydrogen, and nitrogen of the pod valves as applied to the following comparisons: two parental lines MG38 (dehiscent) vs. MIDAS (indehiscent); indehiscent vs. dehiscent introgression lines; and introgression lines with <7.14 vs. ≥7.14% shattering pods.

Groups of introgression lines with different letters have significantly different means (P <0.05; Tukey-Kramer multiple comparison tests).

FIGURE 6 | (A) Relationship between carbon contents and levels of shattering. Green shading, individuals for which shattering was <7.14%; orange shading, individuals for which shattering was ≥7.14%. (B) Relationship between carbon contents and shattering levels excluding individuals with low or no shattering. R 2 given for smoothing spline (λ = 10,000) (red) and linear fit (green).

ANOVA was applied to the following comparisons: Parental lines MG38 (dehiscent) vs. MIDAS (indehiscent); and two groups of 12 introgression lines, one completely indehiscent vs. one highly dehiscent (i.e., with percentage shattering pods >65% of MG38). M.S., Mean Square; Fx/y, F ratio with x and y degrees of freedom for the numerator and denominator, respectively. R2 adj, adjusted R2 ; P, significance level.

NDF, neutral detergent fiber; ADL, acid detergent lignin; ADF, acid detergent fiber.

more mechanical resistance than that of MG38, which appeared to be very fragile instead (Supplementary Figure 9).

### Relationships between Pod Shattering and the Other Plant Characteristics

**Table 6** gives the associations between the levels of shattering (percentage shattering pods per plant) and the other 28 phenotypic traits, of which seven are qualitative and 21 are quantitative. These included morphological, phenological, and productive traits. Overall, shattering was very poorly correlated with all of these plant traits considered. Significant weak associations were detected for three qualitative traits, four quantitative–productive traits, and six precision phenotyping traits that describe pod size and shape (**Table 6**).

For the associations with the qualitative traits, it was observed that: Plants with red stems showed higher shattering (59.6 ± 12.66%) than those with green stems (31.25 ± 1.30%); plants with white flowers showed higher shattering (36.52 ± 1.44%) than those with purple (22.85 ± 2.67%) and light purple (13.66 ± 3.96%) flowers; plants with yellow pods showed higher shattering (35.89 ± 1.44%) than those with striped/yellow pods (17.58 ± 2.80%), with an intermediate position seen for those with striped pods (30.55 ± 6.33%) and yellow/ stripped pods (24.92 ± 5.77%).

For the correlations with the productive traits, valve weight per plant and mean valve weight increased when the shattering level increased, with the opposite for 100-seed weight and Harvest Index at pod level, which decreased when the shattering level increased (**Table 6**). In more detail, the oneway ANOVA between the dehiscent vs. indehiscent introgression lines showed that the former had significantly higher valve weight per plant and mean valve weight than the latter, with increases in the indehiscent introgression lines of 35.4% (P = 0.0085; t-test) and 26.2% (P = 0.0004). In contrast, the opposite was seen for the 100-seed weight and the Harvest Index at pod level, where the indehiscent introgression lines showed an increase of 15.9% (P = 0.0002) and 8.9% (P < 0.0001), respectively, to the dehiscent introgression lines.

The correlations between the levels of shattering with the six precision phenotyping variables that describe the pod morphology were significant (from P < 0.0001 to = 0.021) and all negative (**Table 6**).

### DISCUSSION

### Field Phenotyping of the Shattering Trait in Common Bean

High variations for both levels and modes of pod shattering were recorded. All of the shattering types were distinguishable, which varied from completely indehiscent to "twisting," passing through the two defined "intermediate" states of "fissured" and "shattering but non-twisting" (Lamprecht, 1932). Each introgression line was characterized by counting and classifying the pods into these four categories, with the degree of resistance to manual shattering also independently measured.

As shown by the partition analysis, the best predictor of resistance to manual shattering was the level of shattering TABLE 5 | Correlations between the element compositions and cell-wall fiber contents for the pod valves of the 24 introgression lines, 12 non-shattering, and 12 very high shattering.

NDF, neutral detergent fiber; ADL, acid detergent lignin; ADF, acid detergent fiber; \*P < 0.05; \*\*P < 0.01; \*\*\*P < 0.001.

(percentage of shattering pods per plant), while the mode of shattering (twisting/non-twisting) was less relevant. Moreover, a low threshold of shattering pods per plant (10%) was sufficient to distinguish between the low and medium-high resistant introgression lines. All this suggests that shattering might be controlled by the "switching" of the mechanism of control that determines the abrupt change in the possibility of splitting the pod valves. These data also indicate that when considering both natural or artificial plant populations, genetic studies aimed at deciphering the genetic architecture of the pod shattering trait would benefit from a step-wise approach that comprises the following: (1) comparing indehiscent vs. dehiscent introgression lines (regardless of the degree of shattering); (2) considering only dehiscent introgression lines (regardless of the mode); and (3) considering separately among the dehiscent introgression lines those with twisting and non-twisting pods. Indeed, this approach would allow the genetic basis of the occurrence of shattering (yes/no) and also its tuning (low/high) and mode (twisting/nontwisting) to be described. It should also be noted that the variable of "fissured pods" did not prove useful to predict resistance to manual shattering; this suggests that this trait would be better investigated separately from the others.

### Element Composition and Cell-wall Analysis of the Pod Valves

These shattering and non-shattering genotypes clearly differed in their carbon contents. The contents of carbon, hydrogen, and nitrogen are expected to be stoichiometrically correlated to the amount of organic matter in the tissues (Chiariello et al., 2000), and thus to the cumulative content of carbohydrate, protein, lipid, and all other organic compounds. However, in plants, differences in the carbon content have frequently been correlated to differences in lignin content (Loader et al., 2003). This was also the case for the valves of common bean; indeed, the cell-wall analysis here confirmed that the differences in the carbon contents between the shattering and non-shattering types were mainly correlated with the differences in the lignin contents, in comparison with the other cell-wall components (i.e., hemicellulose, cellulose).

There was an abrupt increase in the carbon content at a level of shattering of ∼7.14%. This value was similar to the threshold of 10% of the shattering pods per plants that explained the largest proportion (65%) of the variance for

showing ventral sheath of pod valves from MG38 (highly dehiscent), and details of dehiscent zone, thick lignified fibers (sclerenchyma), wood cells, dehiscent zone after cracking (arrow). (B) TBO staining, showing ventral sheath of pod valves from MIDAS (indehiscent), with two details of the indehiscence zone. LF, lignified fibers; WC, wood cells.

resistance to manual shattering. These observations suggest that environmental effects might act on the level of shattering, and that the complementation of the whole-plant characterization and chemical element composition analysis can lead to more precise and alternative or complementary phenotyping option.

The comparison of the data in the present study with those from the literature reveal differences between common bean and

soybean. Indeed, in soybean, the high-shattering cultivars were shown to have similar lignin contents (not higher, as observed here for common bean) to the low-shattering cultivars (see **Table 1** of Romkaew et al., 2008). Moreover, in F2 and backcrossed populations between yardlong and wild cowpea, among these three fiber components, the contents of hemicellulose showed the highest correlation with pod shattering (Suanum et al., 2016). All this suggests that there are histological differences between other legumes species and common bean, which appear to be due to differences in the patterns of cell-wall lignification of the pod tissues, or differences in the prevalent fiber type.

### Histological Characterization of the Pod Valves

The data from the histological characterization of pod valves in the present study parallel the observations at the chemical level. Indeed, overall, cell-wall lignification is much more pronounced in the shattering type than the non-shattering type for common bean. Specifically, the ventral sheath of the wild-like genotype (MG38) was characterized by very strong sclerenchymatization of the cells, while the opposite was seen for the domesticated cultivar (MIDAS). This observation is consistent with **Prakken (1934)**, who indicated this anatomical difference at the basis of the presence/absence of pod strings, and on the basis of the shattering/non-shattering phenotypes. Moreover, the presence/absence of pod strings has been shown to be under the control of the St gene (Koinange et al., 1996).

The histological differences between the shattering and nonshattering genotypes for common bean in the present study appear to be more pronounced than those for soybean. Indeed, for soybean, the differences were limited to the dehiscent zone, where excessive lignification of the fiber cap cells was seen in the cultivated non-shattering genotypes, as compared to the wild shattering genotypes (Dong et al., 2014). Furthermore, a major gene, known as SHA1-5, was identified as being responsible for lignin deposition in the fiber cap cells of soybean (2014). The present study did not show any clear histological differences in the common bean dehiscence zone. Albeit it cannot be completely excluded that there were some undetected histological differences here in the dehiscence zone in common bean, these data suggest that the histological basis of pod shattering in bean and soybean are different, at least partially.

Furthermore, the shattering genotypes here had a fibrous and strongly lignified cell layer between the inner and outer parenchyma of the pod wall, while this was not seen for the nonshattering genotypes. This difference was also noted for common bean by Prakken (1934), in their comparison of the "stringy" and "stringless" types. Funatsuki et al. (2014) noted that in cultivated soybean, the differential lignification in the lignin-rich inner sclerenchyma of the pod walls influenced valve twisting and pod shattering. Furthermore, they showed that lignin deposition in this layer was under the control of a major gene, known as PDH1. We note here that in soybean, the difference between the shattering and non-shattering types appears to be in the degree of lignification of this inner sclerenchyma layer (Funatsuki et al.,

Bold, significant associations. M.S., Mean Square; F, F ratio; R 2 adj, adjusted R2 ; r, Pearson correlation coefficient; P, significance level.

2014), while the present study indicates that in common bean this difference is much stronger, with the presence/absence of the lignified layer seen. This suggests that the role of this lignified layer in the shattering might be more relevant (or at least different) in common bean compared to soybean. This might mark another difference between these two closely related crops, of common bean and soybean. Prakken (1934) suggested that in common bean, the control of the traits of "stringlessness" (which depends on the characteristics of the ventral sheaths) and "parchment" (which depends on the layer between the inner and outer parenchyma of the pod wall) was independent, and in both cases was under simple monogenic control. Another study suggested oligogenic control for the stringless trait, with the contribution of either environmental effects or epistatic interactions (Dong and Wang, 2015). Thus, in common bean, the artificial selection might have targeted multiple genes to minimize the seed loss during domestication. This evokes a scenario that arises from the joint consideration of the data obtained in soybean by the independent studies of Dong et al. (2014) and Funatsuki et al. (2014).

Based on the data in the present study, two further conclusive considerations can be made that might be useful to support the identification of shattering genes in common bean. First, it is likely that the genes underlying the shattering trait in common bean are involved in the regulation of the secondary cell-wall deposition or fiber-cell differentiation. This is wellsupported by the data from the chemical analysis and the anatomical–histological investigations here; indeed, fibers are mainly composed of sclerenchymatic cells, that have welldeveloped secondary cell walls. This possibility is also suggested by the data of Suanum et al. (2016), who reported co-localization of QTLs for pod fiber content and pod shattering in back-cross populations between yearlong bean and wild cowpea. Secondly, as the comparison with the literature indicates some chemical and histological differences between soybean and common bean, it might be useful to consider as candidate genes not only those involved in the shattering of soybean, but also those from other phylogenetically more distant species (Dong and Wang, 2015; Li and Olsen, 2016).

### Relationships between Shattering and the Other Plant Traits

The shattering levels were very poorly correlated with the other morpho-phenological traits and productive characteristics of the plants investigated here. However, an interesting consideration arises from the observation that shattering is significantly (albeit poorly) associated with low 100-seed weight, small pod size, and low Harvest Index at pod level. This suggests that pod shattering might have an "energy cost" for the plant (McGinley and Charnov, 1988; Chiariello et al., 2000); i.e., the synthesis of the biomolecules and the creation of the tissues needed for shattering might reduce the resources available for seed and pod development. In agreement with this energy cost hypothesis, the carbon contents of the pod valves were strongly and positively correlated with the levels of shattering. This all suggests that shattering can be better viewed as a syndrome at the pod level. However, as the same data can be explained by pleiotropic effects or linkage drag, more data will need to be collected also in other species to confirm this hypothesis.

## CONCLUSIONS

Pod shattering in common bean was investigated in the present study. With this objective, we set up and adopted a pipeline for phenotypic characterization of this trait. Four main results were achieved: (1) very high shattering levels can be obtained with a high percentage of either twisting or non-twisting pods, or with a balanced combination between these two; i.e., in common bean, the modes of shattering do not have any great impact on the levels of shattering; (2) shattering appears to be controlled by a "switching" mechanism that determines an abrupt change in the ability to split the pod valves; (3) high shattering levels is correlated with high carbon and lignin contents of the pod valves, and with specific histological charaterstics of the ventral sheath and the inner sclerenchymatic layer of the pod wall; and (4) shattering appears to have a "cost.", and it might be more exhaustively described as a "syndrome" at the pod level.

Overall, our pipeline will help with the deciphering of the genetic architecture of shattering in different crops, thus facilitating comparative studies in legumes.

### AUTHOR CONTRIBUTIONS

Designed the project: DR and RP. Managed the project: DR, GA, and RP. Wrote the article: MM and DR. Contributed to the drafting and the critical revision of the article: DR, MM, MR, GA, EBi, EBe, DA, DF, LN, TG, and RP. Contributed plant materials: EBi, EBe, LN, GA, and RP. Performed phenotypic characterization under field conditions MM, DF, TG, and DR. Histological analysis: DA, MM, and DR. Collect data from element composition and cell-wall analyses: MM. Analyzed and interpreted data: DR, MM, MR, RP, and GA. Edited the article: DR, MM, MR, EBi, EBe, GA, and RP. All of the author approve the final version of the manuscript.

### REFERENCES

### ACKNOWLEDGMENTS

Thanks for technical assistance at University of Sassari to: A. Ara, B. Scalas, and M. Pinna for phenotyping under field conditions; M. Deroma for CHN analysis; R. Rubattu for cellwall analysis; G. Becca for histological analysis. This research represents part of the PhD project carried out by MM at the doctoral school of Science and Biotechnology of Agricultural and Forestry Science and Food Production, curriculum Crop Productivity of the University of Sassari (Supervisor: DR; Tutor: GA). MM gratefully acknowledge Sardinia Regional Government that partially funded her PhD scholarship (P.O.R. Sardegna F.S.E. 2007–2013—Obiettivo competitività regionale e occupazione, Asse IV Capitale umano, Linea di Attività l.3.1. Operational Programme of the Autonomous Region of Sardinia, European Social Fund 2007-2013—Axis IV Human Resources, Objective l.3, Line of Activity l.3.1.).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017. 00251/full#supplementary-material

dissemination of soybean. Proc. Natl. Acad. Sci. U.S.A. 111, 17797–17802. doi: 10.1073/pnas.1417282111

**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 © 2017 Murgia, Attene, Rodriguez, Bitocchi, Bellucci, Fois, Nanni, Gioia, Albani, Papa and Rau. 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) or licensor 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.

# Genomic and Phenomic Screens for Flower Related RING Type Ubiquitin E3 Ligases in Arabidopsis

Mirko Pavicic1, 2, Katriina Mouhu1, 2, Feng Wang1, 2, Marcelina Bilicka1, 2, Erik Chovancek ˇ 1 and Kristiina Himanen1, 2 \*

1 Department of Agricultural Sciences, University of Helsinki, Helsinki, Finland, 2 Viikki Plant Science Centre, University of Helsinki, Helsinki, Finland

Flowering time control integrates endogenous as well as environmental signals to promote flower development. The pathways and molecular networks involved are complex and integrate many modes of signal transduction. In plants ubiquitin mediated protein degradation pathway has been proposed to be as important mode of signaling as phosphorylation and transcription. To systematically study the role of ubiquitin signaling in the molecular regulation of flowering we have taken a genomic approach to identify flower related Ubiquitin Proteasome System components. As a large and versatile gene family the RING type ubiquitin E3 ligases were chosen as targets of the genomic screen. The complete list of Arabidopsis RING E3 ligases were retrieved and verified in the Arabidopsis genome v11 and their differential expression was used for their categorization into flower organs or developmental stages. Known regulators of flowering time or floral organ development were identified in these categories through literature search and representative mutants for each category were purchased for functional characterization by growth and morphological phenotyping. To this end, a workflow was developed for high throughput phenotypic screening of growth, morphology and flowering of nearly a thousand Arabidopsis plants in one experimental round.

#### Edited by:

Marcos Egea-Cortines, Universidad Politécnica de Cartagena, Spain

#### Reviewed by:

Federico Valverde, Consejo Superior de Investigaciones Científicas (CSIC), Spain Daniel Hofius, Swedish University of Agricultural Sciences, Sweden

#### \*Correspondence:

Kristiina Himanen<EMAIL_ADDRESS>
#### Specialty section:

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

> Received: 30 September 2016 Accepted: 10 March 2017 Published: 28 March 2017

#### Citation:

Pavicic M, Mouhu K, Wang F, Bilicka M, Chovancek E and ˇ Himanen K (2017) Genomic and Phenomic Screens for Flower Related RING Type Ubiquitin E3 Ligases in Arabidopsis. Front. Plant Sci. 8:416. doi: 10.3389/fpls.2017.00416 Keywords: Arabidopsis, flower, RING E3 ligase, ubiquitin, high throughput, image based phenotyping, phenomics data analysis

### INTRODUCTION

Flowering time control is a complex network that integrates many modes of signal transduction promoting transition from vegetative stage to reproduction and ultimately leading to the development of flower organs. The endogenous changes that signal the beginning of flowering are referred as autonomous pathways (Amasino and Michaels, 2010). Multiple studies have established the major role that photoperiod has in flowering (Piñeiro and Jarillo, 2013). Flowering in Arabidopsis is strongly promoted in long day (LD) conditions but will also ultimately occur under short day (SD) conditions (Steffen et al., 2014). Under LDs, flower induction is dependent on the expression and protein levels of CONSTANS (CO; Suárez-López et al., 2001). Light controls the CO transcription via the circadian clock system, inducing a CO mRNA peak during the latter part of the day (Suárez-López et al., 2001). CO transcription is repressed by CYCLING DOF FACTORs (CDFs; Fornara et al., 2009). Under LDs, the CO mRNA afternoon peak coincides with a blue-light activated complex containing FLAVIN-BINDING, KELCH REPEAT, F-BOX 1 (FKF1) and GIGANTEA (GI), which lead CO transcription repressors CDFs to degradation (Sawa et al., 2007; Fornara et al., 2009; Song et al., 2012). Additionally, the FKF1-GI complex also stabilizes CO protein in the afternoon (Sawa et al., 2007; Song et al., 2012). CO protein degradation is promoted by at least two ubiquitin E3 ligases: HIGH EXPRESSION OF OSMOTICALLY RESPONSIVE GENES 1 (HOS1) and CONSTITUTIVE PHOTOMORPHOGENIC 1 (COP1; Jang et al., 2008; Lazaro et al., 2012). In the morning, red light promotes HOS1 interaction with CO via phytochrome B activation (Lazaro et al., 2012). COP1 mediates CO protein degradation in a complex with SUPPRESSOR OF PHYA-105 (SPA; Laubinger et al., 2006). In the afternoon, blue light inhibits COP1-SPA-mediated CO degradation by activating CRYPTOCHROME 2 (CRY2) interaction with COP1 (Liu et al., 2008). Thus, both CO transcription is up-regulated and CO protein stabilized allowing up-regulation of the mobile flowering signal gene FLOWERING LOCUS T (FT) expression in the phloem during the afternoon under LDs, but not under SDs (Piñeiro and Jarillo, 2013). Also regulation of flower development is likely to involve Ubiquitin Proteasome System (UPS) components (Vierstra, 2009).

The UPS has emerged as a powerful regulatory mechanism that facilitates irreversible transitions between developmental stages, and responses to environmental stimuli by selectively degrading short-lived regulators, such as transcription factors and receptors (Sadanandom et al., 2012). Genetic analyses in plants have proposed that this pathway plays a vital role in hormone regulation, floral homeostasis, stress responses and pathogen defense; however, very few targets have been identified in plants apart from the hormone signaling components (Santner and Estelle, 2010). In the UPS system, the highly conserved 76-amino acid protein, ubiquitin, acts as a covalent molecular tag to signal target proteins for proteasome mediated degradation. Ubiquitin attachment requires three distinct enzymatic activities: E1, ubiquitin activating enzymes; E2, ubiquitin conjugating enzymes; and E3, ubiquitin ligase enzymes. Moreover, the UPS consists of accompanying proteins that modulate target recognition and degradation (such as RAD23, SPA1), deubiquitinating enzymes (DUB1) and the proteasome (26S and 20S structures). According Vierstra (2009) over 6% of the Arabidopsis proteome is potentially involved in UPS. However, the common strategy for functionally addressing the role of all UPS components is still evolving. The ubiquitin E3 ligases are the most abundant UPS components and mediate the important recognition of the target proteins for ubiquitination (Kosarev et al., 2002; Stone et al., 2005). The E3 ligases found in plants belong to one of four subtypes: single subunit Homology to E6-AP C-Terminus (HECT), U-box and Really Interesting New Gene (RING) or multisubunit cullin-RING ligases (Sadanandom et al., 2012). The RING-type E3 proteins are the most abundant among the single subunit E3 ligases (Kosarev et al., 2002; Stone et al., 2005).

To unravel the role of the RING type ubiquitin E3 ligase protein family, we took a reverse-genetics approach to identify the RING E3 ligases that could be involved in regulation of Arabidopsis flowering time and/or flower development. To this end, we first curated the RING E3 protein family, earlier described by Stone et al. (2005), in the most recent Arabidopsis genome release. The Arabidopsis protein sequences were subjected to InterProScan for protein domain search and the number of ubiquitin E3 ligases containing RING domains was established to be 509. Association of these RING protein encoding genes with Arabidopsis flowering and floral organs was done through the Genevestigator transcriptome database (Hruz et al., 2008). To this end, the expression profiles were divided into categories based on their specificity, high expression or enrichment in flower organs and in the developmental stages of Arabidopsis. Several already characterized regulators were identified among these genes, such as the anther dehiscence regulating DAF gene family (Peng et al., 2013), flower size regulating DA2 (Xia et al., 2013) and FRG1 involved in flowering time related DNA methylation (Groth et al., 2014). The wellestablished flowering time regulator COP1 fell just below the cut off criteria due to its wide expression profile. A representative mutant collection for each category was obtained from NASC stock center. Additional candidates were also selected based on literature. The genotypically verified mutant collection was subjected to systematic morphological and growth analysis using an automated imaging based plant phenotyping facility. After the thorough vegetative assessment, the flowering time parameters such as number of leaves at bolting and days to bolting were recorded together with morphological analysis of the flower structures. The phenotypic assessment indicated lines with altered growth, morphology, or flowering time. Furthermore, one of the lines showed growth defects in sepals and petals.

### MATERIALS AND METHODS

### Bioinformatic Screens

To curate putative RING-type ubiquitin E3 ligases in Arabidopsis thaliana genome version ARA11, classification made by Kosarev et al. (2002) and Stone et al. (2005) were used. To this end, the whole Arabidopsis proteome was downloaded from ARAPORT (https://www.araport.org/downloads/), and screened with InterProScan for protein families and domain architecture. To confirm that the newly identified RING domain containing protein sequences indeed represented ubiquitin E3 ligase type RING domains, InterProScan 5 (v5.16-55.0) Gene3D, SUPERFAMILY, ProSiteProfiles, SMART, Pfam, and ProSitePatterns signatures were used. Most of InterProScan tools use Hidden Markov Models (HMMs) to detect conserved domains along protein sequences. HMMs have been developed for conserved protein domains and they define for the software, which and where critical residues should be located along the analyzed protein sequence. From the protein domain collection, the ubiquitin E3 ligase type RING domains were filtered according to the criteria provided by Kosarev et al. (2002) and Stone et al. (2005) for canonical RING domains. Once the RING domains were identified, they were aligned with Jalview using ProbCons algorithm with two rounds of pre-training. The metal ligand binding residues were manually inspected and corrected, and small misalignments were edited. Sequences that failed to meet the criteria of InterProScan search engines were not considered in this study.

### Transcriptomic Database Screens

To associate the curated collection of 509 RING type ubiquitin E3 ligases with flowering the Genevestigator gene expression database software was used (Hruz et al., 2008). The experiments AT-00087, AT-00088, AT-00089, and AT-00090 containing developmental expression data of AtGenExpress initiative microarrays were selected for the analysis (Schmid et al., 2005). In the selected experiments, hybridization probes were available for 393 RING E3 genes out of the 509. From these experiments, the linear expression data was extracted for the developmental stages of developed rosette, bolting, young flower, developed flower, and flower and silique. For flower organs, the gene expression profiles were extracted for categories of shoot apical meristem (SAM), sepal, petal, stamen, and carpels. In these categories, genes were ranked for their at least 2-fold differential expression against the developed rosette. Their relative expression levels were obtained by log2(FC) = log2(FL) − log2(R), where FC is fold of change, FL is flower organ or development stage and R is rosette. The results for each category were sorted by their log2(FC) and all genes with log2(FC) > 1 were considered as up-regulated.

### Candidate Genes Selected by Literature

For the candidate approach, we used interaction networks from BioGRID (http://thebiogrid.org/) and cross-checked them with flowering pathway genes listed in the Flowering Interactive Database FLOR-ID (Bouché et al., 2016) to identify RING E3 ligases interacting with known flowering time regulators CONSTANS (CO), CONSTITUTIVE PHOTOMORPHOGENIC 1 (COP1), and TARGET OF EARLY ACTIVATION TAGGED (EAT) 2 (TOE2).

### Plant Material and Growth Conditions

For functional characterization of the identified top most differentially expressed genes and for the selected candidates, Arabidopsis mutant lines were obtained from the NASC stock center representing CATMA, SAIL, SALK, and GABI-Kat collections (Alonso et al., 2003; Rosso et al., 2003; Schmid et al., 2005; Kleinboelting et al., 2012). Altogether 49 lines were genotyped by combination of segregation analysis and T-DNA PCR with primers listed in Supplemental Table 1. From these, 43 lines represented 30 unique gene accessions (Supplemental Table 1). As a wild type control, Columbia (Col-0) ecotype of Arabidopsis thaliana was used.

For genotyping, plants were grown in vitro on MS media supplemented with the corresponding selection. For phenotyping, seeds were sown directly on soil with 50% peat and 50% vermiculite. Trays were covered with plastic wrap and cold stratified at +4 ◦C for three nights, after which they were transferred to the growth chamber (FytoScope, PSI, Czech Rep.). Seven days after stratification (DAS) the seedlings were transferred to their own pots, placed on the analysis trays and sand was added on top of the peat to prevent growth of any green algae. From the full water saturation of the soil, the water content was let to decrease until 70% and was kept at this level through daily weighing and watering. Growth conditions in the Arabidopsis growth chambers were 16 h light/8 h darkness and 22 ◦C. Relative air humidity of the growth chambers was targeted at 60%. The light intensity was set and controlled at 130 µE (MS6610, Mastech, China).

### Genotyping of the Mutant Lines

Homozygous one locus mutant lines were confirmed by segregation analysis and T-DNA specific PCRs. The PCR primers, T-DNA position and line information were summarized in Supplemental Table 1. The transcript levels of the T-DNA targeted genes were verified by quantitative real-time PCR (qPCR) analysis. The sample material for qPCR was harvested from the tissue indicated by Arabidopsis eFP Browser for each gene expression pattern: if the gene expression pattern indicated at least moderate expression in flower parts during floral development, tips of inflorescences with developing and open flowers were pooled from three to five individual plants. If the expression was in the seeds, young to mature siliques were pooled. If no expression data via eFP was found, new leaves were pooled with developing and open flowers. Three to four biological replications were harvested for each RNA preparation. RNA was extracted using InviTrap R Spin Plant RNA Kit (STRATEC Molecular), complementary DNA was prepared with SuperScript R IV Reverse Transcriptase (Thermo Fisher Scientific), and the qPCRs were performed using Roche Lightcycler R 480 Instrument II (Roche Diagnostics) using LightCycler R 480 SYBR Green I Master (Roche Diagnostics) with primers listed in Supplemental Table 1. Primers were primarily designed to locus downstream of the T-DNA. The fold up values (mutant line against Col-0) were calculated using the 2−11CT method according to Livak and Schmittgen (2001). Reference genes used in this study were the most stable Arabidopsis genes according to Czechowski et al. (2005): TIP41 LIKE (AT4G34270, forward: GTGAAAACTGTT GGAGAGAAGCAA, reverse: TCAACTGGATACCCTTT∧CG CA), AP2M (AT5G46630, forward: TTGAAAATTGGAGTAC CGTACCAA, reverse: TCCCTCGTATACATCTGGCCA) and PTB1 (AT3G01150, forward: TTGAAGGAGTGGAATCTCAC G, reverse: ATGTGCGGAAAGCAGATACC). Significance level of the qPCR were set at 0.1–0.5 FC for knock-down; <0.1 FC for knock-out; and >2 FC for up-regulated (Supplemental Table 1).

### High Throughput Plant Phenotyping Platform

The plant phenotyping facility at the University of Helsinki Viikki campus (http://blogs.helsinki.fi/nappi-blog/) was used for the phenotypic characterization of the selected Arabidopsis mutant collection. The plants were imaged daily by overhead CCD camera for RGB images positioned in a PlantScreenTM analysis chamber with automated plant transportation between the imaging, weighing and watering stations. The RGB images were obtained for 20 plants at the time and stored in central database. The images were pre-processed online as described in Awlia et al. (2016) to allow collecting binary and RGB data for each plant. The obtained binary images were used for calculating growth parameters of area and perimeter. The obtained parameters of area, perimeter and the convex hull were then used for automatic online calculations of morphometric rosette parameters including: roundness1, roundness2, isotropy, eccentricity, compactness, Rotational Mass Symmetry (RMS), and Slenderness of Leaves (SOL) (PlantScreenTM analyzer, PSI, Czech R.). To characterize the general morphology of the mutant lines these nine morphological parameters were grouped into four categories based on their type: raw, circularity, symmetry and center distance, and compared over time (**Figure 1**). Raw parameters were represented by area and perimeter of the rosette and they were calculated by counting pixels of a rosette picture and the edge pixels respectively and transformed to millimeters (**Figure 1A**). The parameters of roundness1 and roundness2 and isotropy represented the circular parameters (**Figure 1B**). The parameter roundness describes rosette area in comparison to perfect circle with same perimeter and is affected by leaf slenderness, petiole length and leaf perimeter. For wild type plant, this parameter usually takes values between 0.1 and 0.5 while a perfect circle has value 1. Roundness value tends to decay overtime due to leaf development that at the same time increases the rosette perimeter. Roundness 2 uses rosette convex hull area and perimeter for its computation and for wild type plants this parameter appears to have values between 0.7 and 1.0 following an oscillating pattern with less steep peaks over time (**Figure 1B**). Isotropy uses the area of a drawn polygon on top of the rosette (**Figure 1B**). Thus, isotropy has a behavior similar to roundness 2 over time, but with less steep peaks and decreasing tendency similar to roundness. The eccentricity and RMS were symmetric parameters (**Figure 1C**). Eccentricity describes how

elliptical the plant rosette is, where a value close to 1 corresponds to a rosette with highly sharp elliptical shape, while a value close to 0 describes a circular shape. Wild type rosette shows a high eccentricity peak that decays over time with a second smaller peak by the end of growth, thus the rosette shape fluctuates between a round and an elliptical shape. On the other hand, RMS describes the symmetry of the plant rosette by making a ratio between the non-overlapping rosette convex hull area and a perfect circle of the same area centered in the plant centroid and the overlapping area of both. RMS shows a similar pattern as eccentricity, but with higher absolute values and a sharper peak. Compactness and SOL were based on the center distance (**Figure 1D**). Compactness is the ratio between the rosette area and the rosette convex hull area. This parameter tells about petiole length and leaf blade width. The parameter SOL explains how sharp the leaf blades are, but it is also affected by the leaf number. SOL was derived from the ratio between squared rosette skeleton and rosette area. Thus, SOL can take values greater than 0 and below 50 in dimensionless units for wild type plants (**Figure 1D**).

### Experimental Design

Ten days old (10 DAS) Arabidopsis plants were subjected to growth and morphological characterization by top view imaging for the following 10 days. One phenotyping round was designed to accommodate a maximum of 960 Arabidopsis plants representing 36 genotypes at a time in six consecutive experimental rounds called F006–F011. The total number of lines analyzed throughout the six rounds was 43. The maximum of 36 genotypes were divided in three batches that were rotating between the growth area and the PlantScreenTM analysis chamber. Each batch consisted of three experimental units of four mutant genotypes randomized with Col-0, each represented by 20 individual plants. One experiment consisted thus of five trays of altogether 100 plants. Each experimental unit had their own Col-0 wild type in randomized block design to normalize for any local differences in the microenvironments of the PlantScreenTM or the growth area. Each line showing any phenotypic responses was analyzed in at least three independent experimental rounds. Lines that did not show differences as compared to the Col-0 wild type were excluded from the subsequent rounds thus resulting in reduced numbers of genotypes included.

### Phenotypic Analysis of Flowering Time and Flower Structures

After the image based growth and morphological measurements of the 20 mutant and 20 Col-0 plants in the PlantScreenTM system, the flowering time parameters leaf numbers at bolting (LAB) and days to bolting (DTB), were manually counted for each of the plant individuals. The number of rosette leaves were counted at the appearance of the flower bud (developmental stage 5.10, Boyes et al., 2001) and the DTB was recorded at the same time. The flowering time phenotypes were observed in two to three independent experimental rounds. Finally, flowers of the main inflorescences were photographed and further dissected for floral organ analysis under stereomicroscope (SteREO Discovery.V20, Zeiss). Microscopic pictures of the inflorescence tips, single flowers, sepals and petals were taken with the attached camera (AxioCam ICc3, Zeiss). The analyzed inflorescences and flowers originated from several independent experiments. Flower developmental stages were determined as in Smyth et al. (1990). To confirm pollen viability pollen grain staining according to the modified Alexander method was performed (Peterson et al., 2010). Anther images were captured using Leica DFC420 C camera attached to an optical microscope.

### Statistical Analysis

The significance of the differences between mutant lines and Col-0 was computed by contrasting two fitted models to the data points using several order polynomials (Mirman, 2014). First, a model was fitted to all data points and then a second model was fitted including the factor genotype (wild type and mutant). These two models were then compared using the Chi square test to determine if the second model explained more variance than the first one beyond the significance (α = 0.05). If the second model was statistically different from the first one, it implied that the compared genotypes were different. These statistical analyses were conducted in R software (https://www.r-project.org/). LAB and DTB analysis were performed by Analysis of Variance using GLM procedure and pairwise comparisons against Col-0 using option Dunnett in the MEANS statement using SAS/STAT© software version 9.4 (SAS Institute Inc., Cary, NC, USA).

### RESULTS

### Curation of the Arabidopsis RING-Type Ubiquitin E3 Ligase Protein Sequences

To re-confirm the published RING-type ubiquitin E3 ligase proteins encoded in Arabidopsis genome the 27,667 Arabidopsis proteins from the latest genome annotation release (ARA11) were scanned for RING domains. Through filtering the signatures in Gene3D, SUPERFAMILY, ProSiteProfiles, SMART, Pfam, ProSitePatterns altogether 509 putative RING domain containing protein sequences were obtained (Supplemental Table 2). RING gene names and descriptions were obtained from Araport using Thalemine tool (Supplemental Table 2). Araport used curated but also automatic gene annotation, therefore many RING domain containing proteins were annotated as RING/U-box protein although they did not contain U-box domain. Similarly, some were annotated as RING/FYVE/PHD zinc finger superfamily proteins. The 509 RING sequences were compared to the previously described RING-type protein sequences (Kosarev et al., 2002; Stone et al., 2005). From the 509 identified RING domain proteins 457 matched with the 490 previously described, thus resulting in 31 non-matching sequences (**Figure 2**). These non-matching sequences were thoroughly analyzed and, 6 of them were found to be merged with other gene models, 10 had no RING domain, 3 were not found in the database, 3 corresponded to pseudogenes, 7 were split and a new locus identifier had been assigned for them, and 2 were transposable elements (Supplemental Table 3A). The 50 additionally identified RING domain proteins were shown to represent diverse RING domains such as, 1 of D type, 4 of

C2 type, 20 of H2, 16 of HC, 2 of S/T, and 7 of V (CH) type, according to the Stone et al. (2005) classification (Supplemental Table 2).

### Differential Gene Expression Data Identifies 122 Flower Related RING Ubiquitin E3 ligases

To associate the RING domain proteins with flowering or flower development two approaches were followed: (1) identifying those with gene expression enhanced or enriched during flower development or in flower organs, and (2) by searching RING proteins interacting with known flower regulators. For the first approach Genevestigator (Hruz et al., 2008) tool was used to rank differentially expressed genes (DEG) of the identified RING genes over Arabidopsis developmental stages and in flower organs relative to their expression in developed rosette (**Figures 3A,B**). In the selected experiments in Genevestigator database, probes were available for 393 of the 509 RING E3 ligases analyzed. The cut off for DEGs was set at 2-fold to be included in the selection resulting in lists of genes of interest for each of the categories. This process was repeated to identify gene expression enrichment at each of the development stages of bolting, young flower, developed flower, and flower and siliques. For the developmental categories altogether 71 DEGs were identified (**Figure 3C**). In addition to the developmental stages, enrichment for shoot apical meristem, sepal, petal, stamen and pistil organs were retrieved and resulted in 109 DEGs (**Figure 3D**). Some of the RING genes were common between these two categories and in total 122 unique RING genes were up-regulated in the flower related processes. The gene identifying AGI codes of these 122 flower related candidates are provided in the Supplemental Table 3B.

For the second approach we identified 6 additional genes of interest through literature study and from interaction networks of CO, COP1, and TOE2 from BioGRID (http://thebiogrid. org/). Based on these interaction screens 5 RING E3 ligases were selected to the study, represented by the following mutant lines; N656705 (AT5G65683), N686069 (AT1G61620), N372291 (AT3G29270), N2037522 and N67002 (AT4G17680), and N742646 (AT2G44410). In addition, a mutant line for COP1, cop1-6, and RED AND FAR-RED INSENSITIVE 2 (RFI2) for which a role in mediating red and far-red light signaling and ubiquitination activity has been shown in vitro, were included (Stone et al., 2005; Chen and Ni, 2006a). This E3 ligase was selected as a candidate since its expression is regulated by circadian clock and rfi2-1 mutant flowers early (Chen and Ni, 2006b). Thus, one mutant allele for RFI2 (N878610) was included in the study. Mutants representing these genes were analyzed together with the flower up-regulated RINGs and were named flower related UPS candidates in the Supplemental Table 1.

### Representative Mutant Collection

For functional characterization of the 122 flower related UPS candidates and those selected based on literature, a mutant collection was obtained from the NASC stock center. The mutants represented lines from CATMA, SAIL, SALK, and GABI-Kat collections (Alonso et al., 2003; Rosso et al., 2003; Schmid et al., 2005; Kleinboelting et al., 2012). Altogether 43 lines were shown to contain T-DNA insertion in one locus, six were doubtful and were omitted from the analysis. To confirm that the T-DNA insertion had interrupted the gene of interest, their altered expression levels were confirmed by qPCR analysis with primers listed in Supplemental Table 1. For 43 accessions representing 30 unique loci from the 122 flower related UPS candidates and the selected candidates a differential gene expression pattern was analyzed. Altogether 19 lines were knock-outs, and 13 knock-down mutants, and for 10 lines upregulation of the gene of interested was observed (Supplemental Table 1). For one line, no differential expression was confirmed and this was excluded from the phenotyping. For 14 lines alleles were available with similar or opposite gene expression patterns.

### Phenotypic Screen of the Mutant Accessions

From the genotypically and qPCR confirmed T-DNA insertion mutant lines, 43 were subjected to phenotypic characterization by top view RGB imaging using the PlantScreenTM system. Image series of each analyzed line were collected daily allowing analyzing the growth and changes in morphology over time. For scoring those lines showing phenotypes, we fitted general additive models (GAM) to each parameter of each analyzed lines (data not shown). Most of the lines showed no differences to their corresponding Col-0 controls. However, three lines were consistently different across the experiments compared to Col-0 in both growth and rosette morphology: csu1-4 (cop1 suppressor 1-4, N686069), sinal7-2 (seven in absentia like 7- 2, N833574) (Peralta et al., 2016) and rha1a-1 (ring-h2 finger a1a-1, N2045046) (**Table 1**, Supplemental Tables 1, 4). The csu1- 4 mutant rosette was clearly smaller than Col-0 and showed a

yellowish coloration (**Figure 4**). The mutant line rha1a-1 seemed to have smaller leaves than Col-0, however, at the end of growth it appeared to have more leaves that resulted in similar final rosette area as compared to Col-0. This line also had shorter petioles and leaf serration. The third line sinal7-2 rosette was clearly larger than Col-0 but did not show major differences in color, shape or number of leaves (**Figure 4**).

To further analyze these three lines, mixed non-linear models were fitted to their data using several order polynomials for parametric analysis of the models. This analysis confirmed the earlier observations of significant changes in growth and development for these lines over time (**Table 1**). Line csu1-4 showed slower growth, reduced rosette area and perimeter compared to Col-0 along the complete measured period (**Figures 5A,D**). For line rha1a-1 the rosette area was very similar to Col-0 being, however, slightly but significantly larger over time probably due to its higher number of leaves (**Figure 5B**, **Table 1**). Although the differences between rha1a-1 and Col-0 were small the statistical model was able to capture those. Conversely, sinal7- 2 showed both area and perimeter larger than Col-0 indicating more vigorous growth (**Figures 5C,F**).

### Circularity Parameters

Morphological data for parameters of circularity that include roundness, roundness 2 and isotropy were also evaluated for these lines (**Figure 1**). Line csu1-4 showed increased roundness over the total period analyzed in comparison to Col-0 (**Figures 6A–C**, **Table 1**). However, csu1-4 roundness curve had similar pattern to Col-0 but shifted to the right (**Figure 6A**). Similar situation was observed for sinal7-2, where the roundness curve shape was almost identical to Col-0 but in this case was shifted to the left, showing lower roundness along the total time period (**Figure 6C**, **Table 1**). Roundness curve of rha1a-1 was neither shifted nor similar to Col-0 curve. This line showed a lower roundness than Col-0 at the beginning of the analysis,

days after stratification.

TABLE 1 | Polynomial order and their respective Chi square probability from ANOVA test for each parameter used in this study.

\*Comparison was performed using an ANOVA test between a base model and a model including the genetic background as factor. Significance codes: \*\*\*p < 0.001; \*\*p < 0.01; \*p < 0.05.

Base model = Parameter ∼ polynomial of Day + Random factor Day and Plant ID.

Model = Parameter ∼ polynomial of Day \* genetic background (Col-0 or knockout line) + Random factor Day and Plant ID.

reaching a stabilization point around 16 DAS (**Table 1**). For Col-0 plants roundness continued decreasing until it become lower than rha1a-1 (**Figure 6B**).

Line csu1-4 showed a similar roundness 2 pattern as Col-0 that is shifted to the right by approximately 2 days (**Figures 6D–F**). Line rha1a-1, showed an oscillating pattern too, however, its roundness 2 values were constantly close to 0.9 with less steep peaks than Col-0, presenting the highest differences between days 12 and 16 (**Figure 6E**, **Table 1**). Similarly, to line csu1- 4, line sinal7-2 presented an oscillating pattern very similar to Col-0, however, this time the curve had shifted to the left by approximately 1 day (**Figure 6F**).

Isotropy showed similar results as roundness and roundness 2, where line csu1-4 and sinal7-2 had similar oscillating pattern as Col-0, but csu1-4 curve is shifted to the right, while the curve for sinal7-2 is shifted to the left (**Figures 6G–I**). Line rha1a-1 showed a constant high isotropy value decreasing over time until reaching Col-0 pattern by day 23 (**Figure 6H**, **Table 1**).

### Symmetry Parameters

The morphological parameters describing symmetry were eccentricity and rotational mass symmetry (RMS) (**Figure 1**). For eccentricity, line csu1-4 showed a similar pattern as Col-0 plants with a large and a small eccentricity peak, but shifted to the right (**Figure 7A**). Line rha1a-1 presented no shift in its curve, but it showed a rather flat peak around days 11 and 15, remaining lower than Col-0 until the end of the analysis (**Figure 7B**, **Table 1**). This result shows that rha1a-1 is less eccentric than Col-0 along the complete analysis. Line sinal7-2 showed also a similar pattern to Col-0 plants with two eccentric peaks, but slightly shifted to the left (**Figure 7C**). For RMS line csu1-4 showed similar pattern as Col-0 plants, but shifted again to the right about 1 day for the highest peak and remained higher than Col-0 in the last days of the analysis (**Figure 7D**). On the other hand, rha1a-1 presented no shift in its curve, but it showed a decrease in the peak around days 11 and 15, decaying faster and remaining lower than Col-0 plants (**Figure 7E**, **Table 1**). Like in eccentricity, sinal7-2 was almost indistinguishable from the Col-0 plants, except for a slight shift to the left captured by the model (**Figure 7F**).

### Center Distance Parameters

The last two morphological parameters analyzed were compactness and slenderness of the leaves (SOL), which

were based on the center distance (**Figure 1**). Here line csu1-4 showed a decay of compactness overtime in a similar way to Col-0 plants, but its curve was shifted to the right (**Figure 8A**). Lines csu1-4 and sinal7-2 presented quite normal compactness curves, while for rha1a-1 the pattern that was less compact than Col-0 plants at the beginning of the analyzed period (**Figure 8B**, **Table 1**). The compactness later rises above Col-0, showing higher compactness values. Like for the previously described parameters, sinal7-2 compactness curve showed slightly lower values than Col-0, except for the last 2 days where Col-0 plants reached sinal7-2 compactness (**Figure 8C**).

Line csu1-4 showed lower SOL values than Col-0, while rha1a-1 and sinal7-2 showed higher SOL values than Col-0 (**Figures 8E,F**, **Table 1**). The main differences in SOL could be observed during the exponential growing phase of the rosette and reaching a plateau at the end of the analyzed period where the differences to Col-0 plants become insignificant (**Figures 8D–F**).

### Flowering Time Phenotypes

Flowering time mutants identified in the screen represented both with reduced and increased leaf numbers at bolting (**Table 2**). Line csu1-4 (AT1G61620) was clearly early-flowering in both experimental replications. AT5G63970, a putative forkhead box protein, mutant line was early flowering in one of two experimental replications. SBP (S-ribonuclease binding protein) family protein (AT4G17680) was late flowering in both experimental replications. As already shown by others, cop1-6 mutant was early flowering in both LAB (7) and DTB (22). In most of the mutant lines, LAB did not differ from Col-0 in all experimental replications, but the trend was observed in both or all. LAB or DTB of sinal7- 2 did not differ from Col-0 in either of the experimental replications.

### Mutation in SINAL7 Causes Flower Growth Phenotypes

Flower morphology of the analyzed mutants was observed under stereomicroscope. The mutant line sinal7-2 was found to produce flower buds of abnormal shape, characterized by presence of cavities in the bud tips (**Figures 9A,B**). These openings were present at one or both sides of the affected buds and were caused by tips of the lateral sepals bending inwards (**Figures 9E,F**). Also medial sepals frequently showed altered morphology: their tips covered the buds to a lesser extent than in Columbia, resulting in their "blunt" appearance. Whereas these phenotypes were present in all 18 analyzed inflorescences of sinal7-2 plants, regardless of the plant age only two out of 13 analyzed wild type inflorescences showed similar sepal features, restricted to the first six flowers on the main stems. Scoring flowers stage late 12–15 (located between

positions 1st and 20th on the main inflorescences) revealed that in 54% of the mutant flowers (43/80) at least one lateral sepal tip was bent inwards—as compared to 6/50, i.e., 12% in Col-0 (the analyzed flowers came from 13 to 9 individual plants, respectively).

Dissecting flower buds at the end of stage 12 revealed that the occurrence of ingrown lateral sepal tips was accompanied by petal wrinkling, as the sepal shape interfered with elongation of the petals (**Figures 9G,H**). Indeed, in some of the mature flowers with bent lateral sepal tips, the petal blades remained wrinkled; in several cases also pistil or stamen shape was affected (**Figures 9C,D,I,J**).

SINAL7 has been shown to mediate ubiquitination of glyceraldehyde-3-phosphate dehydrogenase 1 (GAPC1) enzyme in vitro and to affect its enzymatic activity and subcellular localization in Arabidopsis (Peralta et al., 2016). In plants lacking GAPC1 male sterility was observed (Rius et al., 2008). To investigate whether deficiency of SINAL7 impairs male fertility in the sinal7-2 mutant, pollen viability was inspected according to the modified Alexander method (Peterson et al., 2010). Anthers of 12 mutant and 11 Col-0 flowers in the developmental stages late 12 and 13 were stained (early and late flowers, originating from at least five individual plants per line). However, no difference between the mutant and Col-0 pollen was observed:

anthers of both lines contained almost exclusively viable pollen grains (**Figures 9K,L**).

### DISCUSSION

Genomic knowledge in both model plants and crops is expanding at a fast pace. However, translating the knowledge from sequence to function and thereby from models to applications is hampered by bottlenecks in screening for the phenotypes associated with the genotypes. Here, we set out to conduct a reverse genetic approach (Bolle et al., 2011) by defining a proportion of the RING type ubiquitin E3 ligases to the developmental processes of flowering time control or flower development. To this end, the RING type ubiquitin E3 ligases were first curated in the most recent Arabidopsis genome annotation (ARA11) that had been improved, e.g., by the next generation sequencing techniques (Krishnakumar et al., 2015). Thereby, many gene models had indeed become obsolete, split, merged or their original sequence had changed. We also found that in the annotations there are a considerable number of RING domain containing proteins annotated as RING/U-box genes. RING and U-box share similar functions and are structurally and functionally similar, both are ubiquitin E3 ligases that work as scaffolds between the ubiquitin E2 conjugase and substrate. However, at the amino acid residual level RING and U-box domains are significantly different; in the RING domain the arrangement of cysteines and histidines mediate binding of two zinc ions to stabilize the RING domain, while the U-box domains are stabilized by a set of hydrogen bonds and salt bridges (Wiborg et al., 2008).

Recent studies have revealed complex molecular networks that include ubiquitin E3 ligases in regulation of flowering (Lazaro et al., 2012; Peng et al., 2013; Xia et al., 2013). To start defining the genomic flower related Ubiquitin Proteasome System of RING E3 ligases, we first verified the gene expression patterns of the curated RING genes. RING E3 ligases work at protein level but are likely to be transcriptionally directed to their relevant tissues. From the 509 RING genes, 122 were indeed associated with flowering with enrichment of gene expression in flower organs or during flowering. This observation prompted us to obtain a representative mutant collection for phenotypic evaluation.

To screen for phenotypes associated with the mutant collection, an automated plant phenotyping facility was utilized. To facilitate a phenotypic screening of a large Arabidopsis mutant collection a phenomics workflow was established to analyze simultaneously up to 36 genotypes in a PlantScreenTM imaging system installed at the Viikki campus of the University of Helsinki (http://blogs.helsinki.fi/nappi-blog/). Although T-DNA insertion knock-out mutants do not always impair gene function (Bolle et al., 2011), the high-throughput phenomics screen of altogether 43 genotypes singled out three mutant lines with clear growth and morphology phenotypes, three mutant lines with flowering time phenotypes and only one with altered flower structure.

For the Arabidopsis growth assessment, we analyzed rosette growth from 10 to 20 DAS. The analysis of such longitudinal data is challenging and demands automated statistical analysis and modeling steps. The rosette growth normally follows a sigmoid pattern showing a lag phase represented by slow growth around the first 10 days, accelerating in the middle and slowing down when getting close to the transition from the vegetative to reproductive phase. The best way to model data with sigmoid behavior is by fitting a three parameter logistic regression (3PL) to explain the three stages (Paine et al., 2012; Tessmer et al., 2013; Neilson et al., 2015). However, our analysis time window captured only the lag and the exponential phases, so a 3PL model was not suitable for our data. Therefore, we used polynomials for more flexibility and a better explanation of the data for all the parameters. This was particularly useful for the initial screening of the data of the tens of lines for the complex parameters like roundness, roundness 2, isotropy, compactness and RMS.

Typically, the parameters of roundness 2, isotropy and RMS increase and decrease over time. This behavior is due to the natural cycle of leaf initiation and expansion. At the beginning when the two first true leaves are developed, the rosette has an elliptical shape that becomes more circular when the leaves 3 and 4 appear and start to expand. Because leaves 3 and 4 keep on expanding, while the leaves 1 and 2 have already stopped expanding, the rosette takes an elliptical shape around day 12 (**Figure 4**). This process is repeated each time two new leaves develop and expand, explaining the oscillating behavior of these parameters. The steepness of each peak decrease over time because previously generated leaves expand making the rosette more circular. Thus, recording fluctuations in these parameters allows establishing the developmental timing of leaf initiation and expansion.

Here, three lines showed consistently significant differences in growth and morphology compared to the wild type Col-0. The mutant lines csu1-4 and sinal7-2 showed similar growth curve shapes as Col-0, but shifted to the left or right, respectively, for all morphological parameters. This behavior was explained by their speed of growth over time. If two lines differ in their growth rate but were analyzed only on one particular day after germination, they could show high differences in morphological parameters. Therefore, longitudinal time course analysis of Arabidopsis rosette growth and shape became compulsory for making accurate conclusions about the effect of a mutation also on morphology. On the contrary, the rha1a-1 mutant did not show major differences in growth, but did for morphology. The increased number and serration of rosette leaves in rha1a-1 rendered the rosette perimeter and the skeleton longer, thereby, reducing the roundness and increasing SOL during all time points (**Figures 6B**, **8E**). Furthermore, the increased number of leaves of rha1a-1 prevented its rosette from taking overly elliptical shape, keeping it more circular than Col-0 plants over time (**Figure 4**). This characteristic was translated in

#### TABLE 2 | Number of leaves and number of days to bolting in Arabidopsis mutant lines grown in LDs.

Pairwise comparisons were performed against corresponding Col-0-line using Dunnett's test.

\*Indicates statistically significant difference (α = 0.05).

N = 19–20 in each row.

higher roundness 2, isotropy, compactness and lower eccentricity and RMS (**Figures 6E–H**, **7B–E**, **8B**). Thus, the morphological parameters can be used not only to record developmental timing but also to explain the plant architecture in a numeric manner.

The line showing an early flowering time phenotype was COP1 SUPPRESSOR1 (CSU1).csu1-4 plants flowered three to six leaves earlier than Col-0 grown under LDs (**Table 2**). In addition to early flowering, csu1-4 plants showed vegetative phenotypes: plants were smaller than Col-0 (**Figure 5**), the eccentricity, RMS and roundness2 development started later than Col-0 (**Figures 6**, **7**), and SOL was smaller than in Col-0 (**Figure 8**). CSU1 has been shown to negatively regulate hypocotyls length in the dark, via ubiquitination of COP1 and repression of SPA1 (Xu et al., 2014). Our results indicate that CSU1 may regulate both vegetative and generative development. The line showing a late flowering phenotype, SBP family protein (AT4G17680), flowered one to two leaves later than Col-0 (**Table 2**). This gene was selected for the phenotypic analysis based on its interaction with TOE2. toe2 is late flowering, and toe1 toe2 double mutant represses FT expression (Zhai et al., 2015). Our results suggest that this SBP family protein could be involved in regulation of flowering time possibly through TOE2. Some SBP family members are known to regulate flowering time. Four SBP proteins, BOTRYTIS SUSCEPTIBLE1 INTERACTOR (BOI) and its three homologous repress flowering by repressing FT expression in a CO dependent manner and a CO independent manner via DELLA proteins (Nguyen et al., 2015). This evidence suggests that there might be a connection between SBP proteins and flowering time control.

In the sinal7-2 mutant defects in flower morphology were observed. SINAL7 has been shown to ubiquitinate

FIGURE 9 | Flower phenotypes of the sinal7-2 mutant. Flower developmental stages assigned according to Smyth et al. (1990). Scale bars: 1 mm (A–J) and 100µm (K,L). (A,B) Representative inflorescences of Col-0 (A) and sinal7-2 (B). All flowers and siliques older than stage 12 have been removed. Mutant flower buds contain cavities beneath the bud tip (indicated with white arrows). (C,D) Petals of a Col-0 (C) and a sinal7-2 (D) flower at stage 15. White arrows pointing at the wrinkled mutant petals. (E,F) Adaxial surface of the sepals from a Col-0 (E) and a sinal7-2 (F) flower at stage 15. White arrow pointing at the bending lateral sepal tip of sinal7-2. (G,H) Late stage 12 flower buds of Col-0 (G) and sinal7-2 (H). The medial sepals have been removed to reveal the elongating and wrinkling petals blocked by the ingrown lateral sepals of the mutant. (I,J) Col-0 (I) and sinal7-2 (J) flowers stage 15. (K,L) Representative anthers from Col-0 (K) and sinal7-2 (L) flowers stage 12–13 stained for pollen viability.

glyceraldehyde-3-phosphate dehydrogenase 1 (GAPC1) and to regulate its enzymatic activity and movement to nucleus (Peralta et al., 2016). GAPC1 plays a role in glycolysis, thus regulating carbon metabolism, and it has been also associated with cytoskeleton and mitochondria (Giegé et al., 2003; Anderson et al., 2004). SINAL7 gene is differentially expressed in the gapc1 knockout mutant as well as in u-ATP9 plants with mitochondrial dysfunction (Rius et al., 2008; Busi et al., 2011). Although both gapc1 and u-ATP9 lines showed defects in male fertility (Gómez-Casati et al., 2002; Rius et al., 2008), we did not observe increased number of aborted pollen grains in the sinal7-2 mutant (**Figures 9K,L**), suggesting that SINAL7 mediated GAPC1 regulation does not impact pollen maturation. Although we have not tested if the sinal7-2 mutation influences pollen germination and pollen tube growth, fertility of the mutant did not seem to be strongly compromised. Instead, we observed defects in sinal7-2 flower morphology - cavities in flower buds and wrinkled petals. Sepal curvature is controlled by giant cells in the abaxial epidermis, in which cell expansion is connected with endoreduplication (Roeder et al., 2010, 2012). A couple of mutants have been identified in which reduction of giant cells was accompanied by their sepals bending inwards. Closer examination of sinal7-2 sepal epidermis will show whether the observed bent sepal tips and resulting flower bud cavities (**Figures 9A,B,E–H**) originate from endoreduplication defects, which could suggest a novel role for the SINAL7 protein. Other flower phenotypes of the mutant—wrinkling of petals as well as bending of stamens and pistils (**Figures 9C,D,G–J**)—seem to be a direct consequence of the abnormal shape of sepals posing an obstacle for the developing floral organs during their growth and release from the buds. Nevertheless, at this point it cannot be ruled out that the SINAL7 ubiquitin E3 ligase could be involved in the development of these flower organs in other ways.

Here we showed that automated, imaging based phenotyping platform is an efficient tool to overcome the limiting factors of manual and visual phenotypic measurements of large plant collections. Imaging based platforms also allow deep resolution

### REFERENCES

of the phenotypes and thereby more precise association with the genotypes. Furthermore, the automated plant management and transportation to imaging, facilitates time course experiments. Thereby, recording longitudinal numeric values indicating changes in rosette size and morphology can be utilized in developmental timing of plant growth and development. Here the customized solution of the PSI PlantScreenTM system by top view CCD camera in combination with online data processing was used for high throughput phenotyping of an Arabidopsis mutant collection. The obtained resolution and high throughput, whereby hundreds of plants can be analyzed in the time that normally a handful would be analyzed, is an obvious advantage.

### AUTHOR CONTRIBUTIONS

MP conducted the genomic screens, performed the phenomic screens and statistical analysis of these data. KM performed the QPCR analysis and participated in phenotyping and statistical data analysis, FW designed and revised genotyping of the mutant collection. MB designed and revised the flower phenotype analysis. EC participated in the phenotyping assays and conducted flowering time experiments. KH designed the project as a whole, approved the data and wrote the manuscript.

### FUNDING

This project was supported by The Academy of Finland (Suomen Akatemia #283138, #256094, #250972) and Becas Chile from Comisión Nacional de Investigación Científica y Tecnológica (CONICYT) Chile.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017. 00416/full#supplementary-material

to reduce CONSTANS expression and are essential for a photoperiodic flowering response. Dev. Cell 17, 75–86. doi: 10.1016/j.devcel.2009. 06.015

to regulate seed and organ size in Arabidopsis. Plant Cell 25, 3347–3359. doi: 10.1105/tpc.113.115063

**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 © 2017 Pavicic, Mouhu, Wang, Bilicka, Chovanˇcek and Himanen. 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) or licensor 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.

# Automated Method to Determine Two Critical Growth Stages of Wheat: Heading and Flowering

Pouria Sadeghi-Tehran\*, Kasra Sabermanesh, Nicolas Virlet and Malcolm J. Hawkesford\*

*Department of Plant Biology and Crop Sciences, Rothamsted Research, Harpenden, UK*

Recording growth stage information is an important aspect of precision agriculture, crop breeding and phenotyping. In practice, crop growth stage is still primarily monitored by-eye, which is not only laborious and time-consuming, but also subjective and error-prone. The application of computer vision on digital images offers a high-throughput and non-invasive alternative to manual observations and its use in agriculture and high-throughput phenotyping is increasing. This paper presents an automated method to detect wheat heading and flowering stages, which uses the application of computer vision on digital images. The bag-of-visual-word technique is used to identify the growth stage during heading and flowering within digital images. Scale invariant feature transformation feature extraction technique is used for lower level feature extraction; subsequently, local linear constraint coding and spatial pyramid matching are developed in the mid-level representation stage. At the end, support vector machine classification is used to train and test the data samples. The method outperformed existing algorithms, having yielded 95.24, 97.79, 99.59% at early, medium and late stages of heading, respectively and 85.45% accuracy for flowering detection. The results also illustrate that the proposed method is robust enough to handle complex environmental changes (illumination, occlusion). Although the proposed method is applied only on identifying growth stage in wheat, there is potential for application to other crops and categorization concepts, such as disease classification.

Keywords: image categorization, computer vision in agriculture, automated field phenotyping, automated growth stage observation, Field Scanalyzer, wheat heading stage, wheat flowering time

### 1. INTRODUCTION

An estimated doubling in required crop production is projected by 2,050 in order to meet the demand of the rapid growth human population (Tilman et al., 2011). To achieve this, an approximate 38% increase over current increases in annual crop production rates is required, and on not much more arable land. Further concerns exist around not only achieving this target in a changing climate, but also achieving it sustainably, whereby reducing agricultural inputs to reduce the environmental degradation caused by our agricultural footprint (Tester and Langridge, 2010). With wheat providing >20% of the worlds calorie and protein intake (Braun et al., 2010), the requirement to increase yield and production is widely recognized.

Breeding and precision agriculture, including information-based management of agricultural systems, are fundamental for achieving sustainable increases in wheat productivity and production.

#### Edited by:

*John Doonan, Aberystwyth University, UK*

#### Reviewed by:

*Yuhui Chen, Samuel Roberts Noble Foundation, USA Andrew French, University of Nottingham, UK*

#### \*Correspondence:

*Pouria Sadeghi-Tehran pouria.sadeghi-tehran@ rothamsted.ac.uk Malcolm J. Hawkesford malcolm.hawkesford@ rothamsted.ac.uk*

#### Specialty section:

*This article was submitted to Technical Advances in Plant Science, a section of the journal Frontiers in Plant Science*

> Received: *28 September 2016* Accepted: *09 February 2017* Published: *27 February 2017*

#### Citation:

*Sadeghi-Tehran P, Sabermanesh K, Virlet N and Hawkesford MJ (2017) Automated Method to Determine Two Critical Growth Stages of Wheat: Heading and Flowering. Front. Plant Sci. 8:252. doi: 10.3389/fpls.2017.00252*

**200**

One component critical to both crop breeding and precision agriculture is the monitoring of developmental growth stages, as (i) it helps crop producers understand which phases of wheat development are most vulnerable to biotic and environmental stresses, and (ii) supports precision agriculture by helping making informed-decisions around which treatment should be applied, to what location and when to apply it. Two critical growth stages monitored in crops, including wheat, are heading date and flowering time, as cultivars with appropriate heading time to their target environment and life cycle duration will help maximize yield potential (Snape et al., 2001; Zhang et al., 2008).

The monitoring of heading and flowering stages are still primarily performed by human eye, which is labor-intensive and time-consuming, as these observations need to be performed on up to thousands of cultivars/varieties on a daily or bi-daily basis, given the importance in catching the starting date of these growth stages. Given that manual growth stage monitoring is also subjective, different observers may likely perceive the growth stage of the same plot differently, which introduces human-error into obtained data.

Computer vision offers an effective alternative for growth stage monitoring because of its low-cost (relative to man-hours invested in to manual observations) and the requirement for minimal human intervention. Computer vision has facilitated automation in high-throughput phenotyping, as well as areas of agriculture, such as disease detection (Pourreza et al., 2015), weed identification (Guerrero et al., 2012) and quality control (Alahi et al., 2012; Valiente-González et al., 2014). Despite the efforts of computer vision specialists over the past decades, developing reliable image-based model to identify and categorize images based on visual information is still difficult to achieve and remains an unsolved problem in the computer vision community. The visual recognition of object categories is a natural and trivial task for humans. Humans can recognize objects effortlessly even with changes in an object's appearance, such as viewing direction or a shadow being cast across the object. On the other hand, in computer vision it can be a challenging task to achieve such level of performance due to the difficulties inherent in the problem. Images are quite abstract and subjected to illumination, scale, deformation, background clutter, etc. Moreover, in computer vision, teaching a machine to distinguish and categorize objects is all about teaching it which differences in the image is matters and which don't, by scanning through diverse datasets, which is a computationally exhaustive process.

Computer vision has shown promise in detecting growth stages of crops. For seedling emergence, color segmentation approaches have been applied in maize (Yihang et al., 2014) and oilseed rape (Yu et al., 2013), using images acquired from a digital camera. Some approaches for observing later growth stages, such as heading date and flowering stage, have also been developed. Zhu et al. (2016), developed a method to detect wheat heading stage from RGB images using a two-step coarse-to-fine detection approach. For flowering stages, Guo et al. (2015) used object-recognition to detect flowering stages from rice panicles. Although the approaches by Zhu et al. (2016) and Guo et al. (2015) were effective on a single variety, within small patches of whole canopies, applications that are more versatile and that also can be applied on different varieties on larger scale canopies are required.

This study utilizes a novel visual-based approach to monitor heading and flowering stage of field grown wheat, through the automated learning of the visual consistency between classes of canopy images, in order to identify the critical growing stages of wheat (e.g., whether ears are emerging in canopies). This method searches through an image database to identify and retrieve images containing emerged wheat ears and ears at flowering stages. This visual-based approach is:

### 2. MATERIALS AND METHODOLOGY

The introduced technique is performed in four main steps (**Figure 1**):

Bag of Visual Words (BoVW) proved to be the leading strategy in computer vision applications such as image retrieval and image categorization (Csurka et al., 2004); thus, it is being opted for the presented work. Categorizing digital images, embarks on extracting features and creating a visual vocabulary for the given dataset. It comprises of following states:

However, in this study, several steps are integrated in the process to improve the overall performance compared to Csurka et al. (2004) described in Section 2.3. Our method treats canopy images acquired automatically in the field as a collection of unordered appearance descriptors extracted from local patches;

**Abbreviations:** BoVW, bag of visual words; SVM, support vector machine; SPM, spatial pyramid matching; RBF, radius basis kernel; LLC, local linear constraint; SIFT, scale invariant feature transform; DoG, difference-of-gaussian; LoG, laplacian-of-gaussian; PCA, principle component analysis; SURF, speeded up robust features; KNN, k-nearest-neighbors, UAV, unmanned aerial vehicle.

then, quantizes them into discrete visual words. Each image is defined by a feature vector listing the number of regions which belongs to each cluster and are later used to train a classifier. In addition, the location information is taken into account which is one of the important factors in object recognition scenarios. In the final step, a linear Support Vector Machine (SVM) classifier is used to determine pre-defined classes (e.g., ear emergence, flowering). The experimental results show that the introduced method is capable of automatically identifying key wheat growing stage with high accuracy and efficiency (Section 3).

## 2.1. Field Experiment and Image Acquisition

Six wheat cultivars (Triticum aestivum L. cv. Avalon, Cadenza, Crusoe, Gatsby, Soissons and Maris Widgeon) were grown in the field at Rothamsted Research, Harpenden, UK, sown in Autumn 2015 and maturing in 2016. These cultivars were selected as they had different properties visible to the naked-eye (awns/no awns, differing wax properties, straight/floppy leaves, different ear morphology) (**Figure 2**). All cultivars were sown 20 October 2015, at a planting density of 350 plants/m2 . Nitrogen (N) treatments were applied as ammonium nitrate in the spring, at

rates of 0 kg ha−1 (residual soil N; N1), 100 kg ha−1 (N2) and 200 kg ha−1 (N3) (**Figure 4**).

Widgeon) used, at growth stage Z5.9.