Julien Roy

Jean-Marc Bonneville

Patrick Saccone

Sébastian Ibanez

Cécile H Albert

Marti Boleda

Maya Gueguen

Marc Ohlmann

Delphine Rioux

Jean-Christophe Clément

| Cécile

H Albert

Jean- Christophe Clément

Sébastien Lavergne

Roberto A Geremia
email: robertoageremia@gmail.com

Jean-Christophe Clément

Differences in the fungal communities nursed by two genetic groups of the alpine cushion plant, Silene acaulis

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

of cushion plants revealed two genetic clusters matching known subspecies. The exscapa subspecies was found on both limestone and granite, while the longiscapa one was only found on limestone. Even on similar limestone bedrock, cushion soils from the two S. acaulis subspecies deeply differed in their impact on soil abiotic conditions. They further strikingly differed from each other and from the surrounding bare soils in fungal community composition. Plant genotype variations accounted for a large part of the fungal composition variability in cushion soils, even when considering geography or soil chemistry, and particularly for the dominant molecular operational taxonomic units (MOTUs). Both saprophytic and biotrophic fungal taxa were related to the MOTUs recurrently associated with a single plant genetic cluster.

Moreover, the putative phytopathogens were abundant, and within the same genus (Cladosporium) or species (Pyrenopeziza brassicae), MOTUs showing specificity for each plant subspecies were found. Our study highlights the combined influences of bedrock and plant genotype on fungal recruitment into cushion soils and suggests the coexistence of two mechanisms, an indirect selection resulting from the colonization of an engineered soil by free-living saprobes and a direct selection resulting from direct plant-fungi interactions.

K E Y W O R D S

community genetics, fungal community, gene-for-gene interactions, nurse effect, soil ecosystem engineering

INTRODUC TION

Terrestrial plants interact with aboveground and belowground fungal communities [START_REF] Buée | 454 Pyrosequencing analyses of forest soils reveal an unexpectedly high fungal diversity[END_REF][START_REF] Jumpponen | Massively parallel 454 sequencing indicates hyperdiverse fungal communities in temperate Quercus macrocarpa phyllosphere[END_REF] whose assemblage is influenced both by abiotic factors and by the plant community. Through positive or negative plant-soil feedbacks, fungi, in turn, shape plant community composition [START_REF] Bever | Microbial Population and Community Dynamics on Plant Roots and Their Feedbacks on Plant Communities[END_REF], playing a major role in their evolution [START_REF] Brundrett | Coevolution of roots and mycorrhizas of land plants[END_REF] and in ecosystem functioning [START_REF] Berendsen | The rhizosphere microbiome and plant health[END_REF]. Vascular plants provide dead or living organic matter to most soil fungi; thus, they interact directly with biotrophic fungi but also modify their edaphic environment [START_REF] De Deyn | Plant functional traits and soil carbon sequestration in contrasting biomes[END_REF], hence, impacting the saprotrophic community [START_REF] Millard | Does grassland vegetation drive soil microbial diversity?[END_REF].

Distance can limit fungal distributions since they present biogeographic patterns at the continental [START_REF] Prober | Plant diversity predicts beta but not alpha diversity of soil microbes across grasslands worldwide[END_REF][START_REF] Talbot | Endemism and functional convergence across the North American soil mycobiome[END_REF][START_REF] Tedersoo | Global diversity and geography of soil fungi[END_REF] and regional scales [START_REF] Geremia | Contrasting microbial biogeographical patterns between anthropogenic subalpine grasslands and natural alpine grasslands[END_REF][START_REF] Pellissier | Soil fungal communities of grasslands are environmentally structured at a regional scale in the Alps[END_REF].

Locally, the available pool of fungi then undergoes selection by soil-or plant-related parameters, including aboveground plant species composition. Fungal community composition is correlated with soil moisture [START_REF] Hawkes | Fungal community responses to precipitation[END_REF] or pH [START_REF] Rousk | Soil bacterial and fungal communities across a pH gradient in an arable soil[END_REF], but also with plant legacies, such as soil C/N ratio [START_REF] Prober | Plant diversity predicts beta but not alpha diversity of soil microbes across grasslands worldwide[END_REF] or soil organic matter content [START_REF] Zinger | Contrasting diversity patterns of crenarchaeal, bacterial and fungal soil communities in an alpine landscape[END_REF]. The soil fungal communities are further linked to the present plant communities, and fungal beta diversity (compositional dissimilarity between sites, see [START_REF] Anderson | Navigating the multiple meanings of β diversity: A roadmap for the practicing ecologist[END_REF] is correlated with plant beta diversity [START_REF] Geremia | Contrasting microbial biogeographical patterns between anthropogenic subalpine grasslands and natural alpine grasslands[END_REF][START_REF] Pellissier | Soil fungal communities of grasslands are environmentally structured at a regional scale in the Alps[END_REF][START_REF] Prober | Plant diversity predicts beta but not alpha diversity of soil microbes across grasslands worldwide[END_REF][START_REF] Tedersoo | Global diversity and geography of soil fungi[END_REF]. Also, different plant species give rise to diverging fungal communities from the same starting soil in controlled experiments, likely due to differences in rhizodeposition [START_REF] Mouhamadou | Effects of two grass species on the composition of soil fungal communities[END_REF]. Plant intraspecific genetic variation can also lead to divergent composition of dependent communities [START_REF] Bangert | Genetic assembly rules and community phenotypes[END_REF][START_REF] Bangert | A genetic similarity rule determines arthropod community structure[END_REF]. For endophytic fungi, this may be mediated directly by the differential production of defense compounds as shown in maize [START_REF] Saunders | Evidence for alteration of fungal endophyte community assembly by host defense compounds[END_REF] or in Populus [START_REF] Bailey | Host plant genetics affect hidden ecological players: Links among Populus, condensed tannins, and fungal endophyte infection[END_REF][START_REF] Lamit | Tree genotype and genetically based growth traits structure twig endophyte communities[END_REF]. However, little is known about the influence of plant genetic variation on fungal community composition in soils, which is challenging to disentangle from confounding factors in natural environments. Here, we report on a field survey providing insights into this influence.

Alpine cushion plants offer an interesting system to study fungal community assembly and to assess the role of intraspecific plant genetic variation. These foundation plants locally modify their physicochemical environment [START_REF] Badano | Assessing impacts of ecosystem engineers on community organization: A general approach illustrated by effects of a high-Andean cushion plant[END_REF] and the soil biotic community [START_REF] Roy | Microbes on the cliff: Alpine cushion plants structure bacterial and fungal communities[END_REF], both defining ecological engineering effects (see [START_REF] Jones | Organisms as Ecosystem Engineers[END_REF]. Indeed, once a plant is established in cracks on bare mineral soil, the accumulating litter material below the dense green tissues allows the constitution of a de novo organic soil.

Through the buffering of microclimatic conditions [START_REF] Molenda | The Alpine Cushion Plant Silene acaulis as Foundation Species: A Bug's-Eye View to Facilitation and Microclimate[END_REF] and the improvement of soil organic matter and nutrients status [START_REF] Roy | Microbes on the cliff: Alpine cushion plants structure bacterial and fungal communities[END_REF], cushion plants act as nurse species for microbiota [START_REF] Roy | Microbes on the cliff: Alpine cushion plants structure bacterial and fungal communities[END_REF], soil arthropods [START_REF] Maillet | Compareads: Comparing huge metagenomic experiments[END_REF][START_REF] Molenda | The Alpine Cushion Plant Silene acaulis as Foundation Species: A Bug's-Eye View to Facilitation and Microclimate[END_REF], and other plants [START_REF] Butterfield | Alpine cushion plants inhibit the loss of phylogenetic diversity in severe environments[END_REF]. To our knowledge, however, the extent of the abiotic and biotic soil changes induced by various cushion plant species or subspecies has not been explicitly studied.

Silene acaulis (L.) Jacq is a common circumboreal and alpine cushion plant with slow growth and cushions that can often be older than 300 years [START_REF] Morris | Life history of the long-lived gynodioecious cushion plant Silene acaulis (Caryophyllaceae), inferred from size-based population projection matrices[END_REF]. Its xenogamic reproduction makes that two distinct plant individuals always represent distinct genotypes. Although the evolutionary history of Silene acaulis remains to be fully unraveled, two main subspecies have been proposed in the Alps, S. acaulis exscapa and S. acaulis longiscapa, exhibiting dense-(shorter internodes) and loose-cushion morphologies, respectively [START_REF] Fe R E N C E S Aeschimann | Flora alpina[END_REF][START_REF] Bock | Quelques aspects de la variation infraspécifique chez le Silene acaulis (L.) Jacq. Compt-Rendus 108è[END_REF]. They occur on adjacent territories, probably with genetic barriers, longiscapa being restricted to calcareous bedrock while exscapa is mainly, but not only, found on siliceous bedrock (Sébastien Ibanez, S. I., unpublished obs.). Subspecies determination can, however, be difficult since intermediate cushion architectures may sometimes be found (Figure 1a).

We previously reported that the presence of a cushion affects soil chemical parameters on two bedrocks and that this shift in abiotic properties comes along with a shift in the fungal community; in particular, the fungal turnover from bare soil to cushion soil increased with environmental stress (high elevation and very low pH; [START_REF] Roy | Microbes on the cliff: Alpine cushion plants structure bacterial and fungal communities[END_REF]. In the present study, we first established that the two plant subspecies corresponded to distinct genotypic clusters, which then allowed us to assess the correlation between plant genetic distances and fungal beta diversity. To disentangle this genetic effect from those of the local environment (geology, elevation), and also to measure abiotic soil engineering by the plant, we sampled soils both beneath cushions and in neighboring bare soils. We finally examined the fungal functional guilds across the sampling design to address the distribution, filtering, and recruitment mechanisms.

MATERIAL S AND METHODS

Plant and soil sampling

The sampling area was located in the French Hautes-Alpes along replicated ecotones from alpine to subnival environments and encompassed two neighboring mountain massifs: one mostly calcareous (Cerces) and one siliceous (Combeynot; see Roy et al., 2013 for details). In September 2009, at 19 sites located across six elevation gradients (Figures 1b and2a, Supporting Information Table S1), with each site defining a plant population, we sampled cores of organic soil (5 cm in diameter x 5 cm deep) beneath five S. acaulis cushion plants, and five soil cores on neighboring bare ground. Both types of soil cores were split into two aliquots. One aliquot was frozen for physicochemical characterization. A second one was kept in silica gel until DNA extraction for fungal sequencing. Sites were revisited in July 2010 for plant leaf sampling.

Silene acaulis is a gynodioecious species, but plant gender was not recorded since some individuals were not flowering. Vegetative traits are not distinguishable between plant genders, and pollinators visit both plant genders, thus disseminating the same fungal spores; therefore, we assume that gender will be only a weak source of fungal community variation compared to plant genotype or environment. Plant material was kept in silica gel until DNA extraction for genotyping. Altogether, the sampling consisted of 95 leaf samples and 95 pairs of cushion and bare soil samples. Whereas the Combeynot massif is entirely siliceous (i.e., granite and gneiss), the Cerces massif is mostly composed of carbonated sedimentary rocks (hereafter referred to as limestone), but also contains quartzite (Figure 1b). Even though quartzite and granite are both siliceous bedrocks, the pH of quartzite samples was atypical (Supporting Information Table S1), and the two populations on quartzite were withdrawn from most analyses, leaving nine populations on granite and eight on limestone.

Soil physicochemical analyses

Soil chemical parameters were quantified following the methods described in [START_REF] Clément | The effects of snow-N deposition and snowmelt dynamics on soil-N cycling in marginal terraced grasslands in the French Alps[END_REF]. Briefly, soil samples were sieved at 2 mm and used to estimate soil pH in water, gravimetric soil water content [START_REF] Robertson | Standard soil methods for long-term ecological research[END_REF]

Plant genotyping

Plants genotypes were fingerprinted using AFLP [START_REF] Meudt | Almost forgotten or latest practice? AFLP applications, analyses and advances[END_REF] by the procedure of [START_REF] Vos | AFLP -A new technique for DNA-fingerprinting[END_REF] AFLP profiles were analyzed using Peak Scanner software (Applied Biosystems); the profiles, for which the automated process failed to detect or attribute peak size, were processed manually. The binning and scoring of the profiles were performed using the RawGeno package [START_REF] Arrigo | Evaluating the impact of scoring parameters on the structure of intra-specific genetic variation using RawGeno, an R package for automating AFLP scoring[END_REF] in R software [START_REF] R_Development_Core_Team | R: A language and environment for statistical computing[END_REF]. We kept individuals within the 5%-95% quantile confidence interval around the peak number mean and kept the peaks that corresponded to fragments longer than 100 bp to limit size homoplasy [START_REF] Arrigo | Evaluating the impact of scoring parameters on the structure of intra-specific genetic variation using RawGeno, an R package for automating AFLP scoring[END_REF][START_REF] Pompanon | Genotyping errors: Causes, consequences and solutions[END_REF]. Bins that were not reproducible at least once within a group of replicates (four groups of replicates for each pair of primers and three to four replicates per group) were systematically removed. We calculated the "bin content information" criterion [START_REF] Arrigo | Evaluating the impact of scoring parameters on the structure of intra-specific genetic variation using RawGeno, an R package for automating AFLP scoring[END_REF][START_REF] Pompanon | Genotyping errors: Causes, consequences and solutions[END_REF] for various filtering strategies and found that our conservative strategy retained the most information. A total of 93 individual plants were reliably scored for the presence or absence of 345 reliable and polymorphic loci, which were coded in a binary matrix.

The plant population genetic structure was inferred from AFLP data using STRUCTURE 2.3.2 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF], implementing an admixture model without a priori knowledge of the geographic provenance of individuals. The likelihood of the number of genotypic clusters (K) was estimated for K ranging from 1 to 20 with five runs for each K value with 50,000 burn-in periods followed by 50,000-step Markov chain Monte Carlo (Supporting Information Figure S1). Pairwise genetic distances were assessed with the Jaccard index.

Fungal ITS1 sequencing, clustering, and assignment

The ITS1 DNA region was amplified by PCR using DNA extracted from each soil core as previously described [START_REF] Roy | Microbes on the cliff: Alpine cushion plants structure bacterial and fungal communities[END_REF]. Fungal DNA was amplified using the ITS5 (5′-GGAAGTAAAAGTCGTAACAAGG-3′) and the 5.8S_fungi primers (5′-CAAGAGATCC′GTTGTTGAAAGTT-3′) [START_REF] Epp | New environmental metabarcodes for analyzing soil DNA: Potential for studying past and present ecosystems[END_REF] Bioinformatic processing of sequenced amplicons used the OBITools suite [START_REF] Boyer | Obitools: A unix-inspired software package for DNA metabarcoding[END_REF]. Reads were quality-filtered, and ITS amplicons were reconstructed by merging the 3′ and 5′ pair-end reads using solexapairendnull. We used obigrep to discard merged reads with an alignment score ≤20 (corresponding to a perfect match on ≤5 bases for reads with the highest Phred score) together with those containing ambiguous nucleotides, or errors in the primer sequence or the sample tag. The remaining merged reads were assigned to their sample of origin by their terminal tags using ngsfilter, which trims primers and tags away from the sequence, and identical reads were counted within each sample using obiuniq. Finally, singletons were removed following Coissac, Riaz, and Puillandre (2012), together with sequences shorter than 66 nt (the smallest one in our taxonomic database, see below), and an abundance table was built using obitab.

To cluster sequences into MOTUs, we first aligned them pairwise using a global Needleman-Wunsch algorithm implemented in the Sumatra package (https://metabarcoding.org/sumatra). Each pair yielded a similarity index equal to the longest common subsequence normalized to the alignment length. Indices were then used to build clusters sharing at least 98% similarity [START_REF] Lentendu | Assessment of soil fungal diversity in different alpine tundra habitats by means of pyrosequencing[END_REF] using the nonhierarchical Markov clustering method (MCL; [START_REF] Van Dongen | Graph clustering by flow simulation[END_REF][START_REF] Zinger | Microbial diversity in alpine tundra soils correlates with snow cover dynamics[END_REF].

The taxonomic assignment of MOTUs used a custom ITS1 database derived from the UNITE database (release V7, 2017-01) using ecoPCR, which extracts sequences flanked by two input primers. The most represented sequence of each MOTU was aligned against this ITS1 database using the fasta35 algorithm [START_REF] Pearson | Flexible sequence similarity searching with the FASTA3 program package[END_REF] implemented in the ecotag program, and only environmental sequences identical at a minimum of 80% to a reference sequence were assigned. Possible lifestyles were deduced from the MOTU taxonomy using FUNGuild [START_REF] Nguyen | FUNGuild: An open annotation tool for parsing fungal community datasets by ecological guild[END_REF].

Nine samples with fewer than 900 reads were removed, and MOTUs with less than 10 reads or unassigned to the Fungi kingdom were discarded. The final dataset contained 175 samples and 1,311 fungal MOTUs across 1,016,938 reads. The sequences and metadata of the samples are available from EMBL under ENA entry PRJEB14389.

Statistical analyses

MOTU read counts per sample were transformed to relative abundance. To assess the compositional changes between soil fungal communities, we summarized pairwise Bray-Curtis dissimilarities by non-metric multidimensional scaling (NMDS), using vegan [START_REF] Oksanen | vegan: Community Ecology Package. R package version 2[END_REF]. Environmental fitting was used to visualize the correlation between variation in soil parameters and variations in the fungal community. Dissimilarities between cognate bare and cushion soils were calculated by comparing each cushion soil to each bare soil within each sampling site.

To rank the predictors of community composition in cushion soils, we implemented a variance partitioning approach based on dissimilarity matrices in which fungal pairwise dissimilarities were measured using the Bray-Curtis index, either on fungal frequencies or on their square roots; the latter (Hellinger transformation) minimizes the impact of the most abundant MOTUs.

Plant genetic dissimilarities were calculated using the Jaccard index on AFLP fingerprints. Additionally, Euclidian distances were used for soil abiotic environment (using reduced chemical data combining nitrogen, carbon, water, nitrate and ammonium contents, and pH values), both beneath and outside cushions and for geographic positions. Instead of using a traditional regression analysis of variance, we chose a generalized additive model (GAM) that deals with potential nonlinear relationships between pairwise beta-diversities and environmental distances [START_REF] Ferrier | Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment[END_REF]; see Supplemental data) and thus increases the predictive power. The pure and joint explained variance of each predictor (either plant dissimilarities or environmental distances) was calculated using Nagelkerke's R-square:

The pure fraction represents the effect of a variable after accounting for the effects of all other predictors. We also extended the general formula for two variables to the case of four variables and represented the pure and joint explained variance using a Venn diagram.

To test for significant differences between taxonomic abundances beneath and outside cushions, we fitted two linear mixed models using lme4 [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF]: a null model and a model including, as a fixed factor, the location beneath cushion or in bare soil. The two models included the sampling site as a random factor. The models were then compared in an ANOVA test, and p-values were corrected for multiple testing using the false discovery rate [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF]. Fungal diversity was quantified using the inverted Simpson index.

To identify MOTUs that were recurrently associated with plants,

we grouped samples into six habitats combining the granite/limestone bedrock type, the local plant genetic cluster (see Results), and the beneath/outside cushion location. For each MOTU, we counted occurrences over a 0.1% frequency threshold in each habitat and performed a chi-square test of independence. The p-values were obtained through 10,000 permutations of the original matrix and corrected for multiple testing. We kept those MOTUs that were significantly and positively associated with at least one cushion habitat (p-value < 0.05; Pearson residual > 2) but with no bare soil habitat.

RESULTS

Genetic structure of Silene acaulis populations

Two genetic groups of plants were revealed by population structure inference (Supporting Information Figure S1) and the ordination of distances between AFLP fingerprints (Figure 2b). Most samples had at least 80% of their AFLP markers belonging to a single genotypic group (Figure 2a). This criterion allowed us to define two plant genetic clusters and left only four ambiguous plant individuals. The two genetic clusters matched with subspecies determination (Figure 2a), namely the dense cushion exscapa and the loose-cushion longiscapa (Figure 2), and they were named accordingly.

At each of the 19 sampling sites, only one genetic cluster was found, except for the two sites containing ambiguous individuals (Figure 2a). The longiscapa genetic cluster only occurred on calcareous bedrock (limestone), except for one population found on quartzite (CIIIQ, Figure 2a) for which the abiotic environment was close to that of the calcareous environments (Supporting Information Table S1; Figure 3a). By contrast, the exscapa genetic cluster was found on both calcareous and siliceous (granite and quartzite) bedrocks.

Although exscapa culminated at higher elevations than longiscapa, six populations of exscapa and three populations of longiscapa occurred within the same elevation range (2,500-2,700 m a.s.l.; Supporting Information Table S1, Figure 2a).

Differences in soil abiotic engineering by the two plant genetic clusters

The PCA ordination of soil chemical properties revealed two axes totalizing 79% of the total variability (Figure 3a). Axis 1 was related to nutrient and water contents, and Axis 2 was related to soil pH and nitrate content. Cushions from exscapa genetic cluster induced a strong shift of soil chemical properties along Axis 1 from bare to cushion soils. This shift was smaller for longiscapa genetic cluster (Figure 3a,b). Supporting Information The contents in nitrogen, carbon, water, and ammonium all increased on both granitic and calcareous bedrocks beneath exscapa plants, and soil became less acidic on granite (Supporting Information Figure S3).

For soil beneath longiscapa cushions, the increase in nutrient contents was weaker, although significant for the carbon and ammonium contents (Supporting Information Figure S3).

Differences in the engineering of fungal communities by the two plant genetic clusters

NMDS ordination of fungal community dissimilarities revealed a sharp difference in community composition on the two main geological substrates, which correlated with pH and nitrate contents (Figure 3c). On top of this bedrock-related divergence, fungal communities of exscapa cushions from both mountain ranges segregated away from all other habitats, which correlated with increased soil contents in carbon, total nitrogen, ammonium, and water. For longiscapa, the shift in community composition differed from the one observed for exscapa on limestone (Supporting Information Figure S2). From bare to cushion soils, fungal dissimilarities were strong for both subspecies (Figure 3d). The shifts in community composition included a consistent decrease in fungal diversity beneath the cushions compared to bare soils (Figure 4a). This diversity reduction came along with a strong MOTU turnover, which was modulated differently for the two genetic clusters and the two bedrock types (Figure 4b). F I G U R E 4 Fungal diversity and turnover across habitats. Samples were grouped into six habitats combining the plant genotype ("x" for exscapa and "l" for longiscapa), the bedrock type ("G" for granitic and "C" for calcareous), and the location ("b" for bare soil and "c" for cushion soil). From left to right, n = 40, 43, 12, 11, 22, and 20. (a)

Links between plant genetic distances and fungal beta diversity

The variations in the composition of fungal communities among cushion soils were best explained by the combined variations in plant genotype, soil chemical composition beneath and outside the cushions, and geographic position (Figure 5), which accounted together for 40% of the total variability. This variability was 2% more than when the geographic position was dismissed at the profit of elevation (Supporting Information Figure S4). The largest pure effect (10%) was that of the genetic distances between plants (Figure 5), and large plant genetic distances reflected chiefly differences between individuals of different genetic clusters (Figure 2b). The largest shared part (8%) was explained together by plant genetic and geographic distances (Figure 5), reflecting the geographic separation of the two plant clusters. When using a metric minimizing the impact of the most abundant MOTUs, the pure effects of geographic and genetic distances were equal (Supporting Information Figure S4). This discrepancy suggests that each plant cluster selectively recruited a few MOTUs that rarely occur in bare soils but became highly abundant beneath cushions.

Taxonomic identification of fungal sequences associated with the cushion habitats

Ascomycota largely dominated all fungal communities but decreased beneath cushions vs. bare soils, while Basidiomycota significantly increased (Figure 6a). The latter included a trend for Sebacinaceae to increase (Figure 6b; 2.7% beneath vs. 0.6% outside). In both phyla, the fraction of fungi unassigned at the family level was high in bare soils and decreased strongly in cushion soils (Figure 6b). Remarkably, the abundance of Cladosporiaceae and Dermateaceae (Ascomycota) increased largely and significantly beneath the cushions, cumulating in 19.4% of the reads in cushion soils vs. 6.6% in bare soils. The nearly doubled frequency of Dermateaceae was entirely caused by MOTUs assigned to Pyrenopeziza brassicae (Supporting Information Table S2).

We observed a strong turnover of fungal MOTUs across the six habitats (exscapa cushions on granite and limestone, longiscapa cushions on limestone, and their bare soil cognates). Among

DISCUSSION

We have studied the fungal communities of a foundation cushion plant in a microevolutionary context across contrasting environmental conditions, such as different bedrocks. Genotyping of S. acaulis cushions delineated two genetic clusters, one of which can be found both on siliceous and limestone bedrock, while the other one is restricted to limestone. This unique situation, despite asymmetries in the dataset, allows highlighting the key influence of plant genotypes on soil engineering, both in abiotic modifications and in the filtering of fungal species.

F I G U R E 5 Predictors of variation of the fungal communities hosted by cushions. Pairwise dissimilarities in a set of 64 cushions were used to partition the variance of fungal communities using GAM (see Methods). The pure and joint effects of the different explanatory variables on fungal variation are indicated

The two plant subspecies differentially engineer their soil environment

As previously reported [START_REF] Roy | Microbes on the cliff: Alpine cushion plants structure bacterial and fungal communities[END_REF], cushions engineer a soil enriched in carbon and ammonium, two key soil nutrients that are mainly derived from plant litter; here, we found that this edaphic improvement is clearly stronger for the exscapa cluster. These differences in abiotic soil engineering are likely to derive from differences in cushion architecture. A causal link between the nurse plant stature and its engineering effects has been demonstrated in a dominant coastal shrub [START_REF] Crutsinger | Genetic variation within a dominant shrub structures green and brown community assemblages[END_REF]. The exscapa subspecies displays a denser architecture than the longiscapa one, limiting litter dispersion and certainly providing a higher capacity to limit water loss. Differences in cushion architecture between subspecies most likely have a heritable basis because they correspond to distinct genotypic groups and remain observable over a broad elevation range and on distinct bedrocks. Our sampling, however, does not encompass the whole ecological range of each subspecies and therefore does not explore their entire vegetative trait plasticity [START_REF] Bonanomi | Cushion plant morphology controls biogenic capability and facilitation effects of Silene acaulis along an elevation gradient[END_REF], making it difficult to tease apart the respective influences of phenotypic plasticity and genetic differentiation.

Fungal communities strongly varied from bare soil to cushion soils, pointing to the selective filter exerted by cushions on fungi.

This filtering is highlighted first by a decrease in fungal α-diversity and by strong shifts in fungal composition. These changes were again stronger for the exscapa cluster and seem to be related to the plant subspecies rather than to the environment because they occurred with the same amplitude on the two contrasted bedrocks where exscapa plants were found. S. acaulis cushions improve edaphic properties but can also buffer soil moisture and alter the temperature regimes [START_REF] Bonanomi | Cushion plant morphology controls biogenic capability and facilitation effects of Silene acaulis along an elevation gradient[END_REF][START_REF] Molenda | The Alpine Cushion Plant Silene acaulis as Foundation Species: A Bug's-Eye View to Facilitation and Microclimate[END_REF]. This new microenvironment may indirectly favor the growth of new, especially saprophytic, fungal species. There is also a clear influence of the bedrock, as illustrated by the exscapa cushions, which remained largely different in fungal composition on granite or on limestone.

Bedrock, as seen by bare soil chemistry, contributed little as a pure part to the variability of fungal composition but much more to the fractions shared with a geographic position, owing to the spatial separation of limestone and granite sites.

Indirect and direct mechanisms of fungal recruitment by cushions

For the three cushion habitats (exscapa on granite and limestone, and longiscapa on limestone), comparison to bare soils revealed both a strong turnover for putatively saprophytic and biotrophic fungi. In a single host could respond to different degrees of soil engineering from the plant litter, be it in the content in organic matter, nutrient [START_REF] Hanson | Fungal Taxa Target Different Carbon Sources in Forest Soil[END_REF][START_REF] Treseder | Fungal Traits That Drive Ecosystem Dynamics on Land[END_REF], or water [START_REF] Crowther | Untangling the fungal niche: The trait-based approach[END_REF]. For instance, Geomyces (Ascomycota) possibly bloomed beneath exscapa cushions because they occupy the highest elevation on limestone since this genus includes psychrophilic species able to profit from a variety of organic substrates [START_REF] Hayes | The Geomyces Fungi: Ecology and Distribution[END_REF]. Saprobes could also be host-specific epi-or endophytes turning into decomposers upon plant leaf senescence [START_REF] Osono | Phyllosphere fungi on leaf litter of Fagus crenata: Occurrence, colonization, and succession[END_REF][START_REF] Voriskova | Fungal community on decomposing leaf litter undergoes rapid successional changes[END_REF][START_REF] Zhou | Host-specificity, host-exclusivity, and host-recurrence in saprobic fungi* *Paper presented at the Asian Mycological Congress 2000[END_REF].

Plant genetic distance was the best predictor for fungal variations for major MOTUs; therefore, plant genetic identity directly explained a part of the variance in community composition that was not explained by ecosystem engineering. For less abundant MOTUs, the geographic distance was, however, as strong a predictor as plant genetic distance, possibly reflecting stochastic recruitment processes [START_REF] Dumbrell | Idiosyncrasy and overdominance in the structure of natural communities of arbuscular mycorrhizal fungi: Is there a role for stochastic processes[END_REF]. The plant genetics direct effect likely exerts through the symbiotic and pathogenic fungal communities. Among symbionts, we observed Sebacinae (Basidiomycota), a family of root endophytic and ectomycorrhizal fungi [START_REF] Weiß | Sebacinalesone thousand and one interactions with land plants[END_REF]. The former type of association is more likely since S. acaulis is poorly mycorrhized [START_REF] Cripps | Distribution of Mycorrhizal Types among Alpine Vascular Plant Families on the Beartooth Plateau, Rocky Mountains, U.S.A., in Reference to Large-Scale Patterns in Arctic-Alpine Habitats[END_REF][START_REF] Schmidt | Mycorrhizal and Dark-Septate Fungi in Plant Roots Above 4270 Meters Elevation in the Andes and Rocky Mountains[END_REF].

Among putative pathogenic fungi, we observed opposite plant subspecies specificities within MOTUs belonging to the same genus, Cladosporium [START_REF] Bensch | The genus Cladosporium[END_REF], or to the same species (P. brassicae, Dermateceae; [START_REF] Rawlinson | Taxonomy and biology of Pyrenopeziza brassicae sp.nov. (Cylindrosporium concentricum), a pathogen of winter oilseed rape (Brassica napus ssp. oleifera)[END_REF]. MOTUs assigned to Cladosporium and Dermateaceae have been detected within the stems and leaves of S. acaulis in the Arctic [START_REF] Zhang | Endophytic Fungal Communities Associated with Vascular Plants in the High Arctic Zone Are Highly Diverse and Host-Plant Specific[END_REF]. This strongly suggests that direct interactions between the fungus and living plant tissues account for their recruitment in alpine cushion soils and is consistent with the strong costructure between phyllosphere fungal composition and intraspecific population genetic structure in beeches [START_REF] Cordier | Spatial variability of phyllosphere fungal assemblages: Genetic distance predominates over geographic distance in a European beech stand (Fagus sylvatica)[END_REF] as well as with the phylogenetic signal in the association of tropical trees with fungal pathogens [START_REF] Gilbert | Phylogenetic signal in plant pathogen-host range[END_REF]. These narrow host ranges might result from a local co-evolutionary arm race [START_REF] Chappell | Evolution of host range in Coleosporium ipomoeae, a plant pathogen with multiple hosts[END_REF] or from historical contingencies if the two plant subspecies have distinctive biogeographic origins and each imported their cohort of biotrophic fungi, as described for Carex curvula [START_REF] Geremia | Contrasting microbial biogeographical patterns between anthropogenic subalpine grasslands and natural alpine grasslands[END_REF].

The high abundance of putative pathogens beneath seemingly healthy cushion plants points to a reduced pathogenic virulence that may result from long-term associations. Pathogen virulence is predicted to decrease with increasing host lifespan [START_REF] Van Molken | Virulence in clonal plants: Conflicting selection pressures at work?[END_REF], and plant-fungi interactions may shift from negative to positive with elevation [START_REF] Defossez | Do interactions between plant and soil biota change with elevation? A study on Fagus sylvatica[END_REF]. The cost for the plant to host endophytic fungi has been shown to be weak [START_REF] Kia | Influence of phylogenetic conservatism and trait convergence on the interactions between fungal root endophytes and plants[END_REF]. Fungi associated with cushion plants may, for instance, offset their fitness costs because they reduce herbivore or virulent phytopathogen attacks or increase plant tolerance to abiotic stress [START_REF] Rodriguez | The Role of Fungal Symbioses in the Adaptation of Plants to High Stress Environments[END_REF]. In a situation where nutrients are scarce, there also might be an evolutionary trade-off between the cost of hosting a biotrophic fungus and the accelerated nutrient cycling by the same fungus in the litter. Such a trade-off has been predicted for insect grazers [START_REF] De Mazancourt | Grazing optimization, nutrient cycling, and spatial heterogeneity of plant-herbivore interactions: Should a palatable plant evolve?[END_REF].

Ecological implications

Our results suggest that the plant evolutionary dynamic of cushion plants strongly impacts the fungal composition of alpine ecosystems, possibly through diverging morphological traits. This extended phenotype sensu Whitham [START_REF] Whitham | Community and ecosystem genetics: A consequence of the extended phenotype[END_REF] may have implications for plant succession. Cushion plants, including S. acaulis, help secondary plant species to settle in the landscape [START_REF] Badano | Assessing impacts of ecosystem engineers on community organization: A general approach illustrated by effects of a high-Andean cushion plant[END_REF][START_REF] Bonanomi | Cushion plant morphology controls biogenic capability and facilitation effects of Silene acaulis along an elevation gradient[END_REF]. This settling may occur through reduction in physical stress [START_REF] Bertness | Positive interactions in communities[END_REF], but the presence of antagonist fungi may also be important in this process. The Janzen-Connell hypothesis (see [START_REF] Bever | Microbial Population and Community Dynamics on Plant Roots and Their Feedbacks on Plant Communities[END_REF] for a review) proposes that microorganisms specific to a plant species allow plant coexistence by reducing the fitness of their host and are, thus, at the origin of biodiversity in plant communities. In addition, it was shown that symbiotic root endophytes of cushion plants increase the fitness of the facilitated plants [START_REF] Molina-Montenegro | Fungal endophytes associated with roots of nurse cushion species have positive effects on native and invasive beneficiary plants in an alpine ecosystem[END_REF].

Plant-fungi interactions remain an overlooked component in the natural history and the present functioning of cold-adapted ecosystems as well as their responses to environmental changes.

Sampling along the entire ecological range of S. acaulis subspecies and other cushion plant species, together with the determination of fungi living in soil vs. plant tissues, may reveal the contributions of plant evolutionary dynamics and functional traits on fungal engineering.

ACK N OWLED G M ENTS

This paper is dedicated to Prof. Serge Aubert (1966Aubert ( -2015)), a passionate scientist who was an invaluable source of motivation for this work and who assisted with this project in many ways. Philippe Choler,

Marie

F

  I G U R E 1 Morphology variation in Silene acaulis and sampling design. (a) S. acaulis cushions vary in shape from loose (left) to intermediate (middle) to dense cushions (right). The scale bar unit is 1 cm. (b) Topographical and geological maps showing the sampling sites along the six elevation gradients. Only bedrock types where a plant population was sampled are displayed in the legend

  with minor modifications. Total DNA was extracted from leaf tissues using the DNeasy 96 Plant Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions, and plant DNA was digested for 2 hr at 37°C in a 20-µl mix using 2 U of Mse1 and 5 U of EcoRI (New England F I G U R E 2 Plant genetic groups and their geographic distribution. (a) Proportion of AFLP alleles on plant individuals along the sampling gradients. Colors and letters refer to genotypic clusters: "x," "l," and "a" for exscapa, longiscapa, and ambiguous plants, respectively, as assigned by STRUCTURE (K = 2). The gradient of origin, elevation, and morphology of the sampled plants (triangles for dense and circles for loose cushions) are indicated; letters "C" (calcareous) and "G" (granitic) indicate the bedrock type (except "Q", quartzitic). (b) NMDS plot of the Jaccard distances between AFLP fingerprints of cushion plants. Colors denote genotypic clusters, and shapes denote the bedrock of origin Biolabs, Ipswich, MA, USA). Double-stranded adaptors were ligated to the digested DNA in a 40 µl volume for 2 hr at 37°C using 1 U of T4 DNA Ligase (Roche). Products were then diluted to 1:10, and a preselective PCR (2 min at 72°C, 30 cycles of 30 s at 94°C, 30 s at 56°C, and 2 min at 72°C with a final elongation of 10 min at 72°C) was carried out in a 25 µl volume containing 3 µl of digested ligated products, 1X PCR Buffer II at pH 8.3, 1.5 mM of MgCl 2 , 80 µM of dNTP mix, 0.2 µM of each primer at 10 µM (EcoRI+1 ⁄ MseI+1), and 0.5 U of AmpliTaq Gold DNA Polymerase (Applied Biosystems, Forster City, CA, USA). After performing a 1:20 dilution of the pre-selective PCR products, a selective amplification was carried out in a 12.5 µl volume containing 2.5 µl of diluted pre-selective PCR product, 1X PCR Buffer II, 2.5 mM of MgCl 2 , 80 µM of dNTP mix, 0.2 µM of each FAM-labeled primer (EcoRI+3⁄ MseI+3, EcoRI+3⁄ MseI+3), 8 µg/ml BSA, and 0.5 U of AmpliTaq Gold DNA polymerase (10 min at 95°C; 13 cycles of 30 s at 94°C, 60 s at 65°C to 56°C, and 60 s at 72°C; 23 cycles of 30 s at 94°C, 60 s at 56°C, and 60 s at 72°C with a final elongation of 10 min at 72°C). These primer combinations were chosen based on preliminary tests according to[START_REF] Pompanon | Genotyping errors: Causes, consequences and solutions[END_REF], which resulted in clear banding patterns and sufficient variability. The PCR products of the EcoRI ⁄ MseI combinations were purified using columns of half 5% Sephadex G50 and half Sephacryl S200. The PCR products from each EcoRI ⁄ MseI primer pair were run separately for fragment length analysis; 1.5 µl of the FAM-labeled products was mixed with 10 µl of HiDi formamide and 0.1 µl of Genescan ROX 500 size standard (Applied Biosystems) and then electrophoresed on an ABI PRISM 3130 capillary sequencer (Applied Biosystems).

F

  I G U R E 3 Abiotic and biotic soil engineering. (a) Soil abiotic parameters ordinated by PCA. Arrows represent the correlations of chemical variables to axes magnified three times. (b) Dissimilarities in soil abiotic parameters between cognate bare and cushion soils. (c) Fungal community dissimilarities ordinated by NMDS; vectors indicate correlations of environmental variables with NMDS axes (all significant in a permutation test, p < 0.001). (d) Fungal communities' dissimilarities between cognate bare and cushion soils. Black diamonds indicate means, and distinct gray letters indicate means that differ significantly (Wilcoxon test, p < 0.001)

43

  Figure S5).

  both classes, MOTUs were nearly always recruited by one plant subspecies only, including on the common limestone bedrock. The contribution of cushion soil nutrients to fungal variations owes much to plant engineering, thus highlighting an indirect effect of plant genetics most likely to impact saprobes. Saprobes recurrently associated with F I G U R E 6 Composition of the fungal communities. Bare soils (n = 86) and cushion soils (n = 90) were compared. (a) Phylum abundances: Error bars show SE; letters differ when means are significantly different (p < 0.05). (b) Fungal family abundances. Ascomycota are in the red-to-yellow color range, Mucoromycota appear in blue, and Basidiomycota are in the green range. Families representing <1% of the total reads were pooled. MOTUs without assignment at the family level were pooled by phylum ("uf" prefix); "NAs" indicate MOTUs of unknown phylum. Families whose frequencies vary from bare soils to cushion soils bear a suffix ("_^", increase; "_v", decrease); the significance of the difference is indicated on sectors ("."; "*"; "**"; "***" for p < 0.1, 0.05, 0.01, and 0.001, respectively)

  -Noëlle Binet, and two anonymous referees also provided a number of insightful comments. This research was conducted at the Long-Term Socio-Ecological Research (LTSER) Zone Atelier Alpes, a member of the eLTER-Europe network, and at the Station Alpine J. Fourier, a member of AnaEE France supported by the project Investissement d'Avenir: ANR-11-INBS-0001AnaEE-Services. The Laboratoire d'Écologie Alpine is part of Labex OSUG@2020 (ANR10 LABX56). We are also grateful to the Ecrins National Park for authorizing the fieldwork. We acknowledge Julien Renaud and Charline Giguet-Covex for their help in preparing the geographic database and Figure 1, respectively; Eric Coissac and Frederic Boyer for providing the unpublished Robitools R package; and Laure Gallien, LucieZinger and Damien Georges for their help with data analysis. We thank the CNRS-LECA for providing funding to a group of young researchers to run the "Vertical Ecology" project. The ANR project METABAR (ANR 11 BSV7 020 01) contributed to this project as well, and J. Roy received a PhD grant from the French Research Ministry.


Michael Baumann
email: michael.baumann@imag.fr

Pierre Mozer

Vincent Daanen

Jocelyne Troccaz

Towards 3D ultrasound image based

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Computer-guidance for medical interventions on subcutaneous soft tissue targets is a challenging subject, since the target tracking problem is still not satisfactorily solved. The main difficulties are caused by the elasticity, mobility and inaccessibility of soft tissues. With 3D US a real-time volume imaging technology became available that provides enough spatial tissue information to make image-based tracking possible. Imagebased tracking is essentially a mono-modal image registration problem with a real-time constraint. The primary task is to find the physical transformation T in a transformation space T between two images of the same object. The choice of T depends on the underlying physical transformation (e.g. rigid, affine or elastic) and the requirements of the target application. An extensive review on registration methods is given in [START_REF] Zitova | Image registration methods: a survey[END_REF].

Nowadays, research on mono-modal 3D US registration of soft tissue images focusses on rapid deformation estimation. Most studies in this domain, however make the implicit assumption that the rigid part of the transformation to estimate is either small or known. Confronted with combinations of large rigid transformations and elastic deformations, the proposed solutions fail without rigid pre-registration. For many clinical applications large rigid transformations can be avoided by immobilizing both the patient and the US probe. In the case of interventions without total anesthesia this however causes considerable patient discomfort. Moreover, it is sometimes impossible to fix the US probe, e.g. when the probe serves as a guide for surgical instruments. The respiratory and the cardiac cycle can be additional sources of tissue displacements. In all these cases it is necessary to identify the rigid part of the transformation before carrying out image-based deformation estimation.

Estimation of large rigid transformations is basically a global optimization problem since common similarity measures exhibit search-friendly characteristics (e.g. convexity) only in a small region near the global solution. The computational burden of global optimization in a 6-D rigid transformation space is prohibitive for tracking tasks. [START_REF] Gu'eziec | Medical image registration using geometric hashing[END_REF][START_REF] Eadie | A real time pointer to a preoperative surgical planning index block of ultrasound images for image guided surgery[END_REF] propose to reduce the intra-interventional computation time of global searches by precomputing a feature-based index hash table. During intervention, similarity evaluation is replaced by computation of the geometric index followed by a fast data-base lookup. In the context of US image tracking, this approach has the disadvantage of relying on feature extraction, which often lacks robustness when confronted with partial target images, speckle and US shadows. Also, it cannot reduce the complexity of the optimization problem and pre-computation time is not negligible.

Relatively few investigations involving 3D US image based tracking of soft tissues have been reported. In the context of respiratory gated radiation treatment, [START_REF] Sawada | A technique for noninvasive respiratory gated radiation treatment based on a real time 3D ultrasound image correlation: A phantom study[END_REF] acquire a localized 3D US reference image of the liver or the pancreas in breath-hold state and register it rigidly with the treatment planning CT volume. During therapy, localized US slices of the organ are continuously compared with the reference volume using image correlation to retrieve the planning position of the organ. In [START_REF] Huang | Dynamic 3D ultrasound and MR image registration of the beating heart[END_REF] real-time 3D US images of the beating heart are registered multimodally with a set of 4-D MR images covering the entire cardiac cycle. A localizer is used to initialize the spatial registration process while the ECG signal serves for temporal alignment. The authors achieve precise rigid registration in an overall computation time of 1 second with a mutual information based rigid registration algorithm. In both studies relative rigid movements between probe and target organ are limited to movements caused by the respiratory or cardiac cycles, which are predictable and repeatable to a certain extent.

The target application of this work is 3D transrectal ultrasound (TRUS) prostate biopsy trajectory tracking. Today, prostate biopsies are carried out using 2D TRUS probes equipped with a guide for spring needle guns. With the current standard biopsy protocol, consisting typically of 12 regularly distributed samples, it is impossible to know the exact biopsy locations after acquisition, which makes precise biopsy-based tumor localization, quality control and targeted repeated biopsies impossible. A TRUSbased prostate tracking system would make it possible to project all sample locations into a reference image of the prostate and thus to identify the exact sampling locations.

Image-based prostate biopsy tracking is, however, challenging: (i) the gland moves and gets deformed under the pressure of the TRUS probe. (ii) The patient is neither immobilized nor under total anesthesia. Most patients move significantly during the biopsy procedure. (iii) Since the probe serves also to guide the rigidly attached needle, probe movements are important. Rotations around the principal probe axis of more than 180 • and tilting of up to 40 • are frequent. Also, the probe head wanders over the gland surface during needle placement, which leads to relative displacements of up to 3cm. The global search problem thus fully applies to prostate alignment: tracking a reference on a calibrated TRUS probe cannot solve the problem due to (i) and (ii), and it is not very success promising to minimize similarity measures on biopsy images using only fast down-hill optimizers because of (iii). In this study we propose a solution to the global search problem for TRUS prostate image tracking, which consists in a search space reduction using a probe movement model. We further identify an efficient intensity-based similarity measure for TRUS prostate images and describe a fast multi-resolution optimization framework. Finally, the robustness, accuracy, precision and performance of the method are evaluated on 237 prostate volumes from 14 patients.

Methods

A framework for US image-based tracking

The purpose of a tracking system is to provide the transformation between an object in reference space and the same object in tracking space at a given moment. In the case of image-based tracking, the reference space is determined by the choice of a reference image to which all subsequently acquired images will be registered. In the case of 3D TRUS prostate biospies, it is convenient to acquire a 3D US volume as reference just some minutes before the intervention. Unfortunately, most currently available 3D US systems do not provide real-time access to volume data. They can, however, visualize two or three orthogonal 2D (o2D) slices inside the field of view of the probe in real-time. These slices can be captured using a frame-grabber and used for registration with a previously acquired reference volume [START_REF] Sawada | A technique for noninvasive respiratory gated radiation treatment based on a real time 3D ultrasound image correlation: A phantom study[END_REF][START_REF] Huang | Dynamic 3D ultrasound and MR image registration of the beating heart[END_REF]. Note that compared to 2D US images, o2D planes deliver considerably more spatial information, which potentially makes 3D to o2D registration more robust than 3D to 2D registration. In this work we will evaluate both 3D to 3D and 3D to o2D registration for image-based tracking.

Registration algorithms can be separated into two main classes: intensity-based and feature-based algorithms. As it is challenging to define robust and fast feature extraction algorithms for US images of the prostate, due to the low SNR of US images and the absence of clearly identifiable geometric features in the prostate, this study focuses on intensity-based approaches. Intensity-based measures are known for their robustness in presence of noise and partial image overlaps [START_REF] Zitova | Image registration methods: a survey[END_REF].

Image registration can be modeled as a minimization process of an image similarity measure that depends on a transformation T . There exist robust and fast algorithms for local minimization of image similarity measures. The condition for convergency to the target transformation T is that the optimizer starts from a point inside the capture range of T [START_REF] Shekhar | Mutual information-based rigid and nonrigid registration of ultrasound volumes[END_REF]. However, the capture range of common intensity measures (e.g. the Pearson correlation coefficient (CC) or normalized mutual information (NMI)) is relatively small compared to the transformation space that can be observed for TRUS prostate biopsies. This problem can be attacked from two sides: the first approach is to extend the capture range by improving the similarity measure, while the second method consists in finding a point inside the capture range using a priori knowledge on the probe position.

Several parts of the registration approach require information about the prostate location in the reference image. For our purpose it is sufficient to set an axis-aligned bounding box on the prostate boundaries in the reference image. The bounding box has to be defined by the clinician. No bounding box is needed for the tracking images.

Extending the Capture Range

Similarity Measure: We chose CC as similarity measure since it yields a larger capture range than NMI for mono-modal US registration. Compared to sums of squared distances (SSD), it is insensitive to linear intensity transformations and is capable of detecting inverse correlations. Intensity shifts can occur due to probe pressure variation, while inverse correlations can be observed when evaluating transformations far from the physical solution, in particular for gradient magnitude images.

Multi-resolution pyramid: Optimizing on coarse resolution levels of a gaussian pyramid yields some important advantages: coarse levels are statistical aggregates of the original image which are free of high-frequency noise, in particular speckle noise. Once the optimization on the coarsest level is terminated, the solution will be refined on denser levels, but from a considerably better starting point. This approach not only improves the characteristics of the similarity measure by reducing noise, but also considerably speeds up registration time, as most of the optimization can be performed on low-resolution images.

Attribute-vector approach: The capture range can be extended by combining measures of different aspects of the images to be compared [START_REF] Shen | Hammer: hierarchical attribute matching mechanism for elastic registration[END_REF][START_REF] Foroughi | Elastic registration of 3D ultrasound images[END_REF]. Since there is a strong probability that the similarity measure produces for every aspect a significant minimum near the correct solution, it is possible to amplify and widen the capture range of the solution by combining the measures. Also, it is less likely that noise-related local minima are produced at identical locations, which makes it possible to flatten them out in a combined measure. For this study we chose to evaluate the image intensity and its gradient magnitude (I and J are the images to be compared):

En IJ (T ) := (1 -CC(I, J • T )) • (1 -CC(||∇I||, ||∇J • T ||)) (1)
To improve performance and since gradient intensities are highly random on noisy highresolution images, attribute vectors are only used on low resolution levels of the image pyramid.

Panorama images:

The pyramid-like form of the US beam and the fact that the probe also serves to guide the biopsy needle makes it unavoidable that the gland is often only partially imaged. Hence at least the reference image should contain the entire prostate; otherwise the similarity measure may yield random results when the image overlap gets too small during registration. We therefore acquire three partial prostate volumes using the following protocol: the operator first acquires one image where the prostate is centered in the US beam, and then takes two additional images with rotations of 60 • around the principal axis of the probe. Care is taken to avoid deformation and US shadows. The panorama image resulting from compounding these acquisitions finally serves as reference.

Finding a point in the capture range

Mechanical probe movement model: To estimate large transformations between images, it is necessary to find a point inside the capture range of the similarity measure. Regular sampling of a 6-D rigid transformation space using a very sparse grid size of 10 already requires 10 6 function evaluations, which results in an unacceptable computational burden. The physical constraints exerted by the rectum on probe movements, and the fact that the probe head always remains in contact with the thin rectal wall at the prostate location lead to the following assumptions: 1) the probe head is always in contact with the prostate membrane, 2) the most important rotations occur around the principal axis of the probe, and 3) all other rotations have a rotation point that can be approximated by a unique fixed point F P rect in the rectum.

With these assumptions it is possible to define a probe movement model based on a prostate surface approximation, the probe position in the US image (which is known) and a rotational fixed point in the rectum. As shown in Fig. 1(a), the prostate surface is approximated by a bounding-box aligned ellipsoid. The ellipsoid is modeled using a 2D polar parameterization P R Surf (α,β) . The origin P R Surf (0,0) of the parameterization corresponds to the intersection of the line from the prostate center C P ro to F P Rect . As illustrated in Fig. 1(b), P R Surf (α,β) implements assumption 1) by determining plausible US transducer positions on the prostate surface. Assumption 3) is satisfied by requiring that the principal probe axis must always pass through F P Rect . Finally, a rotation about the principal probe axis implements assumption 2) and thus adds a third DOF (See Fig. 1(c)).

Systematic Exploration

The 3D subspace defined by the probe movement model is systematically explored using equidistant steps. To minimize the computational burden, systematic exploration is performed on the coarsest resolution level. Since the exploration grid points do not change during an intervention, it is possible to precompute and to store all resclices of the panoramic image necessary for the evaluation of the intensity measure. The rotational space around the principal axis of the probe is unconstrained (360 • ), while tilting ranges are limited to the maximum value determined on test data, plus a security margin. The number of steps per dimension are also experimentally determined. The five best results of the systematic exploration are stored with the constraint that all transformations respect a minimum distance between each other. If two results are too close, only the best one is stored. Next, a local search using the Powell-Brent algorithm is performed only on the coarsest pyramid level for each of the five results. The best result of the five local searches is finally used as the start point for a multi-level local optimization. The last level of the final search can be chosen in function of the desired precision and computation time. Note that compared to a single multi-level local search, five local optimizations on the coarsest level are negligible in terms of computation time. The probe is then rotated and translated such that its US origin O US coincides with P R Surf (α, β). In (c), the probe is rotated around its principal axis by an angle λ.

Experiments and Results

The presented method was validated on 237 3D images of the prostate acquired during biopsy of 14 different patients. The imaging device was a GE 3D US Voluson 730 equipped with a volume-swept transrectal probe (GE RIC5-9). All images, except the images used for panorama image creation, were acquired immediately after a biopsy shot. Both 3D to 3D and 3D to o2D were evaluated. All registrations were carried out in a post-processing step. The o2D images used in the tests were not framegrabbed but reconstructed from 3D images. The image resolution was 200 3 . The voxel side lengths varied from 0.33mm to 0.47mm. A five-level resolution pyramid was used for 3D to 3D registration; for 3D to o2D only four levels were used. The final multilevel search was carried out from the coarsest to the third-finest level for 3D to 3D, and to the second-finest level for 3D to o2D registration. A total of 12960 grid points on the movement model were explored during a search run. Registration was carried out on a Pentium 4 with 3GHz.

To measure reproducibility and registration success, 10 registrations were carried out for each volume pair from slightly perturbated start points by adding noise of 2mm and 2 • . This yielded 10 transformations T i that approximate the unknown rigid transformation between the prostate in both volumes. The average transformation T of the T i was computed with the method presented in [START_REF] Gramkow | On averaging rotations[END_REF]. The euclidean distance error

ǫ i E = ||T i • C -T • C||,
with C being the image center, and the angular error ǫ i A , which corresponds to the rotation angle of T -1 i •T , were used to compute the root mean square (r.m.s.) errors ǫ E and ǫ A . A registration was considered successful if ǫ E < 2.0mm and ǫ A < 5 degrees, and if the result T was visually satisfactory when superimposing both volumes in a composite image (See Fig. 2(c)).

Reconstruction accuracy evaluation was more difficult to implement since there is no straight-forward gold standard. In some images, the needle trajectories from previous biopsies were still visible. In these cases, the trajectories were manually segmented, and the angular error between corresponding needle trajectories were used to evaluate rotational accuracy. Also, some patients had significant and clearly visible calcifications inside the prostate. The distances between segmented calcifications were used to determine the translational accuracy. Tab. 1 and Fig. 2 show the results of the evaluations.

3D-3D

3D-o2D The overhead introduced by the systematic model-based exploration accounts for about 25% of 3D-3D , and for 35% of 3D-o2D registration time. The five optimizations on the coarsest level account for about 10% in 3D-3D, and for 20% in 3D-o2D. Panorama image pre-processing and pre-computation of the images for systematic exploration are performed before the intervention and require about one minute of computation time. This study presents a fast and robust rigid registration framework for TRUS prostate images in the context of unconstrained patient movements, of only anatomy-constrained probe movements and of probe-induced prostate displacements. The algorithm yields reproducible results and acceptable accuracy for both 3D-3D and 3D-o2D registration.

The success-rate of 3D-3D registration is very satisfactory, since all failures were either due to significant US shadows caused by only partial contact of the probe head with the rectal wall or by air bubbles in the US contact gel, or to an insufficient US depth with the result that parts of the gland membrane are not visible in the images. In these cases the similarity measure fails because of missing information in the image, and an algorithmic remedy probably does not exist. Additional failures can be observed for 3D-o2D registration, in particular for very small prostates, for which the coronal plane does not contain any prostatic tissue. 3D-o2D registration is also more sensible to poor image quality (e.g. low contrast), to large deformations and to partial prostate images (for which often only one plane contains prostatic tissue). Note that the presented algorithm is not very sensible to bounding box placement precision.

Computation time of local searches could be accelerated using the GPU for image reslicing (which corresponds to approximatively 95% of the computational burden of a similarity measure evaluation), while further optimization of the systematic exploration would require parallelization of the evaluations.

The presented algorithm in particular accurately registers the prostate membranes that are distant to the probe head, and the urethra. The relatively high angular r.m.s. error observed in the needle reconstruction study can be explained with probe-related local deformations that are particularly strong at the needle entry point. We are currently working on a biomechanical gland deformation model that allows for estimation of deformations to improve the accuracy of tissue registration near the probe head.

Fig. 1 .

 1 Fig. 1. Mechanical probe movement model in 2D: (a) shows the computation of the search model surface origin P R Surf (0, 0) from the prostate center C P ro and the (hypothetical) rectal probe fixed point F P Rect . In (b), a 2D polar parameterization is used to determine a surface point P R Surf(α, β). The probe is then rotated and translated such that its US origin O US coincides with P R Surf (α, β). In (c), the probe is rotated around its principal axis by an angle λ.

Fig. 2 .

 2 Fig. 2. Registration accuracy: (a) shows the target image, and (b) the aligned panorama image. In (c) both volumes are superimposed to illustrate registration accuracy for the urethra (arrow), and (d) illustrates the registration accuracy in the upper gland.

Table 1 .

 1 Test results: Numbers in brackets indicate the number of evaluated registrations.

	Registration success	96.7% (237)	87.7% (237)
	Average computation time	6.5s (237)	2.3s (237)
	Angular precision ǫA (reproducibility, r.m.s.)	1.75 • (229)	1.71 • (208)
	Euclidean precision ǫE (reproducibility, r.m.s.) 0.62mm (229)	0.47mm (208)
	Needle trajectory reconstruction (r.m.s.)	4.72 • (10)	4.74 • (9)
	Needle trajectory reconstruction (max)	10.04 • (10)	10.5 • (9)
	Calcification reconstruction (r.m.s.)	1.41mm (189)	1.37mm (181)
	Calcification reconstruction (max)	3.84mm (189)	4.30mm (181)

Acknowledgements: This work was supported by grants from the Agence Nationale de la Recherche (TecSan program, SMI project), from the French Ministry of Industry (ANRT agency), from the French Ministry of Health (PHRC program, Prostate-echo project) and from Koelis S.A.S., France. The clinical data were acquired at the urology department of the Pitié la Salpétrière hospital, Paris.