In cultivated olive (Olea europaea subsp. europaea var. europaea), the cross breeding activities have been delayed by the particularly long generation time (Santos-Antunes et al., 2005), the extended juvenile phase, the high demanding nursery practices, such as the forcing of seedling growth (Rugini et al., 2016) and the time and space needed for plant growing (Picheny et al., 2017). In olive, breeding programs last about 30 years on average (Lavee et al., 2014; Rallo et al., 2016) and have been limited to the empirical selection of a few sporadic intraspecific crosses (Rallo et al., 2008), or to clonal selection (Manaï et al., 2007; Gomes et al., 2008; Trapero et al., 2013; Mousavi et al., 2019), while the timing for the selection of new cultivars in other fruit crops has been greatly reduced, also by the application of new efficient genomic tools (Biscarini et al., 2017; Laurens et al., 2018; Cai et al., 2019). However, the importance of olive cultivation at worldwide level and the new challenges posed by the ongoing climate change, are leading to an ever increasing demand for new cultivars (Gutierrez et al., 2009; Urban, 2015; Bosso et al., 2016).

One of the current limitations to olive productivity is represented by its complex self- and inter-incompatibility system (Saumitou-Laprade et al., 2017a; Alagna et al., 2019), a barrier that may seriously curb yield and restrict the varietal choice for planting to only a few inter-compatible or self-fertile varieties. Fruit set deficiencies due to ineffective pollination are generally underestimated by the farmers, however, it has been demonstrated that supplemental pollination may significantly increase olive production (Ayerza and Coates, 2004), indicating the importance of an effective pollination design of olive orchards.

In olive as well as in other species of the Oleaceae family, such as Phillyrea (Phillyrea angustifolia) and ash (Fraxinus excelsior), an incompatibility system known as diallelic self‐incompatibility (DSI) has been described (Saumitou-Laprade et al., 2010; Vernet et al., 2016; Saumitou-Laprade et al., 2018). It has been hypothesized that two alleles at the DSI locus exist in cultivated olives, S and s, with S dominant over s, with only two possible genotypic combinations (Ss and ss), corresponding to G1 and G2 incompatibility groups, respectively, and by stigma test analysis it was never found G1xG1 or G2xG2 compatibility (Saumitou‐Laprade et al., 2017a), where G1xG2 crosses always generate G1:G2 = 1:1 balanced progenies. All olive cultivars seem to be self-incompatible, even if pseudo-self-fertility might occur for some cultivars in particular conditions (Alagna et al., 2019).

Although the DSI mechanism in olive is known, many important aspects remain to be clarified, such as the location of the incompatibility locus on the genome, the identification of candidate genes controlling this trait, and markers closely linked to incompatibility. The availability of such information will allow for a systematic screening of olive cultivars to identify their group of incompatibility through genotyping with linked markers.

Olive is a diploid species (2n = 2x = 46) with a genome size of approximately 1.4 Gb (Cruz et al., 2016), with a mean C-value of 1.59 pg (1.56 Gb), where more than 30% sequences are represented by tandem repeats (Barghini et al., 2014). Up to now, only intraspecific crosses have been used for olive mapping, and very early studies were performed with dominant markers (de la Rosa et al., 2003; Wu et al., 2004; Aabidine et al., 2010; Khadari et al., 2010). Recently, more dense maps have been produced by the use of codominant markers, such as diversity arrays technology (DArT) (Domínguez-García et al., 2012; Atienza et al., 2014), simple-sequence repeat (SSR) (Sadok et al., 2013), and single-nucleotide polymorphism (SNP) markers (Marchese et al., 2016; İpek et al., 2017; Unver et al., 2017).

SNPs are sequence-tagged markers widely used for association and genetic mapping due to their wide distribution along the genome, high-throughput genotyping, and ease to score (Vezzulli et al., 2008; Deulvot et al., 2010; Lou et al., 2017). Molecular markers linked to the traits of interest can be identified through different strategies, such as genetic linkage mapping based on biparental populations (Curtolo et al., 2017; Ji et al., 2018; Zheng et al., 2018; Sapkota et al., 2019), or through genome-wide association studies (GWAS), based on unrelated individuals (Khan et al., 2013; Nicolas et al., 2016; Elsadr et al., 2019). The generation of high-resolution linkage maps, a prerequisite for gene positional cloning, allows the genetic dissection of quantitative trait loci, assists in comparisons of synteny, and provides marker order for anchoring sequence scaffolds or physical contigs to linkage groups. In perennial fruit crops, the availability of markers tightly linked to traits under selection may strongly facilitate the progress of breeding programs (Minamikawa et al., 2018).

Linkage analysis can be performed on multi-generation families derived by the cross, back-cross, or selfing of homozygous or heterozygous lines. In fruit trees, the use of these progenies is hindered by the lack, in most cases, of homozygous genotypes and by the long generation time (Bai et al., 2007; Jaillon et al., 2007). For this reason, full-sib F1 families deriving from inter- or intra-specific varietal crosses of highly heterozygous parents are generally used in tree species (Bartholomé et al., 2015) and linkage analysis is conducted separately for each parent using a two-way pseudo-testcross mapping strategy (Grattapaglia and Sederoff, 1994). Pseudotestcross mapping has been carried out in many fruit crops, such as apricot (Irisarri et al., 2019), peach (Zeballos et al., 2016), oil palm (Bai et al., 2018), clementine (Ollitrault et al., 2012), grapevine (Zhu et al., 2018), apple (Di Pierro et al., 2016), and forest trees, like poplar (Zhigunov et al., 2017) and oak (Bodénès et al., 2016). A number of markers closely linked to important simple or complex traits has been identified with this strategy (Zhu et al., 2015; Zhai et al., 2016). Marker-trait associations, either for qualitative or quantitative trait loci (QTL), should allow to predict in advance the final breeding value of genotypes, permitting to discard unwanted genotypes immediately after seed germination and to only grow individuals that will later display the traits of interest (Edge-Garza et al., 2015; Migicovsky and Myles, 2017).

In the present work, a F1 progeny derived from the cross of two highly heterozygous and completely self-incompatible cultivars, 'Leccino' × 'Dolce Agogia,' respectively belonging to G1 and G2 groups of incompatibility (Saumitou‐Laprade et al., 2017b), was genotyped by ddRAD sequencing technology in order to construct linkage maps. This crossing was performed because parents show different phenotypes for numerous important agronomical traits. In particular, 'Dolce Agogia' is a vigorous cultivar, with an alternate bearing, resistant to some of the most dangerous olive diseases, such as Verticillium dahliae and Spilocaea oleagina (Iannotta and Scalercio, 2012; Arias-Calderón et al., 2015), whereas 'Leccino' has medium vigor and constant bearing, is susceptible or only partially resistant to the pathogens indicated above (Salman, 2017), but it was recently reported as tolerant to Xylella fastidiosa (Giampetruzzi et al., 2016), the most devastating emerging plant pathogen for the Mediterranean agriculture (Suffert et al., 2009). The maps generated in the present work, represent a significant improvement over the previous ones (Marchese et al., 2016; İpek et al., 2017; Unver et al., 2017), because they include a higher number of SNP markers with a uniform genome distribution and can, therefore, serve as saturated maps for trait mapping. We used the 'Leccino' map to identify the location of the DSI locus, to select sequence scaffolds of 'Leccino' that were anchored to the DSI by linked markers, and to develop and validate SNP markers for the incompatibility trait.

#### MATERIALS AND METHODS

#### Plant Material and Deoxyribonucleic Acid Extraction

A full-sib F1 progeny of 229 individuals was generated from the cross between 'Leccino,' as female parent, and 'Dolce Agogia,' as male parent. The progeny, thereafter referred to as Le×DA, is represented by a first set of 16-year-old seedlings (amounting to 155 individuals) and by a second set of 5 years old seedlings (74 individuals) grown in an experimental field, for a total of 229 F1 genotypes. The progeny showed high phenotypic variability, in particular for the length of the juvenile phase, plant vigor, tree habit, fruit bearing, and cutting's rooting ability (Hedayati et al., 2015). DNA was extracted from leaves using the DNeasy Plant Mini Kit (Qiagen).

#### Parentage Analysis of the Progeny

'Leccino,' 'Dolce Agogia,' and the Le×DA progeny were genotyped with a set of SSR markers for confirming the parentage (Baldoni et al., 2009). PCRs were performed as previously reported (Mousavi et al., 2017) and the amplified fragments were separated on an ABI 3130 Genetic Analyzer capillary sequencer (Applied Biosystems, Foster City, CA). Alleles were called using the GeneMapper 3.7 software (Applied Biosystems, Foster City, CA). Parentage analysis was performed using CERVUS version 3.0.7 program (Kalinowski et al., 2007) to sort out seedlings derived from open pollination.

#### Phenotyping for Incompatibility Groups of the Le×DA Progeny

Those offspring that have reached the mature phase and started blooming (109 out of 229, all belonging to the first set of 16-yearold seedlings), were phenotyped for the incompatibility group. Data on 91 individuals were previously reported (Saumitou-Laprade et al., 2017a), all the others were de novo phenotyped following the same protocol previously applied. In order to minimize errors, phenotyping of all individuals was repeated on previously clonally propagated plants grown in a separate field. Before blooming, flowering twigs of each tree were protected from open pollination with blossom bags. At full blooming, bagged twigs were collected and 20 open flowers per tree were used for the analysis. Sepals and petals were removed from 10 flowers and pistils were placed in two separate plates (one for each of the two pollen donors) onto Brewbacker and Kwack medium. The stigmas of five pistils in one plate were pollinated with pollen of 'Leccino' (G1) and the other five with pollen of 'Dolce Agogia' (G2). From the remaining 10 flowers, pollen was collected and used for pollinating 'Leccino' and 'Dolce Agogia' stigmas. The same protocol applied by Saumitou-Laprade et al. (2017a) was used to observe the growth of pollen tubes on the stigmata.

#### Double Digest Restriction Associated Deoxyribonucleic Acid Sequencing of Parents and Progeny

Genomic DNA was digested with SphI and MboI restriction enzymes, following a double digest restriction associated DNA sequencing (ddRADseq), according to the method proposed by Peterson et al. (2012) and modified by Scaglione et al. (2015). Fragments were added to a ligation reaction containing barcoded adapters, pooled, and then fractioned by agarose gelelectrophoresis. DNA in the size range between 350 and 600 bp was purified using a QIAquick Gel Extraction kit (Qiagen, Venlo, Netherlands). Enrichment PCR was performed with PCR primers that incorporate Illumina hybridization/sequencing sites and index sequences for combinatorial multiplexing. Quality, quantity, and reproducibility of libraries were assessed using a Caliper instrument (DNA High Sensitivity chip). Sequencing was carried out on an Illumina HiSeq 2500 instrument, generating 125-bp paired-end reads.

#### Single-Nucleotide Polymorphism Calling

Reads were aligned using Bowtie 2 software (Langmead and Salzberg, 2012) with default parameters against a whole-genome assembly of 'Leccino' (Muleo et al., 2016). The reference consisted of 509,032 scaffolds, with N50 length of 10,037 bp, amounting to a total of 1.429 giga base pairs (Gbp). Alignments with mapping quality < 4 were removed. Segregating sites were identified using the software Stacks (Catchen et al., 2011) and a bounded SNP model with alpha = 0.05 and upper error (epsilon) of 0.1. Genotypes were called with a minimum coverage of eight reads. Heterozygous genotypes were called within a read coverage ratio of 0.15-0.85 for reference and alternate alleles. Segregating loci were retained if genotypes were called in > 150 progeny.

#### Genetic Mapping

Genetic maps were generated with the double pseudo-test-cross approach for each parent using SNP genotypic data and phenotypic binary data for the incompatibility group. Linkage groups were obtained using the R/qtl module of the R statistical package with a logarithm of odds (LOD) threshold > 10 and a recombination rate of 0.20. Markers were first ordered using MSTMAP (Wu et al., 2008) with default parameters. A Perl implementation of the SMOOTH program (van Os et al., 2005) was used to remove errors within haplotypes. Markers were re-ordered using MSTMAP. Marker distances were calculated using the Kosambi function {map distance in centi-Morgan (cM) equals ¼ ln [(1 + 2r)/(1 − 2r)], where r is the observed recombination frequency}. Markers with distorted segregation were identified using a c2 test (a = 0.05).

#### Comparative Genomic Analyses

Scaffolds of 'Leccino' that were anchored to the DSI locus by flanking markers were aligned with the genome assembly of the wild olive (Olea europaea var. sylvestris, GCA\_002742605.1) using (B)LASTZ. Gene annotation was performed using the highest hit in blastp alignments with the NCBI protein database.

#### Development of a Sequence-Tagged Site Marker for Diallelic Self-Incompatibility

PCR primers were designed on the sequence of scaffold\_6030 using the program Primer3 version 4.0 (forward 5'-3': TTTTGGGTGCGAATTGTCCA, reverse 5 ' - 3 ' : AGGCCACTGTATTTCTAACTCG) for the amplification of a 476-bp fragment spanning two adjacent SNPs that co-segregate with DSI. PCRs were performed using 25 ng of template DNA and Q5 High-Fidelity DNA polymerase (New England Biolabs). The thermal profile consisted of 98°C for 30 s, followed by 35 cycles at 98°C for 10 s, 60°C for 20 s, and 72°C for 30 s, and a final extension at 72°C for 2 min. Amplicon size was checked by 1% agarose gel electrophoresis. Amplicons of the expected size showing unique bands were sequenced using the BigDye Terminator v1.1 Cycle Sequencing Kit (Thermo Fisher Scientific) and an ABI PRISM 3130 XL Genetic Analyzer (Applied Biosystems, Foster City, CA). The obtained sequences were aligned using BioEdit 7.1.7 (www.mbio.ncsu.edu/BioEdit/ bioedit.html), to identify polymorphisms. In order to phase single-nucleotide substitutions in heterozygous individuals and to obtain full sequences from individuals carrying heterozygous small indels, amplicons were cloned by using pGEM-T Easy Vector (Promega) and Escherichia coli XL1 blue strain. DNA from different colonies for each genotype was amplified and sequenced as described above.

### RESULTS

#### Genetic Maps

The parentage analysis allowed to validate the origin of 95% of the 241 seedlings for a total of 229 true-to-type F1 genotypes. The remaining 12 individuals were excluded from further analyses. A total of 16,743 segregating loci were identified in the Le×DA progeny by ddRAD sequencing. The parental maps consisted of 23 linkage groups, including 9,737 RAD loci in 'Dolce Agogia' and 7,006 RAD loci in 'Leccino' (Table 1). Linkage groups were numbered consistently with the chromosome numbering of the wild olive genome assembly (Unver et al., 2017). The male parental map ('Dolce Agogia') included 1,829 genetic bins (i.e., positions on the genetic map with a unique segregation pattern) for a total length of 3,538 cM. The female parental map ('Leccino') included 2,311 genetic bins for a total length of 5,680 cM. The average distance between adjacent genetic bins was 2.46 cM in 'Leccino' and 1.93 cM in 'Dolce Agogia'.

TABLE 1 | Statistics of the parental linkage maps. Chromosome number (Chr), restriction associated DNA (RAD) markers, length (cM), and bins are reported for 'Leccino' and 'Dolce Agogia'.

All linkage groups of the 'Leccino' map were consistently longer that those of the 'Dolce Agogia' map, presumably as a result of a consistently higher recombination rate between homologous chromosomes in megaspore mother cells of 'Leccino' than in pollen mother cells of 'Dolce Agogia.' The shorter map obtained in 'Dolce Agogia' was not explained by runs of homozygosity at the chromosome termini compared to 'Leccino,' except for the upper 34 cM of Linkage Group (LG) 16 in 'Leccino,' which lacked segregating sites in 'Dolce Agogia.' The longer map obtained in 'Leccino' was, on the other hand, not explained by the presence of isolated markers with high likelihood of genotypic errors, which usually locally inflate genetic distances due to inconsistent genotype calls compared to neighboring markers. It was noteworthy that the suppression of recombination was observed in the parent with the higher level of heterozygosity, as revealed by the larger number of segregating sites.

#### Structural Variation Between Parental Genomes and Segregation Distortion

We used 2,305 RAD loci segregating from both parents with an allelic status "ABxCD" for pairing and aligning linkage groups between parental maps. The alignment of the two parental maps did not show inter-chromosomal translocations. Two large inversions were detected in LG 1, involving a segment of about 30 cM, and in LG 8, involving a segment of about 12 cM. We did not detect suppression of recombination around each inversion in any of the parental maps and we mapped segregating sites within the inverted regions in each parental map. In both cases, the parental genomes were therefore homozygous for alternative structural variants, but the haplotypes carrying the inversion were ancient enough to have accumulated segregating sites. Minor rearrangements in marker order were identified on LGs 2, 4, 5, 7, and 17. Extended regions with distorted segregation (a = 0.05) were detected on several linkage groups in both maps, with a higher frequency being observed in 'Dolce Agogia' map (Figure 1). Localized distortion was also detected at markers flanking the lower side of the incompatibility locus (see below) with an excess of marker alleles in phase with the incompatibility recessive allele.

### Identification of the Incompatibility Locus

A stigma test was performed on 106 adult Le×DA individuals. The incompatibility group phenotyping assigned 56 individuals to the G1 group, 47 to the G2 group, and discarded 3 individuals due to uncertain phenotype. Segregation of the phenotype followed the expected 1:1 ratio (c2 = 0.79). A self-pollination test of 'Leccino' and 'Dolce Agogia' showed no growth of pollen tubes on the stigmas, confirming their complete self-incompatibility.

The incompatibility locus was mapped as a Mendelian trait in the 'Leccino' map. The locus was located on LG 18 within an interval of 5.4 cM. Within this interval, the RAD markers on scaffold\_6030 (sizing 42,057 bp) co-segregated with the trait in all 103 individuals of the progeny with a reliable incompatibility

group assignment (Figure 2). The upper border of the locus was supported by a recombination event observed in one individual between the incompatibility group and two RAD markers located in scaffold\_63515 (sizing 3,877 bp) and in scaffold\_26872 (sizing 32,989 bp). The lower border of the genetic interval was supported by five crossing-over events observed in the progeny between the incompatibility group trait and a RAD marker in scaffold\_7600 (sizing 34,918 bp) and by even more crossing-overs with a RAD marker in scaffold\_13712 (sizing 37,403 bp). All these scaffolds of 'Leccino' aligned to a region between coordinates 8,500,000 and 9,100,000 of chromosome 18 in the genome assembly of a wild olive (Unver et al., 2017) (Figure 3).

#### Alignment of Diallelic Self-Incompatibility Markers and Scaffolds With the Wild Olive Genome Assembly and Candidate Genes in the Region

The projection of 'Leccino' scaffolds containing the RAD flanking markers from either side of the DSI locus onto the assembly of chromosome 18 identified a seemingly short physical interval of a few dozen Kbs, comprised between scaffolds 13712 and 7600, on one side, and scaffold 26872 on the other side (Figure 3). The region around scaffolds 13712, 7600, 26872 in the wild olive assembly, corresponding to the chromosomal interval between 8.5 and 8.7 Mb, encodes five predicted proteins (Supplementary Table S1, IDs 1 to 5). However, the inconsistent physical position of scaffold\_7600, which contains a RAD marker co-segregating with DSI and aligns outside of the interval defined by the flanking markers, and the split alignment of the initial 14 Kb of scaffold\_7600 on one side of the locus and the remaining 12 Kb of scaffold\_7600 on the opposite strand at the other side of the locus suggest that the physical interval might be substantially longer. A 450-kb inversion, indicated by a double-head arrow in Figure 3, or an assembly error in the wild olive genome between two flanking sequence gaps would reconcile the genetic marker order in the 'Leccino' map and the assembly of scaffold\_7600 and would define a physical interval for the DSI locus in the wild olive assembly from coordinates 8,720,000 to 9,080,000, with inverted orientation (220 Kb-580 Kb relative coordinates in Figure 3). Under this hypothesis, the DSI region in the assembly of wild olive would encompass 12 additional predicted proteins, bringing the total number of candidate genes to 17 (Supplementary Table S1, IDs 6 to 17). RNA-Seq data also showed that other regions outside of the predicted gene models are transcribed in the physical interval of wild olive.

#### Validation of Sequence-Tagged Site Markers Linked to the Incompatibility Group

We amplified and sequenced a PCR fragment (hereafter referred to as Oe-DSI-locus-fragment-A), spanning two SNPs on scaffold\_6030 that showed co-segregation with RAD loci in 'Leccino' map. The sequencing results of 165 genotypes, which included cultivars and progenies, allowed to identify only two

arrow below indicates the region of order inconsistency between 'Leccino' map and the wild olive genome assembly.

polymorphic sites at positions 63 and 283 bp of the Oe-DSIlocus-fragment-A, which define unambiguously the two groups of incompatibility (Table 2 and Supplementary Table S2). In the parental cultivars Leccino and Dolce Agogia, the amplified fragment showed four alleles differentiated by eleven SNPs (Table 3). The inheritance of four haplotypes segregating from 'Leccino' (S-A/s-a) and 'Dolce Agogia' (s-b/s-c) and the linkage of the S-A haplotype with the genetic determinant of the G1 phenotype was confirmed by analyzing the Oe-DSI-locusfragment-A in the Le×DA progeny (Figure 4). Combining the stigma test results with the genotypic data, the dominant allele S-A and three recessive ones, s-a, s-b, s-c, were identified.

When the same fragment was amplified in 57 olive cultivars with known incompatibility group assignment, an additional haplotype was found, the S-B allele, which represents a new dominant allele distinguishing G1 cultivars (Supplementary Table S2). The five haplotypes carried a total of 12 SNPs and 2 small indels (Table 3). In the 20 G1 cultivars, the dominant S-A allele was detected 13 times and S-B was present in apparent homozygous state in six out of seven cases and in one case with the haplotype s-b, detected by cloning. Both dominant alleles were found in combination with s-a, s-b, and s-c recessive alleles. The presence of only one dominant allele in some cultivars could be due to homozygosity or to the inability to amplify the second allele (null allele). Among the 37 G2 cultivars, the most frequent recessive haplotype was s-b, also showing a high percentage of homozygosity (Supplementary Table S2). A skewed geographical distribution was observed for some alleles between Italy and Spain, the two most represented countries in

TABLE 2 | Incompatibility group phenotypes and corresponding haplotype combinations of Oe-DSI-locus-fragment-A observed in the analyzed olive cultivars.

a Genotype at position 63 bp of Oe-DSI-locus-fragment-A.

b Genotype at position 283 bp of Oe-DSI-locus-fragment-A.

\*The second allele could be considered as null allele.

TABLE 3 | Haplotypes of Oe-DSI-locus-fragment-A and position (base pair distance from the forward primer) of the polymorphisms [single-nucleotide polymorphisms (SNPs) and indels] identified. SNP combinations that identify uniquely each incompatibility group are indicated in bold.

our data set, such the S-B allele, mainly present in Spanish cultivars, and s-c, only present within the Italian ones (Supplementary Figure S1).

#### DISCUSSION

Species with a long juvenile phase require many years to reach the first flowering and fruit setting (Purba et al., 2001; Flachowsky et al., 2011; Yang et al., 2016), delaying the generation of crossbreed, F2, and backcross populations and slowing down the process of genetic improvement. In these cases, the availability of genomic tools to assist the selection of new genotypes becomes mandatory in order to guide the choice of parental combinations, to select genotypes carrying the traits of interest, and to introgress specific alleles into new varieties. Markers linked to the traits under selection, identified either through genetic mapping approaches or by means of genomewide association studies, represent the most powerful tools for breeding in woody perennial crops, offering new opportunities to develop early selection strategies and new ways to integrate variation from different sources (Varshney et al., 2005; Montanari et al., 2013; Bink et al., 2014; Kole et al., 2015; Muranty et al., 2015).

The use of high-throughput sequencing technologies is speeding up the genotyping of mapping populations, providing an unprecedented high number of markers that allows the construction of dense genetic maps (Kujur et al., 2015; Liu et al., 2017). In this work, we have used ddRAD sequencing for genotyping 7,006 and 9,737 segregating sites in a biparental population, which largely exceeded the number of genetic bins that can be resolved in a progeny of 229 individuals. Olive tree is known to be a highly heterozygous species with very high levels of nucleotide diversity observed even in cultivated varieties

'Dolce Agogia' (s-b/s-c) allele combinations and segregation of four haplotypes in the Le×DA progenies.

(Gros-Balthazard et al., 2019). Here we show that the species is characterized also by very high levels of recombination as attested by the lengths of the two parental genetic maps. We observed a significant difference in the length between the male Dolce Agogia (3,538 cM) and the female Leccino (5,680 cM) maps. Heterochiasmy, i.e., the presence of different crossover frequencies in male and female meiosis, has frequently been observed in plants, without a fixed trend of higher frequency in male or female meiosis (Lenormand and Dutheil, 2005). In Arabidopsis, for example, a dramatically higher crossing over rate (575 cM vs. 332 cM) is observed in male than in female meiosis (Giraut et al., 2011). Theory would predict that haploid selection determines heterochiasmy, with the sex experiencing more intense selection during the haploid phase showing lower recombination (Lenormand and Dutheil, 2005). In a highly heterozygous and obligately outcrossing species such as olive, we can expect stronger selection among male gametophytes that could explain the observed difference in genetic map length. An alternative explanation could be provided by a high frequency of chromosomal inversions or other chromosomal rearrangements that suppress recombination. The identification of two large inversions for which the two parental varieties are homozygous for alternative alleles makes us believe that in such a highly heterozygous species, each individual accession may be heterozygous for a much larger set of inversions. As 'Dolce Agogia' appears to be more heterozygous than 'Leccino,' the lower recombination frequency observed could be explained by a higher frequency of heterozygous chromosomal rearrangements that suppress recombination.

In this case, the generation of markers did not represent a limiting factor for trait mapping, whose resolution was rather limited by the size of the progeny subject to phenotyping. Phenotyping of reproductive traits requires adult plants and, in the case of DSI group assignment, is also labor intensive. We mapped the DSI locus to a 5.4 cM genetic interval, which is estimated to correspond to a physical distance of approximately 300 Kb in a chromosome-scale assembly of a wild olive genome (Unver et al., 2017). Gene prediction and gene annotation in the wild olive assembly did not provide any obvious functional candidate, except for transcription factors that might be involved in gynoecium development regulation. Two candidate genes encode proteins putatively related to flower development: the ortholog of Arabidopsis STYLISH 1 encoding a binding protein with nuclear localization that promotes formation of stylar and stigmatic tissues and proliferation of stylar xylem (Gomariz-Fernández et al., 2017; Min et al., 2019), and FAR1 related sequence 5, which is expressed in hypocotyls, inflorescences stems and flowers (Li et al., 2017; Ma and Li, 2018). However, the prediction of proteins with uncharacterized function in the same region, the presence of non-annotated transcribed regions, the expected intraspecific presence/absence variation between genomes, make it necessary to proceed with the complete assembly of the two 'Leccino' haplotypes across the entirety of the DSI locus.

While we confirmed the monogenic nature of the DSI system in olive, as postulated by Saumitou-Laprade et al. (2017a), five haplotypes were identified using a sequence-tagged site (STS) marker in the locus and, via linkage mapping and association mapping, we demonstrated that two of them are in phase with the dominant genetic determinant of the G1 incompatibility group. The STS haplotypes were linked to the DSI genetic determinant, but they do not correspond to the DSI alleles. It was noteworthy that the cumulative frequency of the STS haplotypes in phase with the dominant S allele was roughly half the cumulative frequency of haplotypes in phase with the s allele, although the accessions analyzed in this paper are not a natural population. It is possible to speculate that this condition may reflect a balancing selection for maintaining G1 and G2 genotypes in a cultivated population at frequencies that maximize pollination rate.

In the set of cultivars analyzed in this paper, we found only two dominant S-A and S-B alleles that confer the G1 phenotype and no cultivar was observed carrying both dominant alleles. It was also noteworthy that some haplotypes were more frequent in groups of cultivars typical of specific geographic locations. For instance, the allele s-c showed 18% frequency in the sample and was almost exclusively present in varieties of Central Italy, leading to an excess of homozygous genotypes in that geographic area, whereas the S-B haplotype was mostly present in cultivars from the Iberian Peninsula. Variation of haplotype frequencies at DSI-linked markers in geographically distant populations may simply reflect genetic drift or may indicate adaptive evolution changing the frequency of DSI alleles. Further analyses will help to clarify these issues and to identify additional untapped variation at the DSI locus.

From an evolutionary point of view, the maintenance of a homomorphic DSI in a hermaphrodite species like O. europaea is unexpected and remains to be explained. Indeed, the SI systems are susceptible to rapid invasion by new self-incompatibility alleles, which should experience a strong negative frequencydependent advantage. In Oleaceae species in which DSI was identified, such invasion was not observed (Saumitou-Laprade et al., 2010; Vernet et al., 2016; Saumitou‐Laprade et al., 2017a; Saumitou‐Laprade et al., 2017b; Saumitou-Laprade et al., 2018). The homomorphic DSI system seems to "resist" and this is an intriguing finding (Barrett, 2019). In the absence of obvious evolutionary constraints that could prevent the selection of new SI specificities, the molecular constraint appears to be the best candidate. Therefore, the molecular characterization of the S locus region in Olea and its comparison with other Oleaceae species are of strong interest. The present work is opening the way and should contribute efficiently to solve the evolutionary paradox of the stable homomorphic DSI in Oleaceae.

The STS marker identified in the present work represents a new tool for large-scale screening of thousands of olive accessions from their traditional areas of cultivation across the Mediterranean shores and from new growing areas and will allow to determine their incompatibility group. In addition, it represents a starting point for the identification of the genetic determinant of such a peculiar incompatibility system. A comprehensive understanding of the genetic control of DSI can offer great opportunities to characterize cultivars for their incompatibility group, increase the olive production and guide the orchard plantation design with optimal spatial distribution of inter-compatible varieties.

#### CONCLUSIONS

This work provides markers for a fast and reliable genotyping of olive cultivars for their incompatibility group, offering great opportunities to rapidly screen and identify inter-compatible varieties, planning inter-varietal crosses, and reducing the time for seedling selection. It also represents the initial effort for the identification of the genetic determinants of incompatibility, a starting point for understanding the molecular mechanisms underlying the DSI system in olive.

The genetic maps of 'Leccino' and 'Dolce Agogia' will also serve to identify Mendelian loci or QTLs responsible for other important traits that segregate in the Le×DA progeny. Information yet to be gained on the number and location of the genetic determinants of those traits, along with the DSI-linked ddRAD and STS markers developed in this paper, may pave the way to the application of genomics-assisted breeding in olive.

#### DATA AVAILABILITY STATEMENT

Raw reads of ddRAD-Seq have been deposited in Short Read Archive under the BioProject number PRJNA594490. The matrix of genotypic calls has been deposited in the figshare repository with the DOI 10.6084/m9.figshare.11352068.

#### AUTHOR CONTRIBUTIONS

RM, PS-L, PV, MM and LB conceived the study. SP, RM, NC and LB provided the plant material. RM, SM, MR, PS-L, PV and FA

#### REFERENCES

performed the stigma tests. RM, NC, VP, MR, AF, GM, SS, DS, GDG and MM performed the molecular, sequencing and bioinformatic analyses. FB, LB, RM, SM and VP wrote the first draft of the manuscript. PS-L, PV, DS, GDG and MM contributed to the writing and revision of the manuscript. All the authors agreed on the final version of this work.

#### FUNDING

The research was supported by the European Union's Horizon 2020 Research and Innovation Program Marie Sklodowska-Curie - Before Project (Grant Agreement No 645595), by the EU projects "OLIVE4CLIMATE – LIFE" (LIFE15 CCM/IT/ 000141) and by the Rural Development Program of Umbria Region, 2014-2020 – Measure 16.2.1, INNO.V.O. - Development of alternative varieties to face the new challenges of olive growing, SIAN n.<PHONE_NUMBER>5.

#### ACKNOWLEDGMENTS

We thank the Istituto Tecnico Agrario "Ciuffelli" (ISIS) Todi (PG) for hosting the Le×DA progeny in the field and for their technical assistance.

#### SUPPLEMENTARY MATERIAL

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

SUPPLEMENTARY FIGURE S1 | Geographical distribution of DSI alleles in the Italian and Spanish analyzed varieties.

(Prunus avium L.) is associated with domesticated and bred germplasm. Sci. Rep. 9 (1), 5008. doi: 10.1038/s41598-019-41484-8

downy mildew resistance in Vitis aestivalis-derived 'Norton'. Theor. Appl. Genet. 132, 137–147. doi: 10.1007/s00122-018-3203-6

shape traits in Pisum sativum using SLAF sequencing. Front. Genet. 9, 615. doi: 10.3389/fgene.2018.00615

Conflict of Interest: The authors declared 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 © 2020 Mariotti, Fornasiero, Mousavi, Cultrera, Brizioli, Pandolfi, Passeri, Rossi, Magris, Scalabrin, Scaglione, Di Gaspero, Saumitou-Laprade, Vernet, Alagna, Morgante and Baldoni. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) 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.

# Ploidy-Dependent Effects of Light Stress on the Mode of Reproduction in the Ranunculus auricomus Complex (Ranunculaceae)

#### Fuad Bahrul Ulum1,2, Camila Costa Castro1 and Elvira Hörandl 1\*

1 Department of Systematics, Biodiversity and Evolution of Plants, Albrecht-von-Haller Institute for Plant Sciences, University of Göttingen, Göttingen, Germany, 2 Biology Department, Faculty of Mathematics and Sciences, Jember University, Jember, Indonesia

#### Edited by:

Emidio Albertini, University of Perugia, Italy

## Reviewed by:

Juan Pablo Amelio Ortiz, National Council for Scientific and Technical Research (CONICET), Argentina Eric Javier Martínez, Instituto de Botánica del Nordeste (IBONE-CONICET), Argentina

\*Correspondence:

Elvira Hörandl<EMAIL_ADDRESS>
#### Specialty section:

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

Received: 29 November 2019 Accepted: 23 January 2020 Published: 20 February 2020

#### Citation:

Ulum FB, Costa Castro C and Hörandl E (2020) Ploidy-Dependent Effects of Light Stress on the Mode of Reproduction in the Ranunculus auricomus Complex (Ranunculaceae). Front. Plant Sci. 11:104. doi: 10.3389/fpls.2020.00104 Polyploidy in angiosperms is an influential factor to trigger apomixis, the reproduction of asexual seeds. Apomixis is usually facultative, which means that both sexual and apomictic seeds can be formed by the same plant. Environmental abiotic stress, e.g. light stress, can change the frequency of apomixis. Previous work suggested effects of stress treatments on meiosis and megasporogenesis. We hypothesized that polyploidy would alter the stress response and hence reproductive phenotypes of different cytotypes. The main aims of this research were to explore with prolonged photoperiods, whether polyploidy alters proportions of sexual ovule and sexual seed formation under light stress conditions. We used three facultative apomictic, pseudogamous cytotypes of the Ranunculus auricomus complex (diploid, tetraploid, and hexaploid). Stress treatments were applied by extended light periods (16.5 h) and control (10 h) in climate growth chambers. Proportions of apomeiotic vs. meiotic development in the ovule were evaluated with clearing methods, and mode of seed formation was examined by single seed flow cytometric seed screening (ssFCSS). We further studied pollen stainability to understand effects of pollen quality on seed formation. Results revealed that under extended photoperiod, all cytotypes produced significantly more sexual ovules than in the control, with strongest effects on diploids. The stress treatment affected neither the frequency of seed set nor the proportion of sexual seeds nor pollen quality. Successful seed formation appears to be dependent on balanced maternal: paternal genome contributions.  Pollen quality was in hexaploids better than in diploids and tetraploids. These findings confirm our hypothesis that megasporogenesis is triggered by light stress treatments. Comparisons of cytotypes support the hypothesis that ovule development in polyploid plants is less sensitive to prolonged photoperiods and responds to a lesser extent with sexual ovule formation. Polyploids may better buffer environmental stress, which releases the potential for aposporous ovule development from somatic cells, and may facilitate the establishment of apomictic seed formation.

Keywords: apomixis, single seed flow cytometric seed screening, light stress, meiosis, pollen, polyploidy, Ranunculus, seed formation

## INTRODUCTION

Polyploidy is a heritable trait of obtaining more than two sets of chromosomes in the nuclei (Comai, 2005). A polyploid arises either from intraspecific genome duplication (autopolyploidy) or the merging of the genome of distinct species through hybridization and subsequent genome duplication (allopolyploidy) (Grant, 1981). Polyploidy is quite common in flowering plants, estimated to occur in more than 50% of species (Soltis et al., 2015) and is considered as a major factor in plant evolution (Soltis et al., 2014). Even though polyploidy is potentially obstructed by several disadvantages, e.g., disruption effects of structural enlargement of nuclei, side effects of aneuploidy, and epigenetic mutation, it also provides advantages such as heterosis, gene redundancy, and novel gene combinations. Heterosis favors polyploids that are more vigorous than their diploid progenitors, while gene redundancy protects polyploids from the deleterious effect of mutation (Comai, 2005).

Polyploidisation, with higher DNA content, increases the cell size and promotes diversity of the genome, transcriptome, and metabolome. These improvements imply a greater resistance to environmental change (Schoenfelder and Fox, 2015). Several studies reported a better adaptivity of polyploid plants to abiotic stress conditions, such as salt (Chao et al., 2013), drought (del Pozo and Ramirez-Parra, 2014; Martínez et al., 2018), drought and heat stress (Godfree et al., 2017), cold (Klatt et al., 2018), and light (Coate et al., 2013). The better stress response and adaptation of polyploids to abiotic conditions are probably under epigenetic control (del Pozo and Ramirez-Parra, 2014). Polyploidy changes the methylation profile under stressful environments, as reported, e.g. for Brassica napus after drought (Jiang et al., 2019).

Notably, stress conditions can also influence mode of reproduction, especially apomixis, the asexual reproduction via seed (Nogler, 1984a). Apomixis is widespread in angiosperms (Hojsgaard et al., 2014b), and occurs most frequently in polyploid cytotypes, but occasionally also in diploids (Grant, 1981; Carman, 1997; Hojsgaard and Hörandl, 2019). Gametophytic apomixis, the form of interest here, involves formation of an unreduced embryo sac from an unreduced megaspore via meiotic restitution of the megaspore mother cell (diplospory) or from a somatic cell of the nucellus tissue (apospory) (Asker and Jerling, 1992; Koltunow and Grossniklaus, 2003). Functional seed development through gametophytic apomixis involves three components: (1) apomeiosis (formation of unreduced embryo sac); (2) parthenogenesis (embryo development without fertilization of egg cell); and (3) functional endosperm development with male genome contributions from the pollen (pseudogamously) or independent from pollen (autonomously) (Nogler, 1984a). Male development is usually meiotic, but microsporogenesis is often disturbed, and hence final pollen quality is often strongly reduced (Asker and Jerling, 1992; Izmaiłow, 1996; Hörandl et al., 1997; Mráz et al., 2009). Apomixis is heritable (Ozias-Akins and van Dijk, 2007), and under genetic and epigenetic control (Grimanelli, 2012; Hand and Koltunow, 2014). Natural apomixis is frequently facultative, which means that the plant produces sexual and asexual seeds within one generation, often within the same flower or inflorescence (Bicknell et al., 2003; Aliyu et al., 2010; Cosendai and Hörandl, 2010; Hojsgaard et al., 2013; Schinkel et al., 2016).

Alternation of frequencies of asexual vs. sexual reproduction was observed under abiotic stress conditions, e.g. temperature, drought stress, salt stress, and photoperiod in many different genera (Evans and Knox, 1969; Saran and de Wet, 1976; Quarin, 1986; Gounaris et al., 1991; Klatt et al., 2016; Rodrigo et al., 2017; Klatt et al., 2018). Such a condition-dependent sex is also known from other asexual eukaryotes (Ram and Hadany, 2016). Abiotic stress leads to the accumulation of ROS (Reactive oxygen species) in plant tissues, which triggers oxidative damage, but also can initiate various epigenetic, genetic and hormonal signaling pathways for plant development (Halliwell, 2006; Foyer and Noctor, 2009; Huang et al., 2019). In the germline precursor cells, oxidative stress may increase the level of DNA double-strand breaks (DSBs) as initiator of meiosis. Here meiosis could act as DNA repair system (Hörandl and Hadacek, 2013). The above-mentioned studies on condition-dependent sex in plants support this hypothesis. In polyploids, however, an improved tolerance of stress conditions might decrease the stimulus for meiosis, and consequently trigger the alternative asexual development (Hörandl and Hadacek, 2013). However, a putative differential response of cytotypes to stress conditions with respect to mode of reproduction was so far not investigated.

We use here as a model system three cytotypes of the Ranunculus auricomus complex, a Eurasian polyploid complex with facultative, aposporous and pseudogamous apomixis (Nogler, 1984b; Hojsgaard et al., 2014a). In Central Europe, the R. auricomus complex comprises three closely related and genetically similar sexual progenitor species, and polyploid apomictic hybrids of these taxa (Hörandl et al., 2009; Hodač et al., 2014). One of the hexaploid hybrids (R. carpaticola x cassubicifolius) with facultative apomixis (Hojsgaard et al., 2014a) was used previously for testing the response to light stress. This previous experiment using extended photoperiod enhanced sexual megaspore formation in these hexaploid R. auricomus clones concomitant with oxidative stress (Klatt et al., 2016). In our study, we test the hypothesis that with the light stress treatment, diploids would respond more intensively to stress conditions with higher frequencies of sexual development than higher ploidy levels. Here we extend the treatment of (Klatt et al., 2016) to diploid, lower polyploid (tetraploid), and the same hexaploid plants to observe effects on mode of reproduction in different ploidy levels. To simulate the effect of extended photoperiod on the components of gametophytic apomixis, we study here two developmental steps, namely ovule formation, and seed formation. Since microsporogenesis is meiotic without an alternative asexual developmental pathway, we focus here on pollen quality as a possible factor for successful seed formation. The main aims of this research are to explore with light stress treatments whether ploidy level alters stress response with respect to mode of reproduction, and whether stress response correlates positively to sexual megaspore formation and/or proportions of sexual

### MATERIALS AND METHODS

### Plant Material

seed formation.

We used for the extended photoperiod experiment facultative apomictic plants of the Ranunculus auricomus complex from three different cytotypes. These cytotypes are hybrids that originated from three Central European parental species (R. cassubicifolius, R. carpaticola, and R. notabilis). The diploid plants were synthetic F2 hybrids of R. carpaticola x notabilis and represent sister or sibling individuals from two parental lines; see details of crossing design in Barke et al. (2018). We used these plants because natural diploid apomicts are not known for the R. auricomus complex. The tetraploids were garden offspring of Ranunculus variabilis, which is a putative natural allopolyploid of the R. carpaticola/cassubicifolius lineage and the R. notabilis lineage, and occurs sympatrically with the parental species in Central Europe (Hodač et al., 2014). The hexaploids were garden offspring of Ranunculus carpaticola x cassubicifolius, the same plants as used by Klatt et al. (2016). Hence, all cytotypes are hybrids, and they share the genetic background of closely related parental species (Hörandl et al., 2009). Since the parental taxa and the natural hybrids occur all in the same geographical area and altitudinal zone (Hörandl et al., 2009), we can also assume that they are all pre-adapted to the same natural light conditions. The ploidy level of tetraploids was ascertained using flow cytometry following methods of (Klatt et al., 2016). A list of materials with an identity number and ploidy levels is given in the Appendix (Supplementary Table 1). Plants were cultivated in the old botanical garden of the University of Goettingen from summer to winter for exposure to natural conditions, to stimulate the flower initiation.

### Growth Chamber Setup

The plants were moved into the climate growth chamber when sprouting at the beginning of the spring season.  The 1st year experiment was started from the first week of March 2017; the 2nd year was started from first February 2018. A total of c. 25 plants from each cytotype were grown with 10-h photoperiod (control) and 16 plus 0.5-h photoperiod (stress treatment) following (Klatt et al., 2016). Temperature setup and relative humidity were kept stable at 18°C and 60%, respectively. The light intensity was measured with a photometer (3415F Quantum Light Meter, Spectrum Technologies, Inc, Plainfield, USA) as photoactive radiation (PAR) c. 250 µmol m-2 s -2 (measured at shoot tips).

#### Plant Genotyping

Genotyping by Simple sequence repeats (SSRs) was applied to verify the plant's clonality and the relationships of cytotypes. We conducted SSRs only to tetraploid plants following methods by (Klatt et al., 2016). The SSR data for the other two cytotypes were derived from (Barke et al., 2018) for diploids and (Klatt et al., 2016) for hexaploids. Genomic DNA was performed by extracting dried leaf samples using Invisorb® Spin Plat Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. Multiplex Polymerase Chain reaction (PCR) was conducted at 25 µl volumes, containing 1 µl template DNA, 12.5 Roti®-Pol TaqS Master mix (Carl Roth GmbH + Co. KG, Karlsruhe, Germany), 1 µl Forward Primer, 1 µl Reverse Primer, 0.125 µl MgCl2, 1 µl CAG-Primer (FAM or HEX labeled). PCR reactions were run in a BIORAT™ Thermal Cycler. PRC machine setting was: 94°C for 10 min, then 14 x (denaturation at 94°C for 60 s, annealing at 62°C+ 0.5°C per cycle for 90 s, extension at 72°C for 60 s), followed subsequently by 35 x (denaturation at 94°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 30 s), last extension step at 72°C for 60 s and final storage conditions at 4°C. PCR samples were adjusted before 85 µl formamide (HiDi) was added. This mixture was run in an automatic capillarity sequencer Genetic Analyzer 3130 (Applied Biosystems, Forster City, CA, USA) using Gene Scan 500 Rox (Applied Biosystems) as size standard after a denaturing pretreatment for 3 min at 92°C. Scoring of the electropherograms was done using Genemarker V2.4.2 (SoftGenetics LLC, State College, PA, USA) and exported as a binary matrix presence/ absence of alleles to characterize multilocus genotypes. We applied Neighbour-joining analysis based on Jaccard similarity index in FAMD to test the SSR profiles (Schlüter and Harris, 2006). Branch support values were derived from the majority consensus tree from 1000 bootstrap replicates. FigTree v1.4.2 (Rambaut, 2007) visualized the result.

#### Female Development

Development of embryo sacs was already previously characterized within the R. auricomus complex on both apomictic and sexual species and is quite uniform (Nogler, 1984a; Hojsgaard et al., 2014a; Klatt et al., 2016; Barke et al., 2018): the megaspore mother cell differentiates near the micropyle and undergoes meiosis, resulting in a megaspore tetrad. In sexual development only the chalazal megaspore develops further, and produces after three mitotic divisions a typical 7-celled, 8-nucleate Polygonum type embryo sac (with three antipodals, a binucleate central cell, two synergids, and one egg cell). Apomictic development is characterized by enlargement of a somatic cell in the nucellus which emerges in parallel and aside the megaspore tetrad, and continues embryo sac development into an unreduced Polygonum type embryo sac, whereas all megaspores abort. Embryological analysis of the female development was made at the end of sporogenesis and the beginning of gametogenesis, following Hojsgaard et al., 2014a and Barke et al. (2018). R. variabilis, the only taxon that was analyzed for the first time here, did not show any deviations in timing or type of development. Flower buds were fixed at Formalin: acetic acid: ethanol: dH2O (2: 1: 10: 3.5) (FAA) for 48 h, and stored in 75% ethanol (Hojsgaard et al., 2014a). The flower bud was treated by dehydrating in four steps of 30 min incubation in 1 ml of 70%, 95%, and 100% (two times). Then the flower buds were treated by clearing method in five steps of 30 min in 300 µl of upgrading series of methyl salicylate diluted in ethanol [25%, 50%, 70%, 85%, and 100%; (Young et al., 1979)]. The perianth of selected flower buds was removed, ovaries were dissected and mounted in methyl salicylate on glass slides. Female sporogenesis and early stages of sexual or aposporous gametophyte development were analysed with differential interface contrast (DIC) in a light transmission microscope (Leica DM5500B with DFC 450 Camera, LAS V41 software, Leica Microsystems, Wetzlar, Germany). The determination of sexual and asexual ovules was made by the absence or presence of aposporous initial cells (AIC), respectively (van Baarlen et al., 2002).We excluded ovules with unclear structure and aborted ones. We only considered the data from a plant that had a minimum of five observable ovules. Additional data from (Klatt et al., 2016) were added to increase the N value for the hexaploid cytotype.

#### Seed Set

After we collected the sample for embryological analysis, the remaining flowers were then manually pollinated to increase fertilization rates. In fruiting stages, we bagged a minimum of five peduncles with collective fruits with porous plastic bags to avoid seed loss. We harvested the mature collective fruits and evaluated the proportion of well-developed seeds (seed-set percentage) among ploidies per flower on individual according to Hörandl (2008). Well-developed seeds were stored at room temperature and were used for reproductive pathway analysis.

### Reproductive Pathway of Seed Formation

The reproductive pathway was evaluated by single seed flow cytometric seed screening (ssFCSS) (Matzk et al., 2000). Two steel balls grounded a single seed (Ø 4 mm) in a 2 ml Eppendorf tube in a TissueLyzer II (Qiagen, Hilden, Germany; 30 Hz s-1 for 7 s). Nuclear isolation and staining were attained in two steps using Otto buffers (Otto, 1990). In the first step, nuclear isolation, 200 µl Otto I buffer (0.1 M citric acid monohydrate, 0.5% v/v Tween 20) was added and hand shacked with the ground material for 30 s. The solution was then filtered (30 µm mesh, Celltrics® Münster, Germany) into plastic tubes (3.5 ml, 55 mm x 12 mm, Sarstedt, Nümbrecht, Germany). In the second step, staining, 800 µl otto II buffer [0.4 M Na2HPO2, ddH2O and charged with 3 ng/ml 4',6-diamidinophenyl-indole (Sigma-Aldrich, Munich, Germany)] was added to the filtrate, and the solution was measured directly in Flow cytometer (CyFlow® Ploidy Analyser (Sysmex Partec GmbH, Görlitz, Germany) in the Blue fluorescence (UV LED, gain 365). Histograms were analyzed using CyView™ V.1.6 software (Partec GmbH). The coefficients of variation were less than 8%. The ploidy levels of embryo and endosperm were determined, and peak indices (PI) (mean peak value of the embryo compared to the mean peak of endosperm) were assessed (Supplementary Figure 5). For a Polygonum type embryo sac with two polar nuclei, the peak index for a sexual seed is c. 1.5, while for asexual seeds it can be 2.0, 2.5, or 3.0, depending on the contribution of pollen nuclei to endosperm formation. We observed the following developmental pathways: Sexual, pseudogamous apomixis, autonomous apomixis, and BIII-hybrids (Hojsgaard and Hörandl, 2019). BIII-hybrids arise from an unreduced embryo sac, whereby egg cell and polar nuclei were fertilized.

### Pollen Stainability

Pollen stainability was determined on a minimum of 500 pollen grains per plant from all cytotypes in both chambers by using 10% Lugol's iodine (I2KI) solution, following methods by (Schinkel et al., 2017). The stainability of starch content was used as an indicator of viable pollen under a light microscope (LEICA DM5500B with DFC 450 C camera, LAS V41 software, Leica Microsystems, Wetzlar, Germany) at 400x magnification. The viable pollen grains were indicated by black-stained color, but brownish, reddish, and translucent (empty) pollen was counted as non-viable.

### Statistical Analyses

All data were tested for their normality distribution by Kolmogorov-Smirnov and Shapiro-Wisk test and for their homogeneity of variance with the Levene test. Female development, seed set, reproduction pathway of seed formation, and pollen viability were determined per flower as a percentage and subsequently averaged per plants. The percentage of data were arcsine transformed before statistical analysis. We tested the influence of treatment on mean sexual ovules and seed set among ploidies with General Linear Model (GLM) univariate (Two-way ANOVA) for completely randomized factorial design model, and means were compared according to the least significant difference (LSD) test at 0.05 probability level (pvalue < 0.05). Tukey HSD was performed to the means of sexual ovules to determine the main factors. Nonparametric Kruskal-Wallis and Mann-Whitney U-test were applied to test the influence of treatment on sexual seed formation per ploidy. Boxplots were plotted with untransformed percentage values and show the 25th, and 75th percentile ranges as a box, and the median as a black line: circles are outliers; asterisks are extreme values. All statistical analyses were performed with IBM SPSS Statistic 25 (IBM Deutschland GmbH).

## RESULTS

#### Female Development

The ovule development of all three cytotypes of the R. auricomus complex showed the same pattern of a typical Polygonum type embryo sac (Supplementary Figures 1–4). We had observed 6,505 ovules (c. 18 ovules per flower bud) among cytotypes at megasporogenesis and early megagametogenesis. At this stage, sexual and asexual ovules can be discriminated (Supplementary Figure 4).  During the next step, three cells aborted, and only the chalazal cell remained as functional megaspore. At megagametogenesis stage, the functional megaspore enlarged with the presence of vacuoles and continued with three nuclear divisions, resulting in a total of eight nuclei. Development of sexual ovules was indicated by the absence of any aposporous initial cell (AIC) during megasporogenesis and early megagametogenesis. On the other hand, in asexual ovules, one or more AIC was observed directly near the megaspores at the chalazal pole or near to this area, but at a different optical layer (Figure 1).

#### Effects of Ploidy, Treatment, and Combined Effect of Ploidy/Treatment to the Proportion of Female Development

Extended photoperiod enhanced the proportion of sexual ovules in all three cytotypes of the R. auricomus complex. The mean proportion of sexual ovules significantly increased from control treatment to stress treatment (80.37 (mean) ± 19.38 (sd) % to 99.26 ± 1.26%; p-value < 0.001) in diploid, (57.90 ± 8.79% to 80.29 ± 10.67%; p-value < 0.001) in tetraploids, and 52.61 ± 26.11% to 70.36 ± 20.04%; p-value = 0.006) in hexaploids (Figure 2).ANOVA revealed significant alterations by the main effect photoperiod (pvalue < 0.001) and ploidy (p-value < 0.001), but not by the interrelationship between them (p-value = ns) (Table 1). Tukey HSD revealed significant differences in control treatment between diploids and hexaploids (p-value = 0.047) and in stress treatment between diploids and polyploids (p-value < 0.001) but there is neither a significant difference between tetraploids and hexaploids in the both treatments nor among diploids and tetraploid in the control treatment (p-value = ns) (Supplementary Table 3).

### Seed Set

Extended photoperiod did not influence the proportion of welldeveloped seeds among cytotypes of R. auricomus complex. Our investigation of 83 individuals revealed no significant difference in seed set between plants grown in control and stress chamber (p-value = ns) (Figure 3). Diploid plants under stress treatment produced a higher mean but not significant different proportions of well-developed seeds (mean value = 50.22%) compared to control treatments (mean value = 39.84%; p-value = 0.300). Tetraploid plants under stress treatment produced a mean of 28.97% compared to a mean of 31.09% (p-value = 0.459) under control treatment. Hexaploid plants under stress treatment produced a mean of 43.04% compared to a mean of 42.17% (pvalue = 0.880) under control treatment. A two-way ANOVA revealed only significant differences between the ploidies (p-value < 0.001), but neither a significant effect on treatment nor an interaction effect (p-value = ns) (Supplementary Table 4). Multiple comparison tests revealed that significant differences were observed between diploids and tetraploids (p-value < 0.001; Tukey HSD) and between tetraploids and hexaploids (p-value < 0.001; Tukey HSD) (Supplementary Table 5).

#### Reproductive Pathways of Seed Formation

Extended photoperiod did not enhance the proportion of sexual seeds over ploidies. The mean value of the proportion of sexual seeds was not significantly different between treatments among ploidies (p-value = ns, Mann-Whitney U-test) (Figure 4, Supplementary Table 6). Analysis of 1,468 seeds among ploidies indicated several reproductive pathways in the R. auricomus complex (Table 2). In diploid plants, the majority of seeds was formed sexually while in tetraploid and

FIGURE 2 | Proportions of sexual ovules in the R. auricomus complex plants grown in climatic chamber under prolonged photoperiod (stress) and shortened photoperiod (control). Mean values and statistical significance are given in figure. N = number of individuals. For the test statistics, see Supplementary Table 2.

TABLE 1 | P-values for the two way ANOVAs to determine the interaction effect of stress treatment and ploidy level on the proportion of sexual ovules.

R Squared = 0.574 (Adjusted R Squared = 0.551)

hexaploid plants, asexuality was the most frequent reproduction mode (Figure 4). In diploid sexual seeds, we observed the ratio of embryo to endosperm DNA content of 2C:3C, which is the indication of double fertilization between reduced egg cell with one sperm cell [1(m)+1(p)] and two polar nuclei with the other sperm cell [1(m)+1(m)+1(p)], producing a Peak Index (PI) of 1.5. A few apomictic seeds were observed (two with pseudogamous endosperm and one with autonomous endosperm) only in the stress treatment. The pseudogamous endosperm comes from the development of an unreduced embryo [2(m)] and fertilization of two polar nuclei with one or two reduced or unreduced sperm cells [2(m)+2(m)+1(p) or 2 (p)], with ratios of embryo to endosperm of 2C:5C (PI = 2.5) and 2C:6C (PI = 3.0). Autonomous endosperm develops from an unreduced embryo [2(m)] and unfertilized of two polar nuclei (2Cm+2Cm) with the ratio of embryo to endosperm of 2C:4C (PI = 2.0), which is caused by the absence of paternal genome in seed development.

Tetraploid and hexaploid plants displayed more variation on the mode of seed reproduction. Sexual reproduction mode was present in 39 (6.2%) tetraploid seeds and 36 (7.5%) hexaploid seeds. Pseudogamous endosperm was the most frequent mode of seed formation and appeared in 543 (86.3%) tetraploid seeds and 433 (90.7%) hexaploid seeds. Generally, this mode of reproduction produced a PI value of 3.0. The less frequent forms of pseudogamous endosperm with a PI = 2.5 and PI = 4.0 originated from the contribution of one reduced sperm nucleus or two unreduced sperm nuclei. Autonomous endosperms (PI = 2.0) were the most infrequent mode of seed formation, in a total of four seeds (0.55%) from tetraploids and nine seeds (1.93%) from hexaploids. Another type of reproduction mode, i.e. partial apomixis with an unreduced egg cell fertilized by reduced pollen (BIII-hybrid), was more frequent in tetraploid plants (45 seeds or 12.43%) compared with only one case in hexaploid plants (Table 2).

#### Pollen Stainability

Extended photoperiod did not alter the proportion of viable pollen between treatments. The assessment through 34,348 pollen grains from 67 plants revealed no significant differences in pollen viability between plants of the same cytotype grown in both treatments (p-value = ns; see Supplementary Figure 6). Hexaploids produced a higher mean proportion of viable pollen (mean value = 64.6% in control treatment and 60.7% in stress treatment) compared to diploids (49.9% in control treatment and 52.9% in stress treatment) and tetraploids (50.3% in control treatment and 52.4% in stress treatment). Multiple comparison tests among ploidies revealed that the only significant differences were observed between tetraploid and hexaploid plants (p-value < 0.001; Tukey HSD; Supplementary Table 7).

### DISCUSSION

Mode of reproduction in the facultative apomictic plant is influenced by abiotic stress, e.g. by light (Knox, 1967; Saran and de Wet, 1976; Quarin, 1986; Klatt et al., 2016). However, these studies compared stress and control treatments only within the same cytotype. Under the same conditions, the degree of facultative apomixis is usually related to ploidy level (Delgado et al., 2016; Kaushal et al., 2018). In this study, we presented for the first time developmental patterns among three cytotypes of the R. auricomus complex under stress and control conditions. We tested the hypotheses that prolonged photoperiod enhances only the first component of apomixis, i.e., apomeiotic embryo sac development, with the expectation of a buffer effect of stress in polyploids. The other two apomixis components, i.e. parthenogenesis and endosperm development, were not affected by different photoperiods.

#### Effects of Ploidy, Treatment, and Combined Effect of Ploidy/Treatment to the Proportion of Female Development

Prolonged photoperiod enhanced the proportion of sexual ovules, with a greater effect on diploids but lesser effect on tetraploids and hexaploids. Enhancement on the proportion of sexual ovules after the same type of light stress had been reported before only in the hexaploid cytotype (Klatt et al., 2016). The hexaploids also formed a comparable proportion of sexual ovules

FIGURE 3 | Proportions of well-developed seeds in the R. auricomus complex plants grown in climatic chambers under prolonged photoperiod (stress) and shortened photoperiod (control). Mean values and statistical significance are given in figure. N = number of individuals. For the test statistic, see Supplementary Table 2.

under garden conditions (Hojsgaard et al., 2014a). The three cytotypes of the R. auricomus complex exhibited a similar mode of reproduction as the pairwise comparison of data revealed insignificant differences between ploidies in control treatments. The result of controls and also the high genetic similarity of the three cytotypes (Supplementary Figure 7) make it unlikely that slightly different genetic backgrounds of the cytotypes had influenced the results of our experiments. The proportion of sexual ovules of the diploid cytotype grown in the garden, ranging from 45% to 82% (Barke et al., 2018), was still within the range of our data. These plants represent recently formed synthetic F2 hybrids (Barke et al., 2018) with lower proportions of apospory than in the polyploids that already had established apomixis in the natural source populations. However, despite these more lineage-specific features, differential effects of treatments were observed in all three cytotypes in the early stages of development.

The prolonged photoperiod (16 plus 0.5 h) may have expanded the accumulation of ROS (Reactive oxygen species) in the reproductive tissue, as reported for the hexaploids based on analysis of secondary metabolite profiles (Klatt et al., 2016). Results support the hypothesis that the oxidative lesions might mobilize the meiotic DNA repair system in the megaspore mother cell and trigger meiosis and megasporogenesis (Hörandl and Hadacek, 2013). This stimulus might increase the proportion of functional megaspores as a cellular survival strategy for the germline (Rodrigo et al., 2017), as shown remarkably in our diploids. Differential genetic stress regulation of sexual and apomictic plants was also observed in seedlings of Boechera, and may be important for the bypass of the meiotic pathway (Shah et al., 2016).

In tetraploids and hexaploids, the oxidative stress of prolonged photoperiods might be different. This could be due to altered photosynthetic electron transport capacities (Coate et al., 2013), or to altered secondary metabolite profiles in polyploids and hybrids (Orians, 2000). We speculate that lowered oxidative stress in polyploids might not be severe enough to induce sufficient double strand breaks that would be essential for a correct processing of meiosis (Keeney et al., 2014). Consequently, meiosis and megasporogenesis might be disturbed. Failure of megasporogenesis might release aposporous initial cell (AIC) development. Cell-specific transcriptome studies on aposporous Hieracium subg. Pilosella suggested that contact and cross-talk between AICs and functional megaspores could be the trigger for mitotic development of the former and degeneration of the latter (Juranic et al., 2018). We suppose a similar interaction of AICs and megaspores in the Ranunculus auricomus complex as they always occur together in close neighborhood, and we observed the presence of AICs together with young (2-nucleate stage) meiotic embryo sacs but not at later stages (Supplementary Figures 2C, D and Supplementary Figure 4). The emergence of aposporous initials starts in the Ranunculus auricomus complex mostly at the end of megasporogenesis and is correlated to disturbance of megasporogenesis. The surviving aposporous cells grow faster than the meiotic cell and occupy the mode of megagametogenesis and seed development (Hojsgaard et al., 2014a; Barke et al., 2018). The stress only affects the megaspore but leaves apomeiosis as the surrogate for the sexual pathway (Hörandl and Hadacek, 2013). Alternatively, polyploids with more DNA content have more repair templates for the DSBs, and a higher dose of stress would be required to break the DNA (Schoenfelder and Fox, 2015). Here polyploidy might promote DNA damage tolerance under elevated stress as described (Schoenfelder and Fox, 2015) and buffers stress effects (Hörandl and Hadacek, 2013).

TABLE 2 | Observed reproductive pathways of three cytotypes of the R. auricomus complex.

†Autonomous endosperm.

‡Pseudogamous endosperm, polar nuclei were fertilized by one reduced/unreduced or two reduced/unreduced sperm nuclei.

Cx reflects ploidy based on DNA content: m, maternal genome contribution; p, paternal genome contribution; PI, peak index.

Environmental stress plays a role as an inhibition factor under an epigenetic mechanism that disturbs or interrupts the silencing signal of apomictic-conditioning (Rodrigo et al., 2017). At least in diploid Ranunculus, the treatment might strengthen a signal transduction pathway that promotes switching from apomeiosis to meiosis, as demonstrated in facultative Boechera after drought stress (Mateo de Arias, 2015; Gao, 2018; Carman et al., 2019). In polyploids, the whole duplication genome (WDG) provides the co-loss or co-retention condition, which maintains a constant set of miRNA for basic biological functions (Liu and Sun, 2019). Our data suggested that polyploids respond to the stress via homeostatic regulation in the frequency of apospory vs. megasporogenesis. The high variability of the proportions of sexual ovules among our genetically identical polyploids supports the findings of epigenetic and transcriptional control mechanisms as the background for the phenotypic expression of apospory (Grimanelli, 2012; Schmidt et al., 2014). Our result supports the hypothesis that phenotypic features of apomixis in flowering plants are strongly affected by polyploidy (Delgado et al., 2016; Kaushal et al., 2018) and subjected to epigenetic control (Rodrigo et al., 2017).

#### Effects of Ploidy, Treatment, and Combined Effect of Ploidy/Treatment to the Seed Development and Mode of Reproduction

The prolonged photoperiod affected neither the frequency of seed set, the proportion of sexual seeds, nor the pollen viability. Ranunculus auricomus complex plants generally lose a high seed proportion compared to rates of ovule formation due to their high seed abortion rate, exceeding one-half to two-third (Izmaiłow, 1996; Hörandl, 2008; Hörandl and Temsch, 2009; Klatt et al., 2016; Barke et al., 2018). This failure on seed formation arises at early stages and during the development of endosperm tissue (Barke et al., 2018). The diploid cytotype, which generally reproduces sexually, delivers a better seed set than the higher ploidy levels. In contrast, tetraploids and hexaploids, which are predominantly facultative apomictic, showed a reversed pattern, by increasing frequencies of asexual seeds.

Pollen quality is an external factor influencing the seed set of all cytotypes. The great variation in pollen quality, as observed here, is typical for apomictic plants (Asker and Jerling, 1992). The lower quality of tetraploid pollen was concomitant with a lower seed set of the tetraploids, while the better pollen quality in diploids and hexaploids corresponded to a higher seed set in these cytotypes. For seed formation, the contribution of a male gamete to fertilize the central nuclei is the major requirement for proper endosperm development (Vinkenoog et al., 2003). The diploids keep their sexual ovules growing into sexual seeds in both treatments, while the survival of three apomictic seeds in the stress treatment represented rare exceptions from seed abortion. Similar results have been reported from the garden experiment (Barke et al., 2018). Diploid plants are sensitive to genomic imprinting deviation in the endosperm (Hörandl and Temsch, 2009; Barke et al., 2018), i.e. a 2:1 constant ratio for maternal (m) to paternal (p) genome contribution to endosperm (Spielman et al., 2003; Vinkenoog et al., 2003). The occurrences of genome imbalance in pseudogamously (4m:1p and 4m:2p) and autonomously formed seed (4m:0p) suggested that endosperm imbalance inhibited apomictic seed formation in our diploid cytotype.

On the other hand, in polyploids, the development of sexual ovules aborted to a large extent and was replaced by aposporous initials that completed megagametogenesis. Apomictic seed formation in polyploids is mainly influenced by the competitive capacity of the unreduced embryo sac formation rather than by the light regime during megagametogenesis and seed development (Hojsgaard et al., 2013; Klatt et al., 2016; Hodač et al., 2019). The surviving aposporous initials continue to develop into aposporous embryo sacs, and seeds are formed mostly via parthenogenesis and pseudogamous apomixis. This mode of reproduction is indicated by the parthenogenetic embryo (an unreduced egg cell develops without male gamete fusion) and pseudogamous endosperm (two unreduced polar nuclei fuses with one or two male gametes). Parthenogenesis appears mostly in our asexual polyploid seeds as a significant factor promoting unreduced gametophytes against reduced one and seed formation (Hojsgaard and Hörandl, 2019). A significant number of BIII-hybrids in tetraploids were formed through fertilization of unreduced egg cells as partial apomixis, as it was also occasionally observed in other FCSS studies (e.g. Schinkel et al., 2016; Barke et al., 2018; Klatt et al., 2018). This BIII-hybrid had probably an extremely long period of egg cell receptivity in this cytotype as assumed in diploid Ranunculus (Barke et al., 2018). Additionally, pollen-independent seed development via autonomous apomixis was also a rare event in polyploids. Asexual seed formation via pseudogamy is predominant in most apomictic plants (Mogie, 1992) as observed in our polyploids. The most common developmental pathway, however, used both sperm nuclei, or the unreduced sperm nucleus, for fertilization of polar nuclei, and hence restored the optimal 2m:1p ratio in the endosperm; these pathways result in a peak index of 3.0 in flow cytometric seed screening and represent the major proportions of apomictic seeds in both tetraploids (92%) and hexaploids (88%), see data in Table 2. Unbalanced genome contributions were also observed. Even though the diploids are quite sensitive to genomic imprinting, the polyploids in Ranunculus are more relaxed as expected (Grimanelli et al., 1997; Quarin, 1999). The current theory suggests that epigenetic mutation in polyploids creates relaxation on genomic imprinting during endosperm development (Kaushal et al., 2018). This could be the reason of higher seed set in hexaploid than in tetraploid cytotypes, similar to in hexaploid Potentilla puberula that had higher seed set than the tetraploids (Dobeš et al., 2018). These findings suggest the presence of a buffer effect on genomic imprinting in polyploids.

Our results suggest that the light regime only affects the proportion of sexual ovules, but the effect does not continue on the mode of seed formation. This finding supports the oxidative stress initiation hypothesis (Hörandl and Hadacek, 2013) that light stress affects only female meiosis, but has no relevance to further development. Polyploids express predominantly apospory, probably by improved mechanisms to buffer the abiotic stress, and are able to establish apomictic seed formation. These findings are in line with the general observation that apomixis mostly occurs in polyploid plants, despite the fact that the pathway can occur in diploids as well, albeit in much lower frequencies. Hence, stress resistance of polyploids may indirectly facilitate the establishment of apomixis, but is not necessarily essential for its expression, as proposed by Hojsgaard and Hörandl (2019).

#### CONCLUSIONS

Three cytotypes of facultative R. auricomus complex express the alternation of proportions of asexual ovules into more sexual ovules after prolonged photoperiod. We hypothesize that light stress increases ROS formation that triggers oxidative stress. The oxidative stress might stimulate the meiotic DNA repair system in the megaspore mother cell and suppresses mitotic division, resulting in sexual ovules. The effect of prolonged photoperiod on megasporogenesis was most pronounced in diploids; the lower effect of light stress in polyploids is probably as a consequence of higher stress resistance. In polyploids, high rates of seed abortion left a lower proportion of sexual seeds, whereas in diploids the sexual pathway is still predominant. Seed formation is not influenced by environmental stress conditions, but rather depending on proper endosperm formation. Our findings shed light on the predominance of apomixis occurrence in polyploid plants.

## DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

## AUTHOR CONTRIBUTIONS

FU and EH designed research. FU performed research, analyzed and interpreted data. CC contributed to FCSS and microsatellite analysis. FU wrote the manuscript with contributions of EH.

### FUNDING

This project was funded by The German Research Fund DFG (DFG Hörandl Ho 4395 4-1) to EH and by the Indonesia endowment fund for education, grant no. PRJ-2369/LPDP.3/ 2016 to FU.

### ACKNOWLEDGMENTS

Silvia Friedrichs for nursing the plants; Birthe Barke for help with data interpretation; referees for valuable comments on the manuscript.

### SUPPLEMENTARY MATERIAL

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

#### REFERENCES

Aliyu, O. M., Schranz, M. E., and Sharbel, T. F. (2010). Quantitative variation for apomictic reproduction in the genus Boechera (Brassicaceae). Am. J. Bot. 97 (10), 1719–1731. doi: 10.3732/ajb.1000188

Asker, S., and Jerling, L. (1992). Apomixis in plants (Boca Raton: CRC Press).

in the endosperm of diplosporous apomictic Tripsacum (Poaceae). Sexual. Plant Reprod. 10 (5), 279–282. doi: 10.1007/s004970050098

Grimanelli, D. (2012). Epigenetic regulation of reproductive development and the emergence of apomixis in angiosperms. Curr. Opin. Plant Biol. 15 (1), 57–62. doi: 10.1016/j.pbi.2011.10.002

triggers metabolic reprogramming in facultative apomictic Ranunculus auricomus. Front. Plant Sci. 7, 278. doi: 10.3389/fpls.2016.00278

lovegrass (Eragrostis curvula). PloS One 12 (4), e0175852. doi: 10.1371/ journal.pone.0175852

Conflict of Interest: 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 © 2020 Ulum, Costa Castro and Hörandl. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) 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.