Genetic analysis of anemone-type and single-type inflorescences in chrysanthemum using genotyping-by-sequencing

Flower shape is a key trait of ornamental and commercial importance in breeding programs for chrysanthemum (Chrysanthemum morifolium Ramat.). Understanding the genetic basis of the phenotypic variation seen in inflorescence-related traits will contribute to genetic improvement and to the development of new varieties. In this study, we investigated the genetic determinants of inflorescence traits using an F1 segregating population derived from a cross between two cultivars with different inflorescence types, ‘Puma White’ (anemone-shaped inflorescence) and ‘Dancer’ (single-type inflorescence). Genotyping-by-sequencing identified 26,847 single-nucleotide polymorphisms (SNPs) between 182 F1 progenies and their parents. A genome-wide association study highlighted 17 SNPs mapping to 15 GBS-tags as being significantly associated with three inflorescence traits: flower type, number of ray florets, and disk flower diameter. No single SNP was associated with flower diameter. These SNP-harboring sequences defined ten candidate genes associated with inflorescence traits. We explored the transcript levels for nine of these in flower buds, disk florets and ray florets using publicly available genome and transcriptome data. These results will provide the genetic and genomic foundation to harness important horticultural traits and explore new avenues in chrysanthemum breeding.


Introduction
Chrysanthemum (Chrysanthemum morifolium Ramat.) belongs to the Asteraceae family and is an economically important ornamental plant used as cut flowers or potted plants, comprising a considerable proportion of the global flower industry . Phenotypic diversity in chrysanthemum is illustrated by its multiple flower types with various colors and sizes. Chrysanthemum blooms, or flower heads, are formed by many small disk florets (flowers) surrounded by ray florets (Yoshioka et al. 2010). The flower heads of chrysanthemum can take on various forms described as anemone, Abstract Flower shape is a key trait of ornamental and commercial importance in breeding programs for chrysanthemum (Chrysanthemum morifolium Ramat.). Understanding the genetic basis of the phenotypic variation seen in inflorescence-related traits will contribute to genetic improvement and to the development of new varieties. In this study, we investigated the genetic determinants of inflorescence traits using an F 1 segregating population derived from a cross between two cultivars with different inflorescence types, 'Puma White' (anemone-shaped inflorescence) and 'Dancer' (single-type inflorescence). Genotyping-by-sequencing identified 26,847 single-nucleotide polymorphisms (SNPs) between 182 F 1 progenies and their parents. A genome-wide association study highlighted 17 SNPs mapping to Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10681-022-03124-7. pompon, single, or double inflorescence types, among others (Chen et al. 2009). Anemone-type flower heads are characterized by a prominent flower center composed of elongated and colored disk florets, making them both aesthetically and commercially attractive (Anderson 2006). Single-type flower heads consist of a central and compact disk composed of hermaphrodite florets and are considered the primitive form (Dai and Chen 1997;Dai 2004;Chen et al. 2009).
Crossbreeding has contributed to the development of modern chrysanthemum cultivars (Yang et al. 2019). Cultivated chrysanthemum show great diversity in their morphology and characteristics such as inbreeding depression, self-incompatibility, and high levels of heterozygosity. Accordingly, crosses between two parents with contrasting phenotypic traits of interest will produce F 1 progeny exhibiting extensive phenotypic variation (Yang et al. 2019). As phenotypes are affected by the underlying genotype, the surrounding environment, and genotype × environment interactions, phenotype-based selection is difficult and ineffective. These limitations are currently circumvented with the help of marker-assisted selection (MAS), whereby the genes controlling target traits, or closely linked loci, are used to select lines based solely on their genotype rather than their phenotype. Using MAS, a large collection of progeny derived from crossbreeding can be screened quickly and efficiently irrespective of any effect from the environment. This allows us to explore the influence of a target trait during earlier growth stages than those typically associated with the given trait, for example those influencing flowering time. The identification of the necessary molecular markers has been fueled by genetic analyses such as quantitative trait locus (QTL) mapping and genome-wide association study (GWAS). In chrysanthemum, many studies attempted to determine the genes associated with important traits using traditional molecular approaches such as amplified fragment length polymorphisms (AFLPs), random amplified polymorphic DNA markers (RAPD), inter simple sequence repeats (ISSRs), restriction fragment length polymorphism markers (RFLPs), expressed sequence tag-simple sequence repeats (EST-SSRs), and sequence-related amplified polymorphisms (SRAPs) (Zhang et al. , 2011Fan et al. 2020;Yang et al. 2020). However, the number of available markers for genetic analysis is limited and the resulting genomic regions were large. In addition, genetic work in chrysanthemum is challenging due to its complex genome.
Cultivated chrysanthemums have a hexaploid genome whose stable chromosome number is most frequently (2n = 6x = 54), although some cultivars have been described as aneuploid with a chromosome number ranging from 47 to 67 (Dowrick 1953). It is still unclear whether this event involved autopolyploid or allopolyploid. A recent cytological and molecular analysis revealed segmental allohexaploid in cultivated chrysanthemums (Klie et al. 2014). Our previous study, based on the Ks value distributions, whole-genome triplication (WGT) event (Ks = 0.1), as well as spices divergence (Ks = 0.05). At Ks = 0.05 divergence time, while C. morifolium maintained the original polyploidy status, it is likely that this particular diploid species experienced a rapid diploidization event, including a reduction in chromosome number (Won et al. 2017). Although its gigantic genome size of 6-7 Gb (Van Geest 2017) and high heterozygosity in cultivated chrysanthemum genome status caused in ambiguous of auto or allopolydization evidence for a recent WGD or WGT.
The advent of next-generation sequencing (NGS) technologies has made possible and accelerated the pace of genetic analyses and the development of molecular markers (Varshney et al. 2009). NGS platforms now allow the global identification of millions of single-nucleotide poly-morphisms (SNPs) and insertion-deletion polymorphisms over the genome, which have greatly increased the analytic power of many genetic studies (Varshney et al. 2009;Huang et al. 2010;Phan and Sim 2017). The GBS methodology has been implemented in many plants to explore their genetic diversity, construct linkage maps, map QTLs, and perform GWAS (Nguyen and Lim 2019). Besides simpler diploid species, GBS has also been widely adopted in polyploid species such as wheat (Triticum aestivum) and sugarcane (Saccharum officinarum) (Poland et al. 2012;Balsalobre et al. 2017). A recent example applied GBS to an exploration of genetic diversity and phylogenetic analysis for wild chrysanthemum species using 7,758 SNP markers (Nguyen et al. 2020).
In this study, we used an F 1 population derived from a cross between two cultivars with different inflorescence types, 'Puma White' (anemone-type) and 'Dancer' (single-type) to dissect the genetic basis of inflorescence traits in chrysanthemum. The progeny, consisting of 182 lines, was subjected to phenotypic characterization, genome sequencing by GBS, and GWAS. We then used publicly available transcriptome data to examine the possible function of the candidate genes derived from GWAS on shaping inflorescence traits. The result of this study will provide much needed information to understand the genetic basis underlying ornamental traits and to enhance chrysanthemum breeding efficiency.

Plant materials
A population of 182 F 1 progeny was derived from a cross between two hexaploid chrysanthemum (C. morifolium) cultivars, 'Puma White' and 'Dancer' (Park et al. 2015). To minimize selfing, 'Puma White' was selected as the maternal parent because its anemone-type inflorescence produces very little pollen. All plant materials are currently maintained at the National Institute of Horticultural and Herbal Science, Rural Development Administration, Republic of Korea.

Phenotypic evaluation and data exploration
Fresh healthy stems from the original 182 F 1 population and the parents, were grown in the 128-plug tray for 3 weeks. Rooted cuttings were then transplanted and grown in a polyethylene greenhouse in September 2018 under natural photoperiod to induce flowering in short-day conditions. Four inflorescencerelated traits were investigated in the F 1 population and their parents at the flowering stage: flower type (FT), flower diameter (FD, in cm), number of ray florets (NRF), and disk flower diameter (DFD, in cm). For FT, inflorescences with cushion-like disk florets were considered to belong to the anemone type, while those with short and dense disk florets were defined as single flower heads. Traits were measured using randomly selected inflorescences for each F 1 line as three different F 1 with replicates. Exploration of the phenotypic data was performed and visualized in Microsoft Excel 2016 or with R software (R Core Team 2014). The difference in traits between the two parents was compared using a one-tailed Student's t-test. The normality of the traits in the F 1 population was examined by calculating the skewness and kurtosis of the relevant distributions and by conducting a Shapiro-Wilk normality test. The correlation between traits in the F 1 population was examined by calculating Kendall rank correlation coefficients. The difference in traits between the two inflorescence types of the F 1 population was examined by Wilcoxon rank sum test with continuity correction.

GBS library preparation and sequencing
Total genomic DNA (gDNA) was extracted from young leaves using the DNeasy® Plant Mini Kit (Qiagen) following the manufacturer's protocol. The quality and concentration of gDNA were assessed by electrophoresis on 1.2% (w/v) agarose gels, on a Dropsense96 spectrophotometer (Trinean), and with the Quant-iT PicoGreen dsDNA Assay kit (Invitrogen). gDNA from each sample was digested with the restriction enzyme ApeKI (New England Biolabs) and then used as template for the construction of a GBS library, as previously described (Elshire et al. 2011). After checking the quality of libraries on a 2100 Bioanalyzer instrument (Agilent), the libraries were pooled and sequenced using a NextSeq500 platform (Illumina) as 150 bp reads in single-end mode.
Sequencing data analysis and genotyping GBS data were analyzed using the pipeline Stacks v. 2.5 with default settings (Catchen et al. 2011(Catchen et al. , 2013. Raw reads were demultiplexed, filtered, and trimmed to 100 bp using process_rad-tags (cleaned raw read data) in Stacks and Cutadapt (Martin 2011). Stacks was run in de novo stack formation mode (Catchen et al. 2011). These assembled GBS-tags were replaced as a reference genome and named sequence tag. To track a set of loci or genotypes, the identified variants were finally filtered using VCFtools v. 0.1.15 (Danecek et al. 2011) with the following parameters: minor allele frequency ≥ 5%, mean read depth ≥ 5 reads, and missing genotype < 5%.

GWAS of inflorescence traits
With the set of filtered SNPs defined above, GWAS was performed using TASSEL 5.0 software (Bradbury et al. 2007) with the general linear model (GLM) algorithm. Quantile-quantile (QQ) plots 177 Page 4 of 12 Vol:. (1234567890) indicated that the naïve model consisting of a GLM without any control for population structure (Q) or kinship matrix (K) was a suitable model for analysis. The significance threshold for a trait and marker association was set to P ≤ 0.001. The P values were adjusted by the threshold of the Bonferroni correction for multiple tests (1/n), where n is the total number of SNPs used in the association analysis (Kan et al. 2015;Sun et al. 2017). Significant SNPs located in the same GBS sequence tag (GBS-tag) were considered as belonging to the same locus.
Identification and expression analysis of candidate genes GBS-tag sequences carrying trait-associated SNPs were searched against the genome of diploid Chrysanthemum nankingense (Song et al. 2018) using Basic Local Alignment Search Tool for nucleotide sequences (BLASTN) (Altschul et al. 1997) with an E-value cutoff of 1 × 10 −10 . To determine the transcript levels of candidate genes, transcriptome deep sequencing (RNA-seq) data from C. nankingense were downloaded for six tissues (roots, stems, leaves, flower buds, disk florets, and ray florets) (http:// www. amway abrc. com/ downl oad. htm) (Song et al. 2018). Reads were aligned to the C. nankingense genome and transcript levels were calculated and normalized as fragments per kilobase of transcript per million mapped reads (FPKM) using TopHat and Cufflinks (Trapnell et al. 2009(Trapnell et al. , 2012. Genes with FPKM ≥ 1.0 were considered to be transcriptionally expressed. Transcript levels for each gene were visualized as a heatmap using Z-score normalization in R software (R Core Team 2014). Additionally, candidate proteins were subjected to BLASTP against the Arabidopsis (Arabidopsis thaliana) protein sequences (from the TAIR10 annotation) at the Ensembl Plants website (http:// plants. ensem bl. org).

Variation in inflorescence traits
We characterized the extent of phenotypic variation in inflorescence traits in the 182 F 1 progeny derived from the two parental cultivars (Figs. 1 and 2). The cultivar 'Puma White' has anemone-type inflorescence that are composed of two layers of ray florets at the edge and cushion-like disk florets at the wide flower core (Fig. 1a), while the cultivar 'Dancer' has single inflorescence with one layer of ray florets at the edge and flattened disk florets in a narrow arrangement at the center of the flower (Fig. 1b). We classified their F 1 progeny into 60 anemonetype and 122 single-type inflorescences ( Fig. 1c and Online Resource 1A). FD was the same between the two parents, while NRF and DFD were significantly greater in 'Puma White' than in 'Dancer' (P < 0.001) ( Table 1). In the F 1 progeny, FD ranged from 2.83 to 6.83 cm, with a mean of 4.55 cm; NRF varied from 12 to 66 ray florets, with a mean of 29.24; and, DFD ranged from 1.00 to 3.00 cm, with a mean of 1.73 cm. All three phenotypic traits therefore showed strong transgressive segregation, as the range of phenotypic values extended well beyond that defined by the parental lines (Fig. 2). We next calculated the coefficient of variation for each trait: 16.04% for FD, 26.00% for NRF, and 20.80% for DFD. Shapiro-Wilk normality tests revealed that all three traits significantly deviated from normal distributions (P < 0.05). We observed significant correlations (p < 0.05) between FD and NRF, FD and DFD for the F 1 population (Online Resource 2). A comparison of traits between anemone-type and single-type inflorescence of the F 1 population indicated that DFD was significantly larger in anemone-type inflorescence than in single-type inflorescence (P < 0.001), while FD and NRF were not different between the two types of inflorescences (Online Resource 1B).

Identification of SNPs by GBS
As a prerequisite for genomic analysis, we genotyped all individuals in the F 1 population and their parents by GBS, resulting a total of 688,806,674 raw reads, with an average of 3,743,515 reads per individual (Table 2 and Online Resource 3). We then filtered and preprocessed these raw reads to yield 326,317,886 high-quality reads, or 47.40% of the initial number of raw reads. We then implemented the Stacks pipeline to assemble all clean reads into sequence tags covering 56,486,574 bp in length, to which we mapped the clean reads to identify SNPs. The number of mapped reads per individual ranged from 453,348 to 1,693,058, with an average of 861,374. The Stacks pipeline discovered a total of 240,117 SNPs, which Page 5 of 12 177 Vol.: (0123456789) we narrowed down to 26,849 SNPs for subsequent analysis based on the selection criteria of minor allele frequency of at least 5%, read support from at least five reads, and missing genotype below 5% across the F 1 population.

GWAS and candidate gene prediction
To identify the loci associated with variation in inflorescence traits in chrysanthemum, we performed GWAS using the genotype (26,849 SNPs) and phenotype (FT, FD, NRF, DFD) data from all 182 F 1 population. The association analysis yielded 17 significantly associated SNPs mapping to 15 distinct GBS-tags from the de novo assembly (P ≤ 1.1 × 10 −06 ) (Online Resources 4 and 5). Of these 17 SNPs, 11 mapped to 10 GBS-tags for FT, 3 mapped to 3 GBStags for NRF, and 3 SNPs mapped to 2 GBS-tags for DFD. We detected no significant association between genotype and FD. The GBS-tags GBS-tag31024 and GBS-tag20699 each harbored two SNPs associated with FT and DFD, respectively.
To identify the candidate genes underlying each trait, we queried the genome of diploid C. nankingense, one of the progenitors of C. morifolium (Song et al. 2018) with the 15 GBS-tags of significant SNPs. BLAST searches determined that 12 GBS-tags showed significant hits in the C. nankingense genome (E-value ≤ 2.0 × 10 −26 , percent identity ≥ 92%, and query coverage ≥ 84%), while 3 GBS-tags did not match any sequence in C. nankingense or across the National Center for Biotechnological Information (NCBI) non-redundant (nr) database (Online Resource 6). Of the 12 GBS-tags with a match in the C. nankingense genome, 2 (harboring 3 SNPs) mapped to intergenic regions, but the remaining 10 GBS-tags (and their 11 SNPs) located to genes. Among these 11 genic SNPs, 1 was in a 5' untranslated region (UTR) and 4 were within introns. Two other SNPs introduced a synonymous mutation that did not affect the protein sequences of their associated genes CHR00045339 and CHR00036293, encoding a homolog of ubiquitin-specific protease 16 (UBP16) and UDP-glucosyltransferase 85A1 (UGT85A1), respectively (Table 3). The remaining four SNPs introduced missense mutations in three genes and affected the encoded protein. The SNP in GBS-tag1536:83 mapped to the gene CHR00032674 and resulted in a K46Q amino acid substitution in the encoded aspartic peptidase (AP) anchored to the plasma membrane. The SNP within GBS-tag2771:23 affected gene CHR00074055 by introducing a N79K amino acid substitution in the encoded lipase class 3 family protein. Finally, two SNPs mapping to the same gene (CHR00047948) introduced T131A and G149S amino acid substitutions in the encoded octicosapeptide/Phox/Bem1p protein homolog. arrowheads. Note that flower diameter is the same between the two parents. a Flower diameter. b Number of ray florets. c Disc flower diameter Table 1 Descriptive statistics associated with the variation for flower-related traits based on a set of 182 F 1 population derived from the cross between 'Puma White' and 'Dancer' a Flower diameter (FD, cm); number of ray florets (NRF), disc flower diameter (DFD, cm) b Standard deviation *** Significant between two parents at the 0.001 level according to one-tailed t-test We next examined transcript levels for the 10 candidate genes harboring SNPs using public RNA-seq data from C. nankingense (Song et al. 2018). With the exception of CHR00074055, nine of our candidate  genes appeared to be expressed (expression cutoff of FPKM ≥ 1) (Online Resource 7). The expression of candidate genes involving significant SNP was detected in flower buds, disk florets, and ray florets; each candidate gene also demonstrated differential expression between analyzed tissues (Fig. 3).

Discussion
Molecular genetic analysis and MAS are limited and slow to be implemented in cultivated chrysanthemum (C. morifolium Ramat.) due to its large, complex and hexaploid genome, self-incompatibility, and high level of heterozygosity. Previous genetic work revealed that inflorescence traits are regulated by multiple loci rather than a single major locus in chrysanthemum (Tang et al. 2015). However, several genes were later reported to be involved in shaping inflorescence traits (Chong et al. 2016(Chong et al. , 2019. In this study, we investigated the genetic architecture of inflorescence traits using a set of 182 F 1 segregating progeny derived from a cross between the anemonetype cultivar 'Puma White' and the single-type cultivar 'Dancer'. We used GBS to identify SNPs across the genome in all progeny and their parents, followed by GWAS to detect SNPs significantly associated with each inflorescence trait. The resulting candidate genes will provide valuable resources for chrysanthemum breeding. Genetic analysis using segregating populations is a powerful tool to identify loci responsible for phenotypic variation in target traits. The F 1 progenies used in this study was obtained from a cross between two chrysanthemum cultivars with different inflorescence types and showed substantial variation and transgressive segregation in the flower traits examined, making this F 1 population a suitable mapping population. In the F 1 progeny, more individuals displayed singletype inflorescence than the anemone-type, suggesting that the single inflorescence morphology may be a dominant phenotype in this particular crossing combination, which should be examined further in more diverse germplasms. Progeny with anemone-type flower heads had larger DFD than did single-type progeny, which was consistent with the differences between the parents. By contrast, FD and NRF were similar in the two groups of F 1 progeny with different inflorescence types. We observed a weak significant correlation between the phenotypic traits. A relatively negative correlation was observed between NRF and DFD, which is consistent with the result in Fan et al. (2020). For example, the correlation of − 0.055 was in our study and the value was − 0.45 in Fan et al. (2020). A correlation between FD and DFD was observed within a panel of anemone-type chrysanthemum accessions (Yang et al. 2020). These findings suggest that flower traits are complex in chrysanthemum.
Traditional molecular markers have enabled genetic studies in chrysanthemum such as investigation of inheritance, linkage map construction, QTL mapping, GWAS, and phylogenetic studies, but the number of such markers is insufficient (Zhang et al. , 2011Chong et al. 2019), only reaching hundreds of markers in biparental populations and natural populations of chrysanthemum (Zhang et al. , 2011Li et al. 2016;Fan et al. 2020). With NGS technologies, however, far more SNP markers can be discovered quickly and at low cost even in species lacking a reference genome, which is expected to substantially improve the development of new molecular markers and breeding in chrysanthemum (Su et al. 2019a). GBS identified 7,758 SNPs across 80 wild chrysanthemums accessions (Nguyen et al. 2020) and 26,849 SNPs from 182 F 1 progeny and their parents in the present study. Another NGS method, double digest restriction site-associated DNA (dd-RAD-seq) identified 9,219 SNPs in an F 1 biparental population to dissect the genetic control of petal color in ray florets (Sumitomo et al. 2019). Lastly, specific locusamplified fragment sequencing (SLAF-seq) detected approximately 92,000 SNPs, which helped identify alleles associated with inflorescence traits and waterlogging tolerance in cultivated chrysanthemum accessions (Chong et al. 2016(Chong et al. , 2019Su et al. 2019b).
Inflorescence traits are some of the most important factors that determine the ornamental and commercial value of chrysanthemum. GWAS identified 17 SNPs significantly associated with inflorescencerelated traits mapping to 15 GBS-tags. GWAS is a powerful tool to identify loci associated with specific traits (Yu et al. 2006). In a previous study, GWAS identified DNA markers associated with petal color in a segregating F 1 mapping population of chrysanthemum, in which inheritance is typically rather complex due to its autohexaploid genome (Sumitomo et al. 2019). The random pairing of chromosomes in autopolyploids complicates the calculation of recombination frequencies during QTL analysis, although it is not impossible. Nevertheless, our result demonstrated the applicability of the allohexaploid genome by genetic and gene functional analysis in C. morifolium. Three GBS-tags (GBS-tag50705:95, GBS-tag1280:27, GBS-tag12113:80), specific to C. morifolium were newly discovered based on the C. nankingense genome. A similar approach to our study has been confirmed in other plants as well. In the allotetraploid rapeseed (Brassica napus), GWAS constituted a robust approach to identify candidate genes involved in the regulation of fatty acid composition in a biparental doubled-haploid mapping population (Gacek et al. 2017).
Twelve of the 15 GBS-tags harboring the 17 significant SNPs aligned well to the C. nankingense genome, supporting the notion that the two species are closely related. We hypothesize that the three remaining GBS-tags are specific to C. morifolium. We assigned these 12 GBS-tags to 10 candidate genes potentially associated with inflorescence traits. Publicly available RNA-seq data indicated that all but one of our candidate genes are expressed in flower buds, disk florets, and/or ray florets. Among them, two other significant SNPs affected a gene encoding an octicosapeptide/Phox/Bem1p family protein.
The homologous gene in Arabidopsis was shown to be rapidly upregulated in roots within 30 min of treatment with the phytohormone auxin; its knockout mutant exhibited reduced lateral root growth in the presence of auxin (Pu et al. 2019). One member belonging to this protein family was identified during a GWAS using a collection of barrelclover (Medicago sativa) accessions for its association with a plant branch trait, and its overexpression resulted in an increased number of lateral branches in Arabidopsis (Wang et al. 2020). The octicosapep-tide/Phox/ Bem1p protein family contains Phox/Bem1 (PB1) domains that mediate protein-protein interactions through different charges at their protein surfaces (Korasick et al. 2015). PB1 domains are also present in two classes of auxin signaling regulators, auxin response factors and auxin/indol 3-acetic acid inducible proteins and mediate their interaction (Korasick et al. 2015). We therefore hypothesize that this PB1-containing protein probably regulates inflorescence development and growth in chrysanthemum by facilitating interactions between proteins in the auxin signaling pathway. For the remaining genes containing SNPs in their 5' UTR or introns, or with synonymous mutations (UBP16 and UGT85A1), alternative splicing or differential expression might be involved in the regulation of inflorescence traits. Alternatively, 177 Page 10 of 12 Vol:. (1234567890) tightly linked nearby genes might be responsible for these traits.
This GBS-based association study helped the identification of potential SNPs and genes with phenotypic effects on important horticultural traits in the present study. Of them, ten candidate genes contribute to the inflorescence traits were identified in this study, and awaits the validation and direct evidence of the roles of the candidate genes to support chrysanthemum breeding. The findings of the present study provide the genetic and genomic foundation to harness important horticultural traits and explore new avenues in chrysanthemum breeding.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.