Felipe Pérez De Los Cobos

Eva Coindre

Naima Dlalah

Bénédicte Quilot-Turion

Ignasi Batlle

Pere Arús

Henri Duval
email: henri.duval@inrae.fr

Almond population genomics and non-additive GWAS reveal new insights into almond dissemination history and candidate genes for nut traits and blooming time

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

One of the landmarks of human history was the transition from nomadic hunter-gatherer societies to settled agriculture-based societies. This transition, known as the Neolithic Revolution, marked the beginning of the domestication of wild plant species as cultivated crops [START_REF] Weisdorf | From foraging to farming: explaining the Neolithic revolution[END_REF]. Domestication, and later dispersal and diversification of crops, introduced substantial changes into their genomes, fixing alleles of interest, creating different genetic populations, and adapting these groups to different environmental conditions [START_REF] Gross | Genetic perspectives on crop domestication[END_REF][START_REF] Doebley | The molecular genetics of crop domestication[END_REF]. These changes, accumulated over thousands of years, led to the crops we consume today.

Nowadays, the main actor changing crop genomes is modern breeding. The efficient implementation of breeding strategies requires the correct management of germplasms, optimized genotyping and phenotyping methods, concise knowledge of crop genetic structure, and the study of genetic determinism behind traits of interest [START_REF] Thudi | Genomic resources in plant breeding for sustainable agriculture[END_REF][START_REF] Swarup | Genetic diversity is indispensable for plant breeding to improve crops[END_REF]. Genome-wide association study (GWAS) is a powerful tool to study quantitative traits in plant breeding. It aims to find polymorphic genetic markers (typically SNPs) significantly associated with phenotypic variation [START_REF] Korte | The advantages and limitations of trait analysis with GWAS: a review[END_REF]. However, one of the weaknesses of this technique is that most GWAS models assume that the genotypic variation has an additive effect on the phenotype. This means that non-additive effects, such as dominantrecessive or overdominant interactions, are not included in the models even when they may be relevant for most traits [START_REF] Tsepilov | Nonadditive effects of genes in human metabolomics[END_REF].

Almond [Prunus dulcis Miller (D.A. Webb)] is the most economically important temperate nut tree worldwide. In the period 2011 to 2021, its production increased 54%, reaching 1 684 395 metric tons of kernel [START_REF] Council | Nuts & dried fruit statistical yearbook[END_REF]. It belongs to the Rosaceae family and the Prunus genus with other important crops, including peach, plum, apricot, and cherry.

While research on the almond domestication process is in the early stages, some important insights have been made. The most accepted theory to date is that almond originated from hybridizations with several wild relatives somewhere between the Eastern Mediterranean and Southwest Asia, expanding rapidly to Central Asia and the Western Mediterranean. Many studies using different approaches support this theory, from analyses based on morphology, habitat, and/or coexistence in cultivated areas [START_REF] Ladizinsky | On the origin of almond[END_REF][START_REF] Grassely | Origine et évolution de l'amandier cultivé[END_REF][START_REF] Denisov | Almond genetic resources in the USSR and their use in production and breeding[END_REF], through genomic analyses [START_REF] Zeinalabedini | The origin and dissemination of the cultivated almond as determined by nuclear and chloroplast SSR marker analysis[END_REF][START_REF] Delplancke | Evolutionary history of almond tree domestication in the Mediterranean basin[END_REF][START_REF] Delplancke | Combining conservative and variable markers to infer the evolutionary history of Prunus subgen. Amygdalus s.l. under domestication[END_REF], to archaeobotanic evidence [START_REF] Pérez-Jordà | The emergence of arboriculture in the 1st millennium BC along the Mediterranean's "far west[END_REF][START_REF] Willcox | Early Holocene cultivation before domestication in northern Syria[END_REF][START_REF] Zohary | Domestication of Plants in the Old World. 4th Editio ed[END_REF]. In this sense, many efforts have focused on analyzing the population structure of different almond germplasms [START_REF] Font | Identification of genetic loci associated with quality traits in almond via association mapping[END_REF][START_REF] Guardo | Genetic characterization of an almond germplasm collection and volatilome profiling of raw and roasted kernels[END_REF][START_REF] Pavan | Almond diversity and homozygosity define structure, kinship, inbreeding, and linkage disequilibrium in cultivated germplasm, and reveal genomic associations with nut and seed weight[END_REF]. Nevertheless, these studies have been limited by the geographical origin of the accessions (most accessions came from the same region) or by the low number of markers used. As a result of these drawbacks, our knowledge of the genetic structure of the cultivated almond is still limited. The recent publication of three almond reference genomes [START_REF] Sánchez-Pérez | Mutation of a bHLH transcription factor allowed almond domestication[END_REF][START_REF] Alioto | Transposons played a major role in the diversification between the closely related almond and peach genomes: results from the almond genome sequence[END_REF][START_REF] D'amico-Willman | Whole-genome sequence and methylome profiling of the almond [Prunus dulcis (mill.) D.a. Webb] cultivar 'nonpareil'. G3 Genes, Genomes[END_REF] and the development of a 60 K SNP array [START_REF] Duval | Development and evaluation of an AxiomTM 60K SNP Array for almond (Prunus dulcis)[END_REF] have opened the door to performing genome-wide analyses on almond. So far, there have been three GWASs using genome-wide marker data [START_REF] Guardo | Genetic characterization of an almond germplasm collection and volatilome profiling of raw and roasted kernels[END_REF][START_REF] Pavan | Almond diversity and homozygosity define structure, kinship, inbreeding, and linkage disequilibrium in cultivated germplasm, and reveal genomic associations with nut and seed weight[END_REF][START_REF] Sideli | Genome-wide association analysis and validation with KASP markers for nut and shell traits in almond (Prunus dulcis [mill[END_REF], and several QTLs linked to shell and kernel quality traits were identified. Nevertheless, these studies only focused on additive GWAS models and the variability of the plant material was reduced. Studying a broader germplasm and analyzing nonadditive genotype-phenotype associations would allow the exploration of the origin and historic dissemination of the cultivated almond at the same time that would help to find alleles of interest and QTLs fixed over thousands of years of domestication.

In this study, we explored the genetic structure and genotypephenotype associations in a collection of 243 almond accessions from different origins. For this purpose, we first characterized the genetic diversity of the collection using the almond 60 K SNP array. Then we carried out a GWAS using additive and non-additive models for different traits, including kernel and nut weight, crackout percentage, double kernels percentage, and blooming time. As far as we know, this is the first non-additive GWAS in Rosaceae species. Using candidate gene analysis, we also proposed candidate genes responsible for the main QTLs found in this analysis.

Results

Almond genetic structure defined five ancestral groups

The results from the population structure analysis, additive kinship, phylogenetic tree, and PCA indicated that the most accept-able genetic structure of the 152 landraces included in this analysis is a model with five ancestral populations (Figure 1). According to the population structure analysis (Figure 1D), 80 accessions were part of one of the five ancestral groups, while 72 were considered admixed. The number of accessions was homogeneous between groups, ranging from 12 to 21 accessions (G4 and G2 respectively). Ancestral group G1 was mainly composed by Asian accessions: 12 Chinese, two Iranian, two Turkish, and one Pakistani. This group also included one Greek and one Romanian accession. Ancestral group G2 was entirely formed by Spanish accessions. Ancestral group G3 was formed by 13 Italian accessions and one Greek. In ancestral group G4, nine of 12 accessions were from the USA, along with two French and one Italian accession. Ancestral group G5 had six French accessions, three Tunisian, one Greek, one Jordanian, one Iranian, one Moroccan, and one Spanish accession.

Additive kinship results showed five dense clusters (Figure 1C). These clusters included all the accessions forming the five ancestral groups from the population structure analysis. Clusters formed by G1 and G3 had less connections with other accessions. On the other hand, clusters formed by G2, G4, and G5 accessions had several connections between them and accessions from other countries. G2 accessions, separated in two sub-clusters, were also connected with accessions from Australia and North Africa, among others. G4 and G5 clusters were strongly connected, and the G5 cluster was connected to several accessions from Australia.

The phylogenetic tree had three main clades: C1, C2, and C3 (Figure 1A). C1, mainly formed by Asian accessions, was subdivided in two secondary clades, C1.1 and C1.2. All the Chinese accessions were situated in C1.1. C2 was mainly formed by French, American, and Australian accessions. American and Australian accessions were found in the same secondary clade, C2. Based on the Principal Compound Analysis, the ancestral groups were well separated according to the first two principal components, explaining 9.18% of the overall variation (Figure 1B). Admixture accessions were mostly in the center of the scatter plot with values close to 0 for the two principal components. Genotypes from G4 group were located to the right of the X axis while G1, G2, and G3 groups were to the left of the X axis and separated along the Y axis.

Only three breeding cultivars had a F2 value over 0.4

Landraces and breeding cultivars had low levels of inbreeding using F 2 or F 0.25 values (Figure 2A and2C). In the case of F 2 , used to detect recent inbreeding events, inbreeding was low for both landraces and breeding cultivars. Most of the cultivars had an F 2 value equal to zero, with only nine landraces and seven breeding cultivars exceeding F 2 = 0.1.

In the case of breeding cultivars, there were some exceptions with high levels of inbreeding, such as 'Ayles', 'Amandier rose', and 'Garfi', with F 2 values over 0.4. Using F 0.25 as the inbreeding measure, the inbreeding values of both landraces and breeding cultivars were slightly higher, with a mean F 0.25 around 0.1 in both cases. Only accessions 'Ayles', 'Amandier rose', and 'DPRU 487-A' had a F 0.25 value over 0.4.

Comparing Freq 2 and Freq 0.25 there were large differences (Figure 2B and2D). No SNP had a Freq 2 higher than 0.2, while there were three genomic regions with Freq 0.25 higher than 0.4, in chromosomes 1, 4, and 5. In the case of PD01 and PD05, the genomic regions with high Freq 0.25 contained 24 and 27 genes, respectively, according to the 'Texas' reference genome v2.0 [22]. In contrast, the genomic region with high Freq 0.25 in PD04 contained only one gene, Prudul26A014015.

Linkage disequilibrium decay was between 4259 and 6904 bp

The mean size of the LD block was 6061 bp and ranged from 4259 bp (PD08) to 6904 bp (PD02) (Supplementary Material 1).

All the traits had heritability higher than 0.90 A significant negative correlation was found between NW and CRO (-0.77) (Supplementary Material 1). A weaker correlation was also found between NW and KW (0.49). Genetic variance and the variance due to the environment (residual error) were extracted from the linear mixed models for each trait (Table 1). Heritability for all traits was higher than 0.90.

true positive QTLs were found for the traits under study

In total, 19 associations were detected for the traits under study (Figure 3). Of these associations, six were considered false positives, since the phenotypic data distribution did not match the genotype-phenotype interaction searched. Therefore, 13 QTLs were considered true positives (Figure 3, column 3). These QTLs were named according to the recommendations for standard QTL nomenclature and reporting of the Genome Database for Rosaceae (Table 2). Within these 13 QTLs, only qP-CRO6 had an additive effect, the rest had non-additive effects on the phenotype. By trait, we found three, seven, two and one QTL for NW, CRO, DK and BLO, respectively. No QTL was found for KW. NW and CRO had QTLs located at the same position: qP-NW2 and qP-CRO2 had SPD02_19131932 as top SNP with a recessive effect, while qP-NW3 and qP-CRO3 had SPD03_1356626 as top SNP with a overdominant effect. The variance explained by the QTLs ranged from 9% (qP-NW6) to 49% (qP-CRO2). The combined variance explained for all the QTLs was 45.74% for NW, 71.78% for CRO, 65.87% for DK, and the only QTL found for BLO explained 16.44% of the variance.

Candidate genes controlling crack-out percentage, double kernels percentage, and blooming time

The length of the qP-CRO2 region was 13.8 kb, containing just one gene (Table 3). This gene, Prudul26A013473, was annotated as a NAC domain-containing protein.

The length of the qP-DK7.1 region was 11.5 kb, containing two genes (Table 3) encoding eukaryotic translation initiation factor 3 subunit E and an alpha-carbonic anhydrase domain-containing protein.

The length of qP-DK7.2 was 11.5 kb, containing four genes (Table 3) encoding an ML domain-containing protein, a ribosomal L6 domain-containing protein, and two uncharacterized proteins.

The length of qP-BLO2 was 13.8 kb, containing three genes (Table 3) encoding a cyclic nucleotide-binding domain-containing protein, a beta-amylase, and an uncharacterized protein.

The gene coexpression network (GCN) analysis of Prupe.2G196600, the homologous gene of Prudul26A013473, indicated 25 enriched terms, with nine classified as top enriched terms (TET) (Supplementary Material 2; Figure 4A). Within GObp, we found two TETs: 'plant-type secondary cell wall biogenesis' and 'cell wall organization'. Within GOmf, we found two TETs: 'transferase activity, transferring glycosyl groups' and 'O-acetyltransferase activity'. Within GOcc, we found five TETs: 'Golgi apparatus', 'extracellular region', 'apoplast', 'Golgi membrane', and 'trans-Golgi network'.

The GCN analysis of Prupe.7G052700, the homologous gene of Prudul26A012082, indicated 192 enriched terms, and 28 were classified as TETs (Supplementary Material 2; Figure 4B).

Discussion

Almond genetic structure explains its historical worldwide dissemination

Our results strongly supported the subdivision of these accessions into five ancestral groups. Each group was formed by accessions with a common geographical origin: G2 exclusively by Spanish accessions, while G1, G3, G4, and G5 were mainly formed by Chinese, Italian, American and accessions, respectively. These results are in agreement with Pavan et al. 2021, who, with a more limited germplasm, found four ancestral groups, each formed mainly by accessions from Spain, France, Italy and the USA. Apart from Chinese accessions, G1 included accessions from Pakistan, Iran, Turkey, Greece, and Romania. Due to the diverse origin of the accessions forming this group, they may be considered as part of a more primitive almond pool. This agrees with the first dissemination of almond from its center of origin, spreading throughout south-western Asia, eventually reaching modern Turkey, Greece, and other regions of the Eastern Mediterranean [START_REF] Gradziel | Botany, Production and Uses[END_REF].

After reaching the Eastern Mediterranean, Greeks, and Phoenicians introduced the almond to other adapted areas of the Mediterranean. By the time of the Roman Empire, almond cultivation had spread all along the Mediterranean coast. Ancestral groups G2, G3, and G5 may have been established during this period.

G2 was entirely formed by Spanish accessions. The fact that the Spanish accessions were related to Tunisian and other North African accessions could be explained by the introduction of new North African genetic material during the Arab occupation of the Iberian Peninsula.

Even though almond cultivation was introduced in California by the Spanish missions, our results showed a close genetic relationship between Californian and French accessions. This was noticeable in both phylogenetic tree, kinship, and population structure analysis. In the phylogenetic tree, all Californian accessions were located in only one branch, clustering together with the French accessions. In the kinship analysis, both G4 and G5 clusters were connected by several edges, being the clusters with most connections between each other, indicating a relatively recent common ancestor between these two groups. Finally, the population structure analysis included two French accessions, 'Princesse789' and 'A la Dame (CG14)' within group G4. This genetic relatedness between French and Californian accessions is explained by the introduction of French commercial cultivars in California from 1850 to 1900 [START_REF] Wood | Almond Varieties in the United States[END_REF]. In this sense, 'Princesse789' and 'A la Dame (CG14)' may be among the cultivars introduced in California during that period, or at least, close relatives.

G5 formation also included three Tunisian accessions. This close genetic relatedness may be explained by the exchange of material between the two countries during the French occupation of Tunisia from 1881 to 1956.

Finally, Californian cultivars were introduced to adapted regions of the southern hemisphere, including Chile, Argentina, South Africa, and Australia. This explains the genetic relatedness between Californian and Australian accessions.

The cultivated almond shows signs of inbreeding and domestication

Homozygosity analysis was used to study the inbreeding levels in the studied germplasm. Different ROH lengths allowed us to detect modern (ROH 2 ) and ancient (ROH 0.25 ) inbreeding events. In general, all accessions had a low F 0.25 and F 2 value, regardless of whether they were classified as landraces or breeding cultivars. Nevertheless, there were some breeding cultivars with an F 2 value over 0.4, such as 'Ayles', 'Amandier rose', and 'Garfi'. These results are in agreement with [START_REF] Pérez | Pedigree analysis of 220 almond genotypes reveals two world mainstream breeding lines based on only three different cultivars[END_REF], who concluded that breeding practices could be increasing inbreeding levels.

Looking at Freq 0.25 and Freq 2 , we found no region with high Freq 2 . On the other hand, there were several regions with a high Freq 0.25 , in chromosomes 1, 4, and 5. This may respond to selection pressure during the almond domestication process, fixing in the genome alleles of interest. In the case of chromosome 4, the peak of high Freq 0.25 was formed only by one gene, Prudul26A014015. According to phylomeDB, Prudul26A014015 had two homologous genes in Arabidopsis, AT1G21390 and AT1G76980. AT1G21390 was annotated as EMBRYO DEFECTIVE (EMB) 2170, while AT1G76980 was annotated as patatin-like phospholipase domain protein. EMB genes are related with embryo development in Arabidopsis [START_REF] Meinke | Genome-wide identification of EMBRYO[END_REF]. Since the edible commercial product of the almond is the seed, selective pressure may be exerted on a gene related with embryo and seed development.

Non-additive GWAS, short LD decay and GCN analysis allowed the discovery of several candidate genes for breeding traits

Among the 13 QTLs detected in this study, only one had an additive effect. This indicates that non-additive effects could be the main source of genotype-phenotype interactions in almond. Furthermore, due to the similarity between Prunus genus genomes, this phenomenon could be repeated in other cross-pollinating Prunus species. In peach, a self-compatible Prunus species that still maintains an important level of heterozygosity, this may have predominantly led to the fixation of dominant/recessive QTLs and to the selection for heterozygosis in those that are overdominant.

Another remarkable aspect was that most of the QTLs detected in this study were defined by only one SNP (the top SNP). This is caused by the short LD decay found in the population, affecting this study in two different ways: first, it has facilitated the search for candidate genes, since the genomic regions associated with the traits of interest were small and only a few genes were found in these regions (in the case of qP-CRO2, only one gene was found in the QTL region). Second, this could have limited the detection of some regions of interest. In genomic regions with a lower concentration of SNPs, some regions of interest might have been lost because the distance between SNPs was greater than the LD decay. This could be the case for KW, were no QTL was found.

For the rest of the traits under study, we found several QTLs that could be used in marker assisted selection. However, the mathematical models used to calculate the combined variance explained by these QTLs did not take into account possible epistatic interactions between the QTLs. More research is needed before these QTLs can be used in marker assisted selection.

Among the 13 QTLs detected in this study, only qP-CRO2 has been reported in other studies [START_REF] Pavan | Almond diversity and homozygosity define structure, kinship, inbreeding, and linkage disequilibrium in cultivated germplasm, and reveal genomic associations with nut and seed weight[END_REF][START_REF] Sideli | Genome-wide association analysis and validation with KASP markers for nut and shell traits in almond (Prunus dulcis [mill[END_REF][START_REF] Sánchez-Pérez | Mapping major genes and quantitative trait loci controlling agronomic traits in almond[END_REF][START_REF] Goonetilleke | Genotyping by sequencing in almond: SNP discovery, linkage mapping, and marker design[END_REF]. Only one gene was found in the qP-CRO2 region, Prudul26A013473. This gene, annotated as a NAC TF, was homologous to NST1 in Arabidopsis. NST1 has been reported as a key regulator of the formation of secondary cell walls in woody tissues [START_REF] Zhong | A battery of transcription factors involved in the regulation of secondary cell wall biosynthesis in Arabidopsis[END_REF] and specifically associated with secondary cell wall formation within the enb layer in Arabidopsis seeds [START_REF] Mitsuda | NAC transcription factors, NST1 and NST3, are key regulators of the formation of secondary walls in Woody tissues of Arabidopsis[END_REF]. Furthermore, the GCN analysis of Prupe.2G196600, the homologous gene of Prudul26A013473, gave several enriched terms related to secondary cell wall biogenesis and organization. All this evidence indicates that Prudul26A013473 is the gene responsible for qP-CRO2, having a major role in the transcriptional regulation of the almond endocarp lignification.

For double kernels percentage, we found two different QTLs, qP-DK7.1, and qP-DK7.2. In almond, the presence of double kernels in a shell is due to the development and fertilization of two ovules in the ovary when the secondary ovule does not degenerate [START_REF] Gradziel | Shell seal breakdown in almond is associated with the site of secondary ovule abortion[END_REF]. In the case of qP-DK7.1, only two candidate genes were found in the QTL region, Prudul26A012082 and Prudul26A002885. These genes were annotated as a eukaryotic translation initiation factor 3 subunit E (eIF3e) and alpha-carbonic anhydrase domaincontaining protein, respectively. eIF3e has been reported as essential for embryo development and normal plant cell growth in Arabidopsis [START_REF] Yahalom | Arabidopsis eIF3e is regulated by the COP9 signalosome and has an impact on development and protein translation[END_REF][START_REF] Xia | The Arabidopsis eukaryotic translation initiation factor 3, subunit F (AteIF3f), is required for pollen germination and embryogenesis[END_REF]. The GCN analysis of Prupe.7G052700, the homologous gene of Prudul26A012082, gave several terms related to embryo development (Supplementary Material 2). This indicates that Prudul26A012082 is the gene responsible for qP-DK7.1, having a major role in almond ovule and embryo development.

There were four candidate genes in qP-DK7.2 region, Prudul26A029836, Prudul26A017782, Prudul26A005959, and Prudul26A008330. These genes were annotated as a ML domaincontaining protein, a ribosomal large subunit 6 (RL6) and two uncharacterized proteins, respectively. It has been suggested that RL6, among other ribosomal subunits, is essential for embryogenesis in Arabidopsis [START_REF] Romani | Versatile roles of Arabidopsis plastid ribosomal proteins in plant growth and development[END_REF]. The GCN analysis of Prupe.7G052700, homologous gene of Prudul26A012082, showed several terms related to cell cycle regulation and cell development (Supplementary Material 2): Prudul26A017782 is most likely the gene responsible for qP-DK7.2, having a major role in almond ovule and embryo development.

Within qP-BLO2, three genes were found, Prudul26A000954, Prudul26A028547, and Prudul26A019171. These genes were annotated as a cyclic nucleotide-binding domain-containing protein, an uncharacterized protein and a beta-amylase. Both AT1G01710 and AT4G00520, the homologous genes of Prudul26A000954 in Arabidopsis, were annotated as Acyl-CoA thioesterases. These proteins catalyze the hydrolysis of acyl-CoAs to free fatty acids and coenzyme A. During dormancy breaking in perennial fruit trees, reactive oxygen species (ROS) are produced. One of the pathways that produces these ROS starts from fatty acids, beta-oxidated to monosaccharides, and these monosaccharides produce ROS via mitochondrial respiration or are oxidized via the Pentose Phosphate Pathway [START_REF] Beauvieux | Bud dormancy in perennial fruit tree species: a pivotal role for oxidative cues[END_REF]. The GCN analysis of Prupe.2G169700, the homologous gene of Prudul26A000954, gave several terms related with SWI/SNF complex and BAF60 (Supplementary Material 2). SWI/SNF complexes have been shown to participate in the control of f lower development and blooming time [START_REF] Reyes | The many faces of plant SWI/SNF complex[END_REF]. BAF60 is a SWI/SNF subunit, and induces a change at the high-order chromatin level, repressing the photoperiod f lowering pathway in Arabidopsis [START_REF] Jégu | The BAF60 subunit of the SWI/SNF chromatin-remodeling complex directly controls the formation of a gene loop at FLOWERING LOCUS C in Arabidopsis[END_REF]. Prudul26A000954 therefore appears to be the gene responsible for qP-BLO2, having a major role in blooming time through fatty acids metabolism.

Conclusions

In this study, we carried out genetic structure analysis and nonadditive GWAS in a set of 243 almond accessions. Our genetic results agreed with the archaeological and historical evidence that separate modern almond dissemination into four phases: Asiatic, Mediterranean, Californian, and southern hemisphere. Of the 13 QTLs found for the traits of interest, only one had an additive effect, suggesting that non-additive effects could be the major source of genotype-phenotype interactions in almond and other Prunus species. Based on the fast LD decay and the use of the peach GCN we propose four candidate genes for the main QTLs found in this study.

Materials and methods

Plant material and genotyping

We used a diversity panel of 243 accessions from 21 countries and five continents (Supplementary Material 3). Of the 243 accessions, 161 were maintained in the INRAE collection (43.948611 N, 4.808333 E) and 97 in the IRTA collection (41.170723 N, 1.172942 E), with 78 accessions in common at the two locations. DNA of the 180 accessions from the INRAE and IRTA collections was extracted from leaves according to [START_REF] Antanaviciute | An inexpensive and rapid genomic DNA extraction protocol for rosaceous species[END_REF].

After DNA extraction, samples were genotyped using the 60 K SNP array [START_REF] Duval | Development and evaluation of an AxiomTM 60K SNP Array for almond (Prunus dulcis)[END_REF].

For the remaining 63 accessions, genotype information was obtained from two different sources: 45 resequences were from a previous study [START_REF] Duval | Development and evaluation of an AxiomTM 60K SNP Array for almond (Prunus dulcis)[END_REF] and resequences from 18 accessions were downloaded from NCBI (Supplementary Material 3).

Genotypic data filtering and datasets

SNP calling of samples from the DNA libraries was according to [START_REF] Duval | Development and evaluation of an AxiomTM 60K SNP Array for almond (Prunus dulcis)[END_REF]. Only SNPs present in the 60 K almond SNP array were selected (60 581 SNPs). All the samples were merged, giving a dataset with 243 accessions and 60 581 SNPs. These SNPs were filtered following these criteria: i) call rate per sample higher than 82%, ii) call rate per SNP higher than 90%, and iii) minimum allele frequency (MAF) higher than 5%.

Using this initial dataset (Table 4), we calculated the identityby-state, i.e. the number of SNPs with the same allelic state shared between accessions. Accessions with an identity-bystate higher than 98% were declared clonal groups. In total, 22 clonal groups with two or more accessions were detected. Within each clonal group, the accession with the highest number of SNPs was selected (Supplementary Material 3). The remaining accessions were classified as landraces and breeding cultivars based on pedigree information (Supplementary Material 3).

From the initial dataset, three more datasets were created for each analysis in this study. For the genetic structure analysis, only the 152 accessions classified as landraces were selected. We also created two more data sets for GWAS, including only phenotyped individuals for nut traits and blooming time, with 79 and 167 accessions respectively. After selecting the accessions, the datasets were filtered again following the same criteria described above. For the datasets used in GWAS, we included two more criteria: (iv) SNPs with three genotypic classes, (v) minimum genotypic class frequency higher than 5%. The four datasets used in this study are presented in Table 4.

Genetic structure analysis

A population structure analysis, an additive kinship, a phylogenetic tree and a principal component analysis (PCA) were used to determine the genetic structure of the 152 accessions classified as landraces. The population structure analysis was performed using the LEA R package [START_REF] Frichot | Fast and efficient estimation of individual ancestry coefficients[END_REF]. The number of ancestral groups tested were from 1 to 15 with ten repetitions. An accession was considered to belong to an ancestral group when the coefficient of belonging to that specific group was higher than 60%. If an accession did not belong to any ancestral group, it was considered admixed. Additive kinship was estimated with the rrBLUP R package [START_REF] Endelman | Ridge regression and other kernels for genomic selection with R package rrBLUP[END_REF]. The phylogenetic tree was built using the unweighted pair group method with arithmetic mean algorithm included in the ape R package [START_REF] Paradis | APE: analyses of Phylogenetics and evolution in R language[END_REF], and PCA using the factoMineR R package [START_REF] Lê | FactoMineR: an R package for multivariate analysis[END_REF].

Homozygosity analysis

Runs of homozygosity (ROHs) were analyzed in all 243 accessions using the detectRUNS R package (https://cran.r-project.org/ package=detectRUNS). Two lengths of ROHs were analyzed: ROH 2 higher than 4 163 686 bp and ROH 0.25 higher than 520 461 bp (2% and 0.25% of the almond genome size according to 'Texas' reference genome v2.0 [START_REF] Alioto | Transposons played a major role in the diversification between the closely related almond and peach genomes: results from the almond genome sequence[END_REF], respectively). ROH 2 and ROH 0.25 were detected using a window size equal to 20 SNPs. The maximum gap between SNPs was equal to 1 000 000 bp for ROH 2 and 100 000 bp for ROH 0.25 . For every accession, the overall inbreeding values F 2 and F 0.25 were calculated using ROH 2 and ROH 0.25 , respectively. Finally, we calculated the frequencies Freq 2 and Freq 0.25 with which every SNP was located in a ROH 2 and ROH 0.25 , respectively.

Linkage disequilibrium decay

The squared correlated coefficient, r 2 , was estimated in the 152 individuals classified as landraces using VCFTools v0.1.16 [START_REF] Danecek | The variant call format and VCFtools[END_REF]. As it was calculated individually for every chromosome using a 250 000 bp window, the r 2 was calculated for every combination of SNPs within that window. We used a threshold of 0.2 to set the LD decay which was then represented graphically using a loess regression function with a span of 0.1.

Phenotypic data collection and analysis

As nut traits, we phenotyped nut weight (NW), kernel weight (KW), crack-out percentage (CRO) and double kernels percentage (DK). Each accession was phenotyped between 9 and 12 years in the IRTA collection. From each accession, at least 100 mature fruits were randomly collected. The fruit was considered mature when the mesocarp was fully dry and split along the fruit suture and the peduncle was near to complete abscission. Samples were stored at room temperature for at least 2 weeks. After measuring NW, the shells were cracked to measure the weight of the kernels. All weights were measured using an electronic balance. DK was measuring by counting the number of shells containing double kernels. CRO was calculated according to Equation 1:

CRO = (NW -KW) /NW (1)
Blooming time (BLO) was phenotyped for three consecutive years (2020-22) as Julian days when about 5% of f lower buds were fully open for each tree. This trait was measured in the INRAE collection.

Best Linear Unbiased Prediction (BLUP) for NW, KW, CRO, and BLO was estimated for each genotype using a linear mixed model according to Equation 2:

P ijk = μ + Y ik + G j + e ijk (2)
Where P ijk is the phenotypic value (=BLUP) of the kth repetition of the jth genotype in the ith year, μ is the mean value of the phenotypic trait, Y ik is the fixed effect of the kth repetition of the ith year, G j is the random genotypic effect of genotype j, and e ijk is the residual error of the model.

BLUP for DK was estimated for each genotype using a linear mixed model according to Equation 3:

P jk = μ + G j + e jk ( 3
)
where for equation 2, P jk is the phenotypic value of the kth repetition of the jth genotype, μ and G j have the same meanings as in equation 1, e jk is the residual error.

For every trait, broad-sense heritability (h 2 ) was estimated as:

h 2 = σ 2 G σ 2 G + σ 2 ε n
where σ 2 G is the genotype variance, σ 2 ε is the residual variance and n is the mean number of measures.

Genome-wide association study (GWAS)

We explored additive and non-additive genotype-phenotype associations in two different datasets: nut traits and blooming time datasets, with 79 and 167 accessions respectively. For this purpose, we transformed these genotypic datasets as follow. For additive effects, the three possible genotypes of a biallelic marker with a reference allele (a1) and an alternative allele (a2), were written in numeric representation as 1 (a1a1, homozygous for the reference allele), 0 (a1a2, heterozygous) and -1 (a2a2, homozygous for the alternative allele). For dominant effects, genotypes a1a1 and a1a2 have the same effect in the phenotype, so a1a1 and a1a2 were codified as 1 and a2a2 as -1. For recessive effects, genotypes a1a2 and a2a2 have the same effect in the phenotype, so a1a2 and a2a2 were codified as -1 and a1a1 as 1. Note that the dominant and recessive transformations correspond to a dominant-recessive genotype-phenotype interaction, but we had to differentiate the effects of a dominant reference allele or a dominant alternative allele. For overdominant effects, genotypes a1a1 and a2a2 have the same effect in the phenotype, so genotypes a1a1 and a2a2 were codified as 1 and a1a2 as 0 (Supplementary Material 1) [START_REF] Tsepilov | Nonadditive effects of genes in human metabolomics[END_REF].

The mixed model from rrBLUP R package [START_REF] Endelman | Ridge regression and other kernels for genomic selection with R package rrBLUP[END_REF] was used in this study. BLUPs were used as phenotypic data for each trait. For every model, we used three different corrections: including the additive kinship (K), the population structure (Q) or both (K + Q):

Y = μ + Xβ + Qv + Zu + ε ( 4
)
where Y is the vector of phenotypic values, μ the overall mean, X the allelic state matrix, β the allelic effect of each SNP 4, Q is the structural matrix estimated by the LEA R package, v is an effect vector estimated by the model and used as a fixed effect, Z is an incidence matrix linking observations to the vector u that is a polygenic random effect with a covariance structure defined by the kinship (K as previously estimated) u ∼ N (0;2KVg), and ε the residual effect. The choice of each correction was based on the adjustment of the p-values obtained to a uniform distribution as expected under the null hypothesis. The corrected Bonferroni threshold at 5% was used to identify significant association between phenotypic data and genotypic markers.

Before considering any significant genotype-phenotype association found as a QTL, we used visual analysis to confirm that the phenotypic data distribution matched the genotype-phenotype interaction searched (e.g., if an association found with the additive transformation matched an additive phenotypic distribution). We considered any significant genotype-phenotype association matching its phenotypic data distribution as a true positive QTL.

For the QTLs considered as true positives, we assumed the Simple model's R 2 as the variance explained for those QTL. We also calculated the combined variance of the QTLs detected for every trait. In this case, we assumed as the combined variance explained the R 2 of a linear regression using all the top SNPs detected for a trait.

Candidate gene analysis

For every trait, we selected the QTL with the highest -log10(pvalue). If the QTL selected had a -log10(p-value) higher than 6.5, we performed a candidate gene analysis. In the case of DK, we used candidate gene analysis on the two QTLs we found, as both had a -log10(p-value) higher than 6.5.

Every QTL region was defined using the position of the top SNP and the estimated LD decay for every chromosome. The beginning of the QTL was defined as the top SNP position minus the estimated LD decay and the end of the QTL was defined as the top SNP position plus the estimated LD decay. We determined the number of genes located in the QTL region using the 'Texas' reference genome v2.0 [START_REF] Alioto | Transposons played a major role in the diversification between the closely related almond and peach genomes: results from the almond genome sequence[END_REF]. Any gene located less than 2000 bp from the QTL region was included in it, as we considered that the regulatory region of that gene was situated within the QTL region. Then we searched the homologous genes from peach (Prunus persica) and Arabidopsis (Arabidopsis thaliana).

To obtain more information on the function of the most suitable candidate genes, Prudul26A013473, Prudul26A012082, Prudul26A017782, and Prudul26A000954, we used the PeachGCN v1 [START_REF] De Los Cobos | First large-scale peach gene coexpression network: a new tool for predicting gene function[END_REF] for gene coexpression network analysis. We first determined the homologous of our candidate genes in peach, then extracted coexpressing genes. Finally, we performed an enrichment analysis of the coexpressing subnetworks using Gene Ontology (GO) and Mapman ontologies [START_REF] Ashburner | Gene ontology: tool for the unification of biology[END_REF][START_REF] Thimm | Mapman: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes[END_REF]. The significance threshold was held at q-value <0.05. Enriched terms annotating at least 2% of the genes in the coexpressing subnetworks were classified as top enriched terms (TET).

Figure 1 .

 1 Figure 1. Genetic structure analysis. A) Phylogenetic tree. B) Principal components analysis C) Additive kinship. Edges with absolute weight less than 0.05 are not represented. D) Population structure analysis.

Figure 3 .

 3 Figure 3. Genome-wide association analysis. Row A) Recessive Q model for Nut weight. Row B) Over-dominant K + Q model for Nut weight. Row C) Additive Q model for Crack-out percentage. Row D) Dominant Q model for Crack-out percentage. Row E) Recessive Q model for Crack-out percentage. Row F) Over-dominant K + Q model for Crack-out percentage. Row G) Dominant K + Q model for Double kernels percentage. Row H) Recessive Q model for blooming time. Column 1) Manhattan plots. Column 2) Q-Q plots. Column 3) Boxplots of the true positive QTLs.

Figure 4 .

 4 Figure 4. Top enriched terms of the most suitable candidate genes. A) Prupe.2G196600 B) Prupe.2G169700 C) Prupe.7G052700 D) Prupe.7G092200.

  Downloaded from https://academic.oup.com/hr/article/10/10/uhad193/7282279 by INRA Centre Val de Loire user on 19 February 2024

  1, but separated in two tertiary clades, C2.1.1 and C2.1.2. French Downloaded from https://academic.oup.com/hr/article/10/10/uhad193/7282279 by INRA Centre Val de Loire user on 19 February 2024

							Pérez de los Cobos et al. | 3
	A		B	0.20
		0.4
				0.15
	Inbreeding (F)	0.2 0.3	Frequency	0.10
		0.1		0.05
		0.0
		Landraces	Breeding cultivars	PD01	PD02	PD03	PD04	PD05	PD06	PD07	PD08
		Group					Chromosome
	C		D
				0.5
		0.4
				0.4
	Inbreeding (F)	0.2 0.3	Frequency	0.2 0.3
		0.1
				0.1
		0.0		0.0
		Landraces	Breeding cultivars	PD01	PD02	PD03	PD04	PD05	PD06	PD07	PD08
		Group					Chromosome
	Figure 2. Homozygosity analysis. A) Inbreeding using F2. B) Freq2 results. C) Inbreeding using F0.25. D) Freq0.25 results.
	accessions were in two different secondary clades, C2.2 and C2.3,
	along with some Tunisian, Italian and Spanish accessions. The
	Spanish, Italian and Tunisian accessions were in C3. And while
	the Spanish and Tunisian accessions were found in the same
	secondary clade, C3.1, most of the Italian accessions formed
	another secondary clade, C3.2.

Table 1 .

 1 Partition of variance of the traits under study

		KW	NW	CRO	DK	BLO
	Genotype	0.024	2.05	136.62	0.0098	59.91
	Error	0.023	0.43	9.27	0.0038	19.59
	Number of repetitions	12.24	12.33	12.24	9.45	4.05
	Heritability	0.93	0.98	0.99	0.96	0.93

NW: nut weight, CRO: crack-out percentage, DK: double kernels percentage, BLO: blooming time, VE: Variance explained, CVE: Combined variance explained.

Table 2 .

 2 Summary of the QTLs identified, indicating the trait and the QTL, the genetic effect, the correction model used, the closest SNP and its chromosome location, the p-value, the variance explained and the combined variance explained

	Trait	Name	Effect	Correction	Top SNP	Chr	-log10(p-value)	VE	CVE
	NW	qP-NW2	Recessive	Q	SPD02_19131932	2	6.07	33.71%	45.74%
		qP-NW3	Overdominant	K + Q	SPD03_1356626	3	6.13	18.77%
		qP-NW6	Overdominant	K + Q	SPD06_24556800	6	6.1	9.04%
	CRO	qP-CRO6	Additive	Q	SPD06_28386604	6	7.19	22.66%	71.78%
		qP-CRO5.1	Dominant	Q	SPD05_2529870	5	6	23.22%
		qP-CRO1.1	Recessive	Q	SPD01_31737488	1	5.75	47.5%
		qP-CRO2	Recessive	Q	SPD02_19131932	2	8.28	49.13%
		qP-CRO1.2	Overdominant	K + Q	SPD01_39980467	1	7.76	33.44%
		qP-CRO3	Overdominant	K + Q	SPD03_1356626	3	7.37	27.9%
		qP-CRO5.2	Overdominant	K + Q	SPD05_14094224	5	6.23	23.57%
	DK	qP-DK7.1	Dominant	K + Q	SPD07_7765332	7	7.46	30.7%	65.87%
		qP-DK7.2	Dominant	K + Q	SPD07_11278816	7	7.54	47.57%
	BLO	qP-BLO2	Recessive	Q	SPD02_17409544	2	6.99	16.44%	-

K: kinship; Q: population structure, NW: nut weight, CRO: crack-out percentage, DK: double kernels percentage, BLO: blooming time, VE: Variance explained, CVE: Combined variance explained.

Table 3 .

 3 List of candidate genes

	QTL	Candidate genes	Prunus Persica homolog	Arabidopsis thaliana homolog	Uniprot function prediction
	qP-CRO2	Prudul26A013473	Prupe.2G196600	AT3G61910 AT2G46770	NAC domain-containing protein
	qP-	Prudul26A012082	Prupe.7G052700	AT3G57290	Eukaryotic translation initiation factor 3 subunit E
	DK7.1	Prudul26A002885	Prupe.7G052800	AT4G21000 AT4G20990	Alpha-carbonic anhydrase domain-containing protein
	qP-	Prudul26A005959	Prupe.7G092000	AT5G06470	Uncharacterized protein
	DK7.2	Prudul26A029836	Prupe.7G092100	AT5G06480 AT3G11780	ML domain-containing protein
		Prudul26A017782	Prupe.7G092200	AT2G18400	Ribosomal_L6 domain-containing protein
		Prudul26A008330	Prupe.7G092300	-	Uncharacterized protein
	qP-	Prudul26A000954	Prupe.2G169700	AT4G00520 AT1G01710	Cyclic nucleotide-binding domain-containing protein
	BLO2	Prudul26A028547	Prupe.2G169800	AT3G49050	Uncharacterized protein
		Prudul26A019171	Prupe.2G169900	AT2G45880	Beta-amylase

Table 4 .

 4 Description of the four datasets used

	Dataset	N • accessions Accessions included	N • SNPs
	Initial	243	All	54 112
	Structure	152	Classified as landraces	53 985
	Nut traits	79	Phenotyped for nut traits	22 928
	Blooming time 167	Phenotyped for blooming time 16 804

Downloaded from https://academic.oup.com/hr/article/10/10/uhad193/7282279 by INRA Centre Val de Loire user on 19 February 2024

Acknowledgements

We acknowledge Dr. Vincent Segura for his incisive comments reviewing this paper, substantially improving it. We acknowledge financial support from the grants SEV-2015-0533 and CEX2019-000902-S, PID2019-110599RR-I00, and PID2020-118612RR-I00 (Better Almonds) funded by MCIN/AEI/ 10.13039/501100011033; the Grant RTA2017-00084-00-00 funded by MCIN/AEI/ 10.13039/ 501100011033 and by 'ESF Investing in your future'; the Grant PCI2019-103670 funded by MCIN/AEI/ 10.13039/501100011033 and co-financed by the European Union; and to the project ANR-18-PRIM-0001 (FREECLIMB) funded by the French National Research Agency. Finally, the grant of F.P.C. PRE2018-086724 funded by MCIN/AEI/ 10.13039/501100011033 and by 'ESF Investing in your future'. This work has been carried out within the framework of F.P.C.'s PhD in Plant Biology and Biotechnology at the Autonomous University of Barcelona.

Data Availability

All resequences used in this study can be found at NCBI under the following accession numbers: SRR3141040, SRR3141049, SRR3141057, SRR3141065, SRR3141073, SRR3141083, SRR3141098, SRR3141113, SRR3141181, SRR3141192, SRR3141204, SRR3141229, SRR4036105, SRR4045224, SRR4045225, SRR4045227, SRR4045228, SRR4045229.

Author Contributions

F.P.C.: conceptualization, analysis, and writing the original draft. E.C.: conceptualization, analysis and writing the original draft. N.D.: laboratory tasks and analysis. B.Q.T.: supervision, reviewing and editing. I.B.: supervision, reviewing, editing and project management. P.A.: supervision, reviewing, editing and project management. I.E.: conceptualization, supervision, reviewing, editing, project management and seeking funding. H.D.: conceptualization, supervision, reviewing, editing, project management and seeking funding.

Conflict of Interest

The authors declare that they have no conf lict of interest.


Gianpietro Consolaro

Nelson Lossing

Denis Barthou

Zhen Zhang

Corinne Ancourt

Cedric Bastoul

ParameTrick: Coefficient Generalization for Faster Polyhedral Scheduling

Keywords: polyhedral scheduler, polyhedral optimization, integer linear programming

Polyhedral Schedulers have been widely used for loop optimization in general-purpose compilers and, more recently, deep-learning compilers. State-of-the-art scheduling algorithms define a vast space of loop transformations and find an optimal solution according to pre-designed cost functions. However, this Integer Linear Programming (ILP) approach can sometimes have complexity issues. The resolution time of ILP grows rapidly as the size of the problem increases. The complexity of the kernel to optimize, in terms of the number of statements, loops, and dependencies, has a significant impact on solving time. Additionally, the performance of ILP solvers can be severely impacted by big coefficients in the ILP, which mostly come from the loop bounds in the original kernel.

To tackle this issue, this paper introduces a technique called "ParameTrick". This technique is used during the scheduling pipeline to replace some large coefficients with parameters, which are new variables in the ILP model. The approach is analyzed to determine its impact on the solving time, as well as the transformations that are lost and how to limit the number of unfeasible transformations. Results show that ParameTrick leads to faster solving times while the space of dropped desirable solutions is limited.

Introduction

The polyhedral model has been widely used in optimizing compilers. In recent years, it has become fundamental in deep learning frameworks such as Pytorch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] (using Tensor Comprehension [START_REF] Vasilache | Tensor Comprehensions: Framework-Agnostic High-Performance Machine Learning Abstractions[END_REF]) or MindSpore [START_REF] Chen | Deep Learning and Practice with MindSpore[END_REF] (using AKG [START_REF] Zhao | AKG: Automatic Kernel Generation for Neural Processing Units Using Polyhedral Transformations[END_REF]).

Polyhedral optimization can expose parallelism and improve data locality by exploiting different architectural features.

Polyhedral optimization relies on an algebraic representation, using polyhedrons (set of affine constraints) to express concepts such as domain (iteration space), dependencies (ordering relations between different iterations), and scheduling (complete execution order of the domain). This representation makes applying classic loop transformations (skewing, tiling, permutation, etc.) possible using linear algebra transformations.

The main challenge of this technique is to find a transformation that fully exploits architectural features to maximize performance. Polyhedral schedulers (such as Pluto [START_REF] Bondhugula | A Practical Automatic Polyhedral Parallelizer and Locality Optimizer[END_REF][3], islscheduler [START_REF] Verdoolaege | Isl: An Integer Set Library for the Polyhedral Model[END_REF], Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. I. One-dimensional time[END_REF]) use ILP (integer linear programming) to automatically find performant loop transformations that preserve semantics. The complexity of the scheduler consists of defining meaningful affine cost functions that guide the search for an optimal transformation. While this part of the research remains open, we tackle another aspect of the complexity of solving ILP problems in this paper. As polyhedral optimization is part of bigger complex compilation pipelines, it is required to minimize the time to find optimal transformations. ILP solving time can significantly increase with an increase in the number of constraint variables or the coefficient values. In polyhedral optimization, the number of constraints and variables grows with the number of statements and for-loops of the input program. Additionally, big loop-bound values and complex data accesses can result in larger ILP coefficients.

This paper presents ParameTrick, a simple preprocessing step that reduces complexity due to coefficients related to large loop-bounds. It may decrease the number of constraints in ILP problems and restrict the exploration space. Thanks to this trick, we will show how the compilation time can be consistently decreased while minimally limiting the transformation space. Our results first show the compilation speedup obtained in some of the PolyBench [START_REF] Pouchet | Polybench: The polyhedral benchmark suite[END_REF] dataset cases and some real cases from the MindSpore [START_REF] Chen | Deep Learning and Practice with MindSpore[END_REF] pipeline. Furthermore, we show how, thanks to this technique, we can find solutions to problems that would normally become untreatable.

In Section 2, we will briefly describe polyhedral optimization. In Section 3, we describe ParameTrick technique and its limits. In Section 4, we show our result analysis on Polybench (Section 4.1) and some tests coming from MindSpore (Section 4.3). Finally, a summary and some future ideas are proposed in Section 5.

Background

Polyhedral optimization is a technique applied for loop transformations, and it can be divided into three main entities: polyhedral representation, polyhedral Scheduler, and polyhedral code generation.

• The polyhedral representation is an algebraic representation of the kernel to optimize, composed of domain, dependencies, and scheduling functions. • The polyhedral Scheduler is an algorithm that takes the polyhedral representation as input and tries to find a new execution order (transformation) of the loop iterations. This is an automatic algorithm based on ILP formulation • The code generation (not discussed in this paper) focuses on generating the code that respects the transformation found by the polyhedral scheduler.

Polyhedral Representation

The polyhedral representation chooses to represent the forloop-based computations using polyhedrons (expressed as a set of constraints).

In the representation, we can find three main components.

Domain

The domain expresses exactly the set of iterations executed by the kernel. This is represented as a polyhedron, where each constraint is an affine combination of iterators and (if needed) parametric loop bounds. The Domain for a given statement S can be mathematically represented as follows:

D 𝑆 ( ì 𝑁 ) =        ì 𝑖𝑡 𝑀 𝑆 • ì 𝑖𝑡 ì 𝑁 1 ≥ 0        where ì
𝑖𝑡 is the vector of iterators surrounding the statement, ì 𝑁 is a vector of parametric constants (normally used for the loop bounds), and 𝑀 𝑆 is a matrix of integer coefficients defining the polyhedron.

Dependency

Polyhedral data dependencies describe ordering relationships between iterations. When accessing multiple times the same memory location in different iterations, if at least one of the accesses is a write, we have a dependency: a dependency 𝛿 𝑆→𝑅 specifies that the statement S must be executed before the statement R to preserve the semantic.

The dependency is defined on a set of iterations of the two statements. The set of these constraints composes a polyhedron similar to the Domain:

𝛿 𝑆→𝑅 =            ì 𝑖𝑡 𝑆 , ì 𝑖𝑡 𝑅 𝑀 𝑆→𝑅 • ì 𝑖𝑡 𝑆 ì 𝑖𝑡 𝑅 ì 𝑁 1 ≥ 0            where ì 𝑁 is a vector of parametric constants, ì 𝑖𝑡 𝑆 ∈ 𝐷 𝑆 ( ì 𝑁 ) (it- eration of Statement S), ì 𝑖𝑡 𝑅 ∈ 𝐷 𝑅 ( ì 𝑁 ) (iteration of Statement R).

Scheduling Functions

The scheduling function maps each statement iteration to a specific multi-dimensional date. The scheduling defines a total order using a lexicographic representation. Given a statement S, the scheduling function Θ 𝑆 is defined as a combination of scheduling dimension 𝜙 𝑆,𝑖 :

Θ 𝑆 : D 𝑆 ( ì 𝑁 ) ì 𝑖𝑡 → N 𝑚 ↦ → (𝜙 𝑆,0 ( ì 𝑖𝑡) ... 𝜙 𝑆,𝑚-1 ( ì 𝑖𝑡))
where m is the number of scheduling dimensions, and 𝜙 𝑆,𝑖 are defined by:

𝜙 𝑆,𝑖 ( ì 𝑖𝑡) = 𝑇 𝑆,𝑖 • ì 𝑖𝑡 ì 𝑁 1 (1)
where 𝑇 𝑆,𝑖 is the transformation vector, with ì 𝑖𝑡 being the vector of iterators surrounding the statement S and ì 𝑁 the vector of parametric constants of the kernel.

Polyhedral Scheduler

The polyhedral scheduler objective is to find a loop transformation, a scheduling function transformation that reorders the iteration execution, preserving the semantics while optimizing a pre-designed cost function. The cost functions normally focus on performance-related objectives, such as data locality and parallelism.

The Scheduling problem is expressed using ILP representation. Most of the known for-loop transformations can be represented simply through linear transformations of the initial schedule. Different state-of-the-art schedulers, mainly differ because of the different approaches (iterative like Pluto [START_REF] Bondhugula | A Practical Automatic Polyhedral Parallelizer and Locality Optimizer[END_REF][3], isl-scheduler [START_REF] Verdoolaege | Isl: An Integer Set Library for the Polyhedral Model[END_REF], and [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. I. One-dimensional time[END_REF], or one-shot such as [START_REF] Pouchet | Loop Transformations: Convexity, Pruning and Optimization[END_REF]), and different cost functions

The polyhedral scheduler aims to determine the optimal Θ 𝑆 , specifically the 𝜙 𝑆,𝑖 (0 ≤ 𝑖 ≤ 𝑚 -1) that compose it. The variables of the ILP problem are composed of the coefficients of the transformation matrix 𝑇 𝑆,𝑖 shown in Equation 1, plus some cost functions if necessary.

Among the cost functions and constraints defined in state of the art, the validity constraint is the fundamental one because it guarantees the transformation legality (semantic preservation).

Validity Constraint

The validity constraint (introduced by Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. I. One-dimensional time[END_REF]) ensures that 𝑇 𝑆,𝑖 will preserve the semantics of the program. It is defined for each dependency 𝛿 𝑆 ->𝑅 , and it constraints T to be scheduled after S.,

The constraint for the full scheduling function Θ 𝑆 and Θ 𝑇 can be represented as follows:

( ì 𝑖𝑡, ì 𝑖𝑡 ′ ) ∈ 𝛿 𝑆→𝑅 ⇒ Θ 𝑅 ( ì 𝑖𝑡 ′ ) ≻ Θ 𝑆 ( ì 𝑖𝑡)
The ≻ stands for lexicographically greater. Considering that each Θ is composed by a set of 𝜙 𝑖 (0 ≤ 𝑖 ≤ 𝑚 -1), the validity constraint expressed on the single dimensional component 𝜙 𝑖 becomes:

( ì 𝑖𝑡, ì 𝑖𝑡 ′ ) ∈ 𝛿 𝑆→𝑅 ⇒ 𝜙 𝑅,𝑖 ( ì 𝑖𝑡 ′ ) ≥ 𝜙 𝑆,𝑖 ( ì 𝑖𝑡) (2)
that, thanks to the definition given in Equation 1, can be written as:

𝑇 𝑅 𝑖 * ì 𝑖𝑡 𝑅 ì 𝑁 𝑅 1 -𝑇 𝑆 𝑖 * ì 𝑖𝑡 𝑆 ì 𝑁 𝑆 1 ≥ 0 (3)
The two vectors 𝑇 𝑆 𝑖 and 𝑇 𝑅 𝑖 are the variables of our ILP problem, but as we may notice, the constraint is not linear. Farkas lemma must be applied to obtain linear constraints, as shown in [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. I. One-dimensional time[END_REF].

ParameTrick

ParameTrick consists of a simple simplification applied during the Scheduling process. As shown in section 2, the polyhedral Scheduler solves ILP systems, composed of constraints such as the Validity one, coming from the dependencies. The complexity of the ILP depends on several factors:

• number of dependencies: the number of dependencies impacts the number of constraints. The more dependencies we have, the more constraints will compose the ILP. • number of statements and their dimensionality: these two factors directly impact the number of variables of the ILP. • the numerical complexity of the constraints: this depends on the numerical complexity of the dependencies. Big coefficients in the original dependencies may bring many complex constraints, especially considering that Farkas lemma and Fourier Motzkin Elimination are applied to the Validity constraint and other cost functions.

Kernel f o r ( i = 0 ; i <= 5 3 7 ; i + + ) { S : c [ i ] = b ; f o r ( j = 0 ; j <= 5 3 7 ; j + + ) { R : d [ i ] [ j ] = c [ i ] ; } } Domain 𝐷 𝑆 = 1 0 -1 537 𝑖 𝑆 1 ≥ 0 𝐷 𝑅 = 1 0 0 -1 0 537 0 1 0 0 -1 537 𝑖 𝑅 𝑗 𝑅 1 ≥ 0 (a) Initial Kernel

and Domains

Dependency 𝛿 𝑆→𝑅 ParameTrick addresses the third factor, reducing the numerical complexity of the ILP problem. It replaces all the dependency polyhedron's large constant coefficients (loop bounds) with new parameters (that will be added as new columns). This reduces the presence of large coefficients in the dependencies but increases the number of variables in the ILP.

1 0 0 0 -1 0 0 537 0 1 0 0 0 -1 0 537 0 0 1 0 0 0 -1 537 1 -1 0 0 -1 1 0 0 𝑖 𝑆 𝑖 𝑅 𝑗 𝑅 1 ≥ 0 Validity Constraint 0 0 537 -1 1 0 0 0 0 -1 1 0 -537 537 537 -1 1 0 -537 537 0 -1 1 0 𝑡 _𝑖 𝑆 𝑡 _𝑖 𝑅 𝑡 _𝑗 𝑅 𝑡 _1 𝑆 𝑡 _1 𝑅 1 ≥ 0 (b) Before ParameTrick Dependency 𝛿 𝑆→𝑅 1 0 0 0 0 -1 0 0 1 0 0 1 0 0 0 0 -1 0 1 0 0 0 1 0 0 0 0 -1 1 0 1 -1 0 0 0 -1 1 0 0 0 0 0 0 1 0 𝑖 𝑆 𝑖 𝑅 𝑗 𝑅 𝑁 1 ≥ 0 Validity Constraint 0 0 1 -1 1 0 0 0 0 0 0 -1 1 0 0 0 -1 1 1 -1 1 0 0 0 -1 1 0 -1 1 0 0 0 0 0 0 0 0 -1 1 0 𝑡 _𝑖 𝑆 𝑡 _𝑖 𝑅 𝑡 _𝑗 𝑅 𝑡 _𝑁 𝑆 𝑡 _𝑁 𝑅 𝑡 _1 𝑆 𝑡 _1 𝑅 1 ≥ 0 (c) After ParameTrick (N = 537)
If the transformation contains a parametric shifting (where the parameter has been introduced by ParameTrick), we finally substitute it with the original coefficient.

The pre-processing step significantly impacts the solving time of the ILP because it simplifies the constraint system. Additionally, this trick ensures that only actual dependencies are present because we use the real bounds during the dependency analysis, allowing for precise dependency detection. If we substituted the loop bounds before the dependency analysis, we could detect some extra dependencies.

This trick simplifies the ILP solver work but can lead to the exclusion of some transformations. This strongly depends on the constraints that we add for each parameter introduced. The naive method would be to add no information about the parameters. However, adding more constraints (for example, specifying that a parameter is greater than 0 or the relation between different parameters) allows us to find more solutions in general but increases the compilation time.

In Fig. 1, we show a basic kernel with the corresponding domain (a), dependency, and Validity constraint deriving from dependency (b). As visible, even from a trivial dependency as the one in the example, we can see that the corresponding validity constraint is numerically complex. In Fig. 1(c), we show the corresponding dependency after applying ParameTrick (substituting the loop bound 537 with a parameter N ) and the derived Validity constraint. The numerical complexity of the matrices (both the dependency one and the constraint one) decreased. On the other side, the number of variables increased. We need to add the parameters introduced by our technique for the dependency. For the validity constraint, we need to add, for each statement, a variable for each parameter introduced. These variables are used to obtain parametric shifting.

Input Kernel f o r ( i = 0 ; i < 1 5 ; i + + ) { f o r ( j = 0 ; j < 1 0 ; j + + ) { S : A[ i ] [ j ] = 0 ; } } Collapsed Kernel f o r ( i = 0 ; i < 1 5 0 ; i + + ) { S : A[ i / 1 0 ] [ i %10] = 0 ; } Initial Schedule Θ 𝑆 = [𝑖, 𝑗] Collapsed Scheduling Θ 𝑆 = [10𝑖 + 𝑗]
Even though the transformation looks trivial, it can strongly impact the total Scheduling compilation time. The majority of the scheduler time is spent inside the ILP solver. Simplifications as ParameTrick can tremendously decrease the solving time, especially when the input kernel contains many dependencies that could become untreatable when building the ILP problem.

Limitations

ParameTrick objective is to decrease the complexity of the ILP solver, but in reality, this can also reduce the space of possible transformations. In other words, ParameTrick narrows down the set of scheduling transformations that the scheduler can find, similar to how input kernels with only parametric loop bounds would limit the scheduler search space.

The main transformation that is dropped is loop collapse. Focusing on Fig. 2 example, we can see that the collapse scheduling transformation consists of a skewing where the skew coefficient is the initial loop bound.

If we apply ParameTrick, using, for instance, 𝑁 = 15 and 𝑀 = 10, we cannot find internally the same collapse. This should correspond to:

Θ 𝑆 = [𝑀 • 𝑖 + 𝑗]
but 𝑀 •𝑖 is not an affine transformation in the form described in Equation 1.

When we substitute the loop bounds with parameters, we lose several transformations because of the loss of precise information. To preserve some of these transformations, we can introduce simple domain constraints, for instance, 𝑁 > 0, and relational constraints like 𝑁 > 𝑀. However, even with these constraints, it may not always be enough, especially with more complex dependencies. In some cases, it may be necessary to know the exact difference between the parameters or other relations, but this would increase the complexity of the problem, thereby losing some of the benefits of our approach.

Results

PolyBench

In this section, we present the results of our experiments on PolyBench [START_REF] Pouchet | Polybench: The polyhedral benchmark suite[END_REF] benchmark, where we applied ParameTrick to improve the compilation time. We used PolyTOPS [START_REF] Consolaro | PolyTOPS: Reconfigurable and Flexible Polyhedral Scheduler[END_REF] scheduler, specifying the Pluto-like configuration (based on Pluto [START_REF] Bondhugula | A Practical Automatic Polyhedral Parallelizer and Locality Optimizer[END_REF][3]), using the ILP solver provided by isl-0.25 for the experiments.

PolyBench is a benchmark containing 30 kernels from different domains, such as data mining, linear algebra, stencil applications, etc. This benchmark has been extensively used for comparisons in the state-of-the-art schedulers because it contains the most common characteristics of different domains. Notice that, in our benchmarks, we excluded adi, nussinov, deriche, and ludcmp because they contain a reversed loop. Our scheduler, similarly to Pluto, cannot schedule such a case, and it would simply fall back to the initial schedule.

All the PolyBench cases are, by default, parametric. To show the benefits of ParameTrick in our experiment, we applied a preprocessing that substitutes the parametric loop bounds with the corresponding constants (using the values provided by the Large dataset defined in Polybench). This means that the whole pipeline, including the final code generation, is applied to the non-parametric version of the programs. Moreover, ParameTrick introduces parameters only internally to the schedule, but we substitute them with the corresponding constant in case of parametric shifting.

In the results, we decided to show three different cases:

• Original: for this case, we did not apply ParameTrick.

• p-trick: in this case, we applied ParameTrick, adding the inequalities specifying that our parameters are positive.

• p-trick-extra: for this evaluation, we use ParameTrick, adding inequalities on the positivity of the parameters and constraints specifying the relationship (greater or smaller) between the different parameters.

For the experimental evaluation, we used a workstation with Intel Xeon E5-2683 CPU (x86 64), with 2 sockets with 16 cores each (2 threads for each core). 80 MiB of L3 cache. The compiler is gcc-10.5 or dimensionality. We can notice tremendous speedups in cases like 3mm, cholesky, doitgen, floyd-warshall, gemver, lu.

Compilation Time

In some cases, we can notice a slowdown, which can be explained by the fact that the kernels are really simple (and the compilation time is paltry). In these cases, the numerical simplification is unimportant, while adding variables (by adding parameters) has a negative impact.

Furthermore, we can highlight how p-trick performs generally better than p-trick-extra. The last method introduces some extra information about the relation between different parameters, describing the dependency polyhedron more precisely. This allows us to find some solutions that cannot be found without this extra information, but it increases the complexity of the problem.

Execution Time

In the previous subsection, we showed the strong impact of ParameTrick on the compilation time. For a complete analysis, we need to analyze the impact of the scheduling transformation that is found to understand what kind of opportunities are lost.

For these results, we simply compared the execution time of the different transformations found in the three different modalities (original, p-trick, and p-trick-extra). Our plutolike scheduler uses the Proximity cost function and Validity constraint. Similarly to Pluto, it applies some post-processing, such as tiling (using a default size of 32), intra-tile optimization, and wavefront skewing when necessary. The final code is generated using Cloog [START_REF] Bastoul | Code Generation in the Polyhedral Model Is Easier Than You Think[END_REF] using OpenMP to generate parallel loops.

For cases like doitgen, durbin, floyd-warshall, gramschmidt, symm, cholesky the transformation found is different because, without ParameTrick, the scheduler would apply a loop collapse (or, more in general, a skewing depending on the loop bound). For durbin in particular, we can notice a relevant slowdown. In this case, without ParameTrick we can apply tiling and find parallelism, which is not possible when using ParameTrick. For doitgen, similarly to durbin, without ParameTrick, we can find external loop parallelism and tile more dimensions (obtaining bigger loop bands). For correlation, on the contrary, we obtain a significant speedup from some interchanges and different shifting choices. For floyd-warshall and cholesky instead, we obtain good speedup, even if the solution would be suboptimal if we do not apply ParameTrick.

For cases as correlation, 2mm, and atax, p-trick-extra can find the same original solution, but p-trick cannot because there is a shifting that cannot be applied without knowing the relation between the different parameters. Even though the different solution is sub-optimal in general, it is interesting to notice that in correlation, the sub-optimal solution (considering the Proximity cost function [START_REF] Bondhugula | A Practical Automatic Polyhedral Parallelizer and Locality Optimizer[END_REF]) is in practice much faster in execution time. For 3mm, verifying if the transformation found is the same was impossible because the original version takes an indefinite time. The scheduling transformation remains the same for all the other cases (where we obtain a speedup of 1).

Dataset Size Analysis

This section analyses the relationship between dataset size (i.e., loop bounds sizes) and compilation time (excluding ParameTrick). To conduct the tests, we begin with the large dataset sizes as defined in polybench and then increase the loop bounds by a factor indicated on the axis (for instance, 2xlarge is equivalent to large, but with each loop bound multiplied by 2).

In Fig. The size analysis requires further analysis to be explained in detail because, as we can see, the compilation time growth changes depending on the input, suggesting that some dependencies (and some particular polyhedrons) are more variant to the increase of the loop bounds. Furthermore, we think GMP internal representation can play a fundamental role in these slowdowns, but we leave this analysis for the future since the GMP overhead is given by isl. To obtain a precise analysis, this will require a deep look into isl implementation, gathering information about GMP usage for each use case.

We want to remark that we used isl instead of FPL for this analysis because, in this case, isl outperforms FPL. We noticed that while FPL scales better than isl when the ILP number of variables increases, isl can handle the increase of numerical complexity better than FPL.

MindSpore AKG

The ParameTrick is not limited to simple benchmarks such as PolyBench and can be applied in real applications such as AKG [START_REF] Zhao | AKG: Automatic Kernel Generation for Neural Processing Units Using Polyhedral Transformations[END_REF] in MindSpore [START_REF] Chen | Deep Learning and Practice with MindSpore[END_REF]. MindSpore is a deep-learning framework that can be used, similarly to Pytorch and Tensorflow, to define deep-learning models.

AKG is the deep-learning compiler used to lower and optimize the operators of the model. Its optimization relies on polyhedral optimization. Minimizing the compilation time without losing essential optimization opportunities is a key objective in deep learning compilers. A method such as ParameTrick has a lot of value in such an environment because it can tremendously reduce the full compilation time of deep learning models. The tradeoff is minimal, losing almost no important transformations, especially considering the simplicity of deep learning operators, which rarely require complex transformations.

Table 2 shows the benefits obtained using ParameTrick to some realistic use cases, showing noticeable speedups. In these experiments, we used FPL-16.0.6 [START_REF] Pitchanathan | FPL: Fast Presburger Arithmetic through Transprecision[END_REF] as ILP-solver to show that the impact of ParameTrick is preserved even on a more recent solver. The cases shown come from the composition of real deep learning operators from some of the most used deep learning models, such as Resnet, Bert, and Transformers. For some cases, the compilation time is a few hundred milliseconds, but considering that a deep learning model can be composed by thousands of these operators, speedups such as the one showed in 2 are quite impressive. In this scenario, the final scheduling transformation is the same, enabling or disabling ParameTrick, so no further performance analysis is presented.

Conclusion

In this paper, we discussed a simple technique, ParameTrick, that can be applied to mitigate the compilation time spent during polyhedral scheduling. Our method temporarily substitutes numerical loop bounds with parametric ones during the scheduling phase. If the final scheduling transformation contains some parametric solutions, we substitute the parameters with the original constant. This simplifies the numerical complexity of the ILP formulated, achieving great speedup in the overall compilation time for several input cases from PolyBench and, similarly, for some cases from the AKG pipeline. The drawback of this technique is the reduction of the space of possible transformations that may lead to suboptimal solutions. We showed that in PolyBench, our technique can lose some good optimization opportunities, but in most cases, it does not affect the execution performance while drastically reducing the compilation time. Furthermore, thanks to our technique, we can deal with input problems that would be untreatable otherwise, as shown for 3mm from PolyBench and in two cases coming from the AKG scenario. Our work shows the benefits of ParameTrick and indirectly highlights the different optimization opportunities lost when using such a technique or when the input kernel contains parametric bounds.

For future work, we strongly believe that a deep analysis of the impact of different loop bounds on the different ILP solvers' compilation time would be extremely interesting to understand better the benefits of ParameTrick. Moreover, some of the PolyBench examples highlighted how, in some cases, we can obtain better solutions (in execution time) that are suboptimal considering the Proximity cost function. In the future, we would like to investigate these cases to understand the limitations of Pluto's cost function and how to overcome them.

Figure 1 .

 1 Figure 1. Simple kernel and respective domain representations (a), the dependency and the corresponding validity constraints after Farkas-lemma and Fourier Motzkin Elimination (b). Finally, the dependency and the validity constraints after applying ParameTrick (c)

Figure 2 .

 2 Figure 2. Loop Collapsing example: in this figure, we can see the loop collapsing transformation (top) and the corresponding schedule transformation (bottom)

Figure 3 .

 3 Figure 3. Compilation time (seconds), with (in red) and without (blue) ParameTrick, changing the dataset size (loop bounds). For the experiments, isl ILP solver has been used.

Table 1 .

 1 Table. 1 shows the compilation time speedup applying ParameTrick on PolyBench test cases. The most relevant speedups are coming from cases with complex dependencies in general or a high number of statements PolyBench

			Compilation	Execution
	Case	Original Time (ms)	Speedup (original / p-trick)	Speedup (original / p-trick-extra)	Speedup (original / p-trick)	Speedup (original / p-trick-extra)
	3mm	>600000	>3636.32	>176.67	n.a.	n.a.
	correlation	2201.4	2.31	0.09	25.12	1
	cholesky	9032.07	151.43	147.72	2.05	2.23
	lu	14774.03	263.04	291.1	1.61	1.61
	floyd-warshall 11820.16	391.19	386.67	1.27	1.27
	bicg	19.04	0.91	0.73	1	1
	covariance	887.05	5.63	2.41	1	1
	fdtd-2d	650.69	6.86	5.47	1	1
	gemm	22.49	1.01	0.59	1	1
	gemver	3282.38	107.28	106.12	1	1
	gesummv	29.96	0.93	0.93	1	1
	heat-3d	751.84	4.41	4.2	1	1
	jacobi-1d	30.97	1.23	1.22	1	1
	jacobi-2d	351.02	5.33	4.82	1	1
	mvt	11.35	0.91	0.94	1	1
	seidel-2d	51.89	0.97	0.9	1	1
	syr2k	22.25	0.99	0.79	1	1
	syrk	20.93	0.92	0.73	1	1
	trisolv	23.15	1.09	1.08	1	1
	trmm	31.14	1.02	0.61	1	1
	2mm	2523.95	35.93	5.53	0.8	1
	gramschmidt	274.84	1.53	0.32	0.79	0.79
	symm	322.79	4.55	3.5	0.78	0.78
	atax	59.94	2.12	1.28	0.76	1
	doitgen	3984.24	34.59	7.6	0.08	0.08
	durbin	195.04	2.79	2.67	0.00002	0.00002

original compilation time (left) and execution time (right) and speedup using two different versions of ParameTrick. p-trick when enabled with basic constraints. p-trick-extra when enabled with additional constraints between parameters. The highest speedups among the two methods are in bold. The timing is in milliseconds (ms). For these results, isl-solver has been used to solve the ILP problem.

Table 2 .

 2 [START_REF] Uday Bondhugula | PLuTo: A Practical and Fully Automatic Polyhedral Parallelizer and Locality Optimizer[END_REF], we can see that for 2mm, covariance, and gemver, the compilation time is greatly impacted by the loop-bounds size. ParameTrick has the advantage of being completely invariant to this factor, remaining constant for all the possible Compilation time for some AKG[START_REF] Zhao | AKG: Automatic Kernel Generation for Neural Processing Units Using Polyhedral Transformations[END_REF] input cases. Original compilation time, compilation time using ParameTrick (p-trick), and the speedup. FPL-16.0.6[START_REF] Pitchanathan | FPL: Fast Presburger Arithmetic through Transprecision[END_REF] has been used as ILP-solver. The transformation found is the same when enabling/disabling ParameTrick for all the cases, so no performance analysis has been provided.loop bounds values. On the other hand, we can see that cases like fdtd-2d are not really impacted by the loop bounds variation.

	Case	Original Time (ms)	Time (ms) (p-trick)	Speedup (p-trick)
	batch_norm	>600000	9764	>61.45
	two2fractal_v1	96519	78	1237.42
	two2fractal_v2	105	51	2.05
	two2fractal_v3	344	28	12.29
	maxpool_grad_v1	5583	2333	2.39
	force_grad	381	180	2.12
	max_pool_grad_v2 >600000	529	>1134
	hpl_cholesky	9291	121	76.78
	hpl_lu	29396	97	303.05


G Capson-Tojo
email: gabriel.capson-tojo@inrae.fr

D J Batstone
email: d.batstone@awmc.uq.edu.au

T Hülsen
email: t.huelsen@awmc.uq.edu.au

Using a mechanistic model to further develop purple phototrophic bacteria systems for resource recovery

Keywords: resource recovery, purple phototrophic bacteria, wastewater treatment, anoxygenic photosynthesis, purple non-sulfur bacteria, mechanistic modelling

Introduction

Processes based on purple phototrophic bacteria (PPB) are a promising option for resource recovery from wastewater. The photoheterotrophic metabolism of PPB (using light as energy source) allows them to simultaneously remove organics and nutrients (e.g., N and P), recovering them as value-added products (e.g., biomass (to be used as single-cell protein source or as fertiliser), or polyhydroxyalkanoates) at high biomass yields (up to 1 g CODbiomass•g CODremoved -1 ) [START_REF] Capson-Tojo | Purple phototrophic bacteria for resource recovery: Challenges and opportunities[END_REF].

The single mechanistic model available to represent nutrient recovery by PPB is the Photo-anaerobic model (PAnM) [START_REF] Puyol | A mechanistic model for anaerobic phototrophs in domestic wastewater applications: Photo-anaerobic model (PAnM)[END_REF]. The PAnM is accurate to represent laboratory-scale conditions working in a controlled environment, but it cannot predict the behaviour of outdoor PPB reactors treating complex wastewater streams. Consequently, the PAnM has been recently extended to account for: changing environmental conditions, light attenuation, the complex nature of waste streams, microbial synergies and competitions, inihibitory processes, and the diverse metabolic capabilities of PPB [START_REF] Capson-Tojo | Expanding mechanistic models to represent purple phototrophic bacteria enriched cultures growing outdoors[END_REF]. All these factors are essential for the accurate simulation of potential full-scale PPB-based processes for resource recovery.

After calibration and validation with data from a flat-plate demonstration-scale photobioreactor (PBR), the resulting extended model (ePAnM) has been used here to estimate optimal operational and design parameters, and to evaluate scenarios that could foster the implementation of PPB-based systems (e.g., intermittent night aeration). With a clear lack of experimental data from outdoors PBRs, these results might be tremendously useful for the future implementation of PPB reactors.

Methodological Approach

The model was calibrated and validated using data from a 10 m demonstration-scale PPB PBR, with a volume of 1 m 3 (the first of its kind). A period of 33 days, working at varying operating/environmental conditions, was selected (see [START_REF] Capson-Tojo | Expanding mechanistic models to represent purple phototrophic bacteria enriched cultures growing outdoors[END_REF] and Capson-Tojo et al. (2022) for the plant description and the calibration procedure).

The resulting model was used to simulate scenarios to determine optimal operational and desing parameters (e.g., hydraulic retention time (HRT) and reactor width). Simulations assumed an ideal influent with constant composition, with excess of nutrients, and with volatile fatty acids (VFAs) concentrations in the influent of 500-3,000 mg COD•L -1 . Simulation time was 30 days, assuming a reactor equivalent to the aforementioned PBR. Anaerobic, continuously illuminated conditions were assumed (near infrared intensities of 5-200 W•m -2 ), with continuous feeding (pH controlled at 7) and temperatures of 10-45 ºC. Tests at different HRTs assumed a width of 0.08 m and those at varying widths assumed an HRT of 0.5 days. MATLAB (R2015a) was used.

A scenario that could improve the system performance by allowing continuousoutdoors operation was also simulated. This approach (proposed in Hülsen et al. ( 2022b)) combines anaerobic-photosynthetic growth during daytime and aerobicheterotrophic growth during night hours. HRTs of 1 and 2 days were tested, under the same illumination profile used for model validation (real day-night variations). The reactor characteristics were those used in the simulations at different HRTs (influent of 1,500 g COD•L -1 ). At night, dissolved oxygen (DO) levels were 0-5.7 mg DO•L -1 .

Results and Discussion

Figure 1 shows the total VFA removal efficiencies (used as performance crtiteria) from the simulations at different HRTs and reactor widths. The results show that, at optimal temperatures and light intensities (e.g., 25-35 ºC and >25 W•m -2 ) HRTs above 0.5 d must be kept if PPB biomass washout is to be avoided (Figures 1A-B), thus keeping an efficient removal performance. For outdoors systems, this implies that a semicontinuous reactor fed only during daytime will have to run at a minimum HRT of 1 d. This is in agreement with previous limits presented in the literature [START_REF] Capson-Tojo | Expanding mechanistic models to represent purple phototrophic bacteria enriched cultures growing outdoors[END_REF] and with values calculated using common PPB growth rates [START_REF] Capson-Tojo | Purple phototrophic bacteria for resource recovery: Challenges and opportunities[END_REF]. Non-optimal temperatures (10, 45 ºC) resulted in longer required HRTs, as they resulted in slower kinetics. Reducing the light intensities did not jeopardise the process performance until very low values (i.e., below 25 W•m - 2 ). This is obviously also related to the reactor width, which affects light attenuation 8th IWA Water Resource Recovery Modelling Seminar Stellenbosch, South Africa, 15-18 Jan 2023 Affiliation: International Water Association (IWA) Using a mechanistic model to further develop purple phototrophic bacteria systems for resource recovery 3 of 4 through the reactor. Simulations at different reactor widths show that, at common average light intensities and influent concentrations (e.g., 75 W•m -2 and 3,000 mg COD•L -1 ), the performance falls drastically at values over 15-30 cm (assuming 2-sided illumination). This value is also in agreement with numbers from the literature, where it was proven that light attenuation through the reactor is mostly caused by light absorption by pigmented biomass, and that this effect is more pronounced in PPB than in microalgal cultures [START_REF] Capson-Tojo | Light attenuation in enriched purple phototrophic bacteria cultures: implications for modelling and reactor design[END_REF]. Lower incident intensities obviously resulted in shorter optimal widths, and higher influent concentrations lead to more biomass accumulation, favouring light attenuation and reducing the allowable widths.

The simulations results from the day-anaerobic/night-aerobic operation at different HRTs (1 and 2 d) are presented in Figure 2, where the predicted concentrations of VFAs, nutrients, and PPB biomass under different aeration regimes are provided.

Figure 2 Simulation results showing the concentrations of (A) acetate (SAC), (B) other VFA (SVFA), (C) inorganic nitrogen (SIN), (D) inorganic phosphorus (SIP), and (E) PPB biomass (XPB) in an outdoors flat plate PBR working at an HRT of (A1-E1) 1 day, and of (A2-E2) 2 days. Three aeration strategies are compared: no aeration, intermittent aeration during night-time, and continuous aeration. The reactor was continuously fed with an influent containing 600 mg COD•L -1 of total VFAs. A natural illumination profile was assumed, with an incident NIR fraction in the reactor of 100%.

While continuous aeration resulted in PPB outcompetition (and lower nutrient removal due to dominance by heterotrophs) and no aeration resulted in no removal of organics and nutrients during night-time (lack of light), intermittent aeration merged the best of both options. It provided total VFA removal continuously, while keeping an average nutrient removal similar to that achieved in non-aerated reactors (particularly at an HRT of 1 days, thanks to the wash-out of aerobic grazers and PPB aerobic growth).

Conclusion

The ePAnM was used to determine crucial optimal operational and design parameters in PPB flat plate PBRs, concluding that: (i) a minimum HRT of 0.5 d must be kept if PPB biomass washout is to be avoided and, (ii) reactor widths below 15-30 cm should 8th IWA Water Resource Recovery Modelling Seminar Stellenbosch, South Africa, 15-18 Jan 2023 Affiliation: International Water Association (IWA) Using a mechanistic model to further develop purple phototrophic bacteria systems for resource recovery 4 of 4 be used (implying 7.5-15 cm depths for ponds). These values (still unknown in the literature) might be starting points for coming techno-economic analyses, which will determine the feasibility of PPB-based systems for resource recovery. The proposed daily-anaerobic night-anaerobic approach might be an option to improve the system performance, allowing continuous operation, but slightly sacrificing biomass yields due to lower aerobic biomass yields. Further experimental validation is needed.

Opinion

As a technology still in its infancy, there is a current debate around which cultivation strategy (e.g., open ponds or PBRs) should be used in PPB Defenders of using ponds base their reasoning mostly on results from lab-scale systems and on costanalysis using data from microalgae reactors [START_REF] Alloul | Operational Strategies to Selectively Produce Purple Bacteria for Microbial Protein in Raceway Reactors[END_REF][START_REF] Alloul | Volatile fatty acids impacting phototrophic growth kinetics of purple bacteria: Paving the way for protein production on fermented wastewater[END_REF][START_REF] Puyol | The photo-anaerobic raceway for domestic sewage: cradle-to-cradle design in wastewater treatment[END_REF]. Although this approach has merits (and obviously we could be wrong), we think that PBRs might be the right approach to pursue (if not the only feasible one). Fundamental differences between PPB and microalgae (reviewed in Hülsen et al. (2022)) impose that open ponds for PPB need to be too shallow (<15 cm), and thus too expensive. Using results from the presented simulations, we give numbers to the maximum depth that ponds could have, further confirming their limitations. Although dedicated economic analyses should be performed using experimental data and the results presented here, PBRs appear as the sole potentially feasible alternative.

  Stellenbosch, South Africa, 15-18 Jan 2023 Affiliation: International Water Association (IWA) Using a mechanistic model to further develop purple phototrophic bacteria systems for resource recovery 2 of 4

Figure 1

 1 Figure 1 Simulation results showing the total VFA removal efficiencies in a flat plate PBR at different (A) HRTs and temperatures (ºC), (B) HRTs and light intensities (W•m -2 ), (C) reactor widths and light intensities (W•m -2 ), and (D) reactor widths and influent concentrations (mg COD•L -1 ).