# Desert truffle genomes reveal their reproductive modes and new insights into plant-fungal interaction and ectendomycorrhizal lifestyle

José Eduardo Marqués-gálvez, Shingo Miyauchi, Francesco Paolocci, Alfonso Navarro-ródenas, Francisco Arenas, Manuela Pérez-gilabert, Emmanuelle Morin, Lucas Auer, Kerrie Barry, Alan Kuo, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

José Eduardo Marqués-gálvez, Shingo Miyauchi, Francesco Paolocci, Alfonso Navarro-ródenas, Francisco Arenas, et al. Desert truffle genomes reveal their reproductive modes and new insights into plant-fungal interaction and ectendomycorrhizal lifestyle. New Phytologist, 2021, 229 (5), pp.2917-2932. 10.1111/nph.17044. hal-03082553

HAL Id: hal-03082553

https://hal.univ-lorraine.fr/hal-03082553

Submitted on 18 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Desert truffle genomes reveal their reproductive modes and new insights into plant-fungal interaction and ectendomycorrhizal lifestyle

José Eduardo Marqués-Gálvez<fnref n="1" />,²* ID, Shingo Miyauchi²* ID, Francesco Paolocci³ ID,
Alfonso Navarro-Ródenas<fnref n="1" /> ID, Francisco Arenas<fnref n="1" /> ID, Manuela Pérez-Gilabert⁴ ID, Emmanuelle Morin² ID,
Lucas Auer², Kerrie W. Barry⁵ ID, Alan Kuo⁵, Igor V. Grigoriev⁵,⁶ ID, Francis M. Martin²† ID,
Annegret Kohler²† ID and Asunción Morte<fnref n="1" />† ID

Author for correspondence:
Asunción Morte
Email: amorte@um.es

Received: 10 April 2020
Accepted: 9 October 2020

New Phytologist (2021) **229**: 2917–2932

doi: 10.1111/nph.17044

**Key words:** arid environment, desert truffles, drought stress, ectendomycorrhizal symbiosis, MAT genes, mycorrhiza, plant–mirobe interactions.

## Summary

* Desert truffles are edible hypogeous fungi forming ectendomycorrhizal symbiosis with plants of Cistaceae family. Knowledge about the reproductive modes of these fungi and the molecular mechanisms driving the ectendomycorrhizal interaction is lacking.

* Genomes of the highly appreciated edible desert truffles *Terfezia claveryi* Chatin and *Tirmania nivea* Trappe have been sequenced and compared with other Pezizomycetes. Transcriptomes of *T. claveryi* × *Helianthemum almeriense* mycorrhiza from well-watered and drought-stressed plants, when intracellular colonizations is promoted, were investigated.

* We have identified the fungal genes related to sexual reproduction in desert truffles and desert-truffles-specific genomic and secretomic features with respect to other Pezizomycetes, such as the expansion of a large set of gene families with unknown Pfam domains and a number of species or desert-truffle-specific small secreted proteins differentially regulated in symbiosis. A core set of plant genes, including carbohydrate, lipid-metabolism, and defence-related genes, differentially expressed in mycorrhiza under both conditions was found.

* Our results highlight the singularities of desert truffles with respect to other mycorrhizal fungi while providing a first glimpse on plant and fungal determinants involved in ecto to endo symbiotic switch that occurs in desert truffle under dry conditions.

## Introduction

The so-called desert truffles are a group of edible hypogeous fungi that establish mycorrhizal symbiosis with annual and perennial shrubs belonging to Cistaceae, a plant family adapted to arid and semiarid areas (Kovács & Trappe, 2014; Roth-Bejerano et al., 2014). These areas are characterized by an aridity index AI <0.5—that is, the ratio between the annual precipitation and potential evapotranspiration (United Nations Educational, Scientific and Cultural Organization, 1979)—and by poorly fertile soils with sandy texture and low inputs of organic matter (Bonifacio & Morte, 2014). Under these environmental conditions, mutualistic symbionts promote growth and survival of their host plants. For example, *Terfezia claveryi* mycorrhizas promote *Helianthemum almeriense* survival during drought periods by altering its physiological and

nutritional parameters (Morte et al., 2010), such as increased phos-
phorus, nitrogen (N), and potassium uptake (Morte et al., 2000).
As a result of this interaction, a fine-tune regulation of both plant
and fungal aquaporin expression (Navarro-Ródenas et al., 2013;
Marqués-Gálvez et al., 2020) and hydrogen peroxide content in
roots (Marqués-Gálvez et al., 2019) also occur. The type of mycor-
rhiza established by these fungi is known as an ectendomycorrhiza
(EEM). EEMs are characterized by the co-occurrence of an intercel-
lular Hartig net, intracellular hyphas penetrating the cortex cells,
where they form coil-like structures, and a thin and disordered fun-
gal mantle surrounding the colonized roots (Morte et al., 1994; Yu
et al., 2001). The occurrence of either intercellular or intracellular
mycorrhizal structures depends on several factors. In vitro condi-
tions, high auxin, high phosphate, and/or high water content favour
the intercellular mycorrhizal type, whereas field conditions, low
auxin, low phosphate, and/or low water availability favour the intra-
cellular mycorrhizal type (Gutiérrez et al., 2003; Zaretsky et al.,
2006a; Navarro-Ródenas et al., 2012, 2013).

<footnote n="1">¹Departamento de Biología Vegetal (Botánica), Facultad de Biología, Universidad de Murcia, Campus de Espinardo, Murcia 30100, Spain; ²INRAE, UMR 1136, Interactions Arbres/ Microorganismes (IAM), Centre INRAE GrandEst - Nancy, Université de Lorraine, Champenoux 54280, France; ³CNR-IBBR, Istituto di Bioscienze e Biorisorse, UOS di Perugia, Perugia 06128, Italy; ⁴Departamento de Bioquímica y Biología Molecular-A, Universidad de Murcia, Campus de Espinardo, Murcia 30100, Spain; ⁵US Department of Energy Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, CA 94598, USA; ⁶Department of Plant and Microbial Biology, University of California, Berkeley, Berkeley, CA 94598, USA</footnote>

<footnote>*These authors contributed equally as first authors.</footnote>

<footnote>†These authors contributed equally as senior authors.</footnote>
Among desert truffles, *Terfezia* and *Tirmania* ascocarps (fruit-
ing bodies) are well-known delicacies and largely marketed
mainly because of their taste and nutritional values (Kagan-Zur
*et al.*, 2014). *Terfezia claveryi* (Honrubia *et al.*, 2002; Morte
*et al.*, 2008) and *Terfezia boudieri* (Slama *et al.*, 2010; Kagan-Zur
*et al.*, 2014), both associating with *Helianthemum* spp. as host
plants, are the only ones that have been successfully cultivated,
the former in Spain and the latter in Tunisia and Israel. Since
most climate models point to increased temperature and lower
precipitation rates for the next decades, areas covered by arid and
semiarid ecosystems are expected to increase (Schlesinger *et al.*,
1990; Lavee *et al.*, 1998; Huang *et al.*, 2016). Thus, desert truf-
fles are becoming a promising, alternative crop in regions such as
the Mediterranean basin and the Middle East, where desertifica-
tion will cause an alarming increase in dry lands (Morte *et al.*,
2017). *Helianthemum almeriense* Pau is a drought deciduous
shrub that is well represented in the Mediterranean basin. This
species establishes EEM symbiosis with several desert truffles in
the wild, including *T. claveryi*. The symbiotic relationship
between *H. almeriense* and *T. claveryi* has been widely used as a
model for the research of desert truffle symbiosis at a basic and
applied level—see Morte *et al*. (2017) for review. However, at
the present time, the information on the molecular mechanics
involved in this particular symbiosis is scarce.

Several studies have revealed the specific genomic signature of
ectomycorrhizal, arbuscular mycorrhizal, or ericoid mycorrhizal
fungi (Martin *et al.*, 2008, 2010; Tisserant *et al.*, 2013; Kohler
*et al.*, 2015; Chen *et al.*, 2018; Martino *et al.*, 2018; Murat *et al.*,
2018; Morin *et al.*, 2019; Sun *et al.*, 2019; Venice *et al.*, 2020),
but similar studies have not yet been carried out on EEM fungi.
These genome and transcriptome analyses have identified hall-
marks defining the mycorrhizal lifestyle. Abundant and highly
upregulated specific effector-like small secreted proteins (SSPs)
together with convergent losses of secreted carbohydrate-active
enzymes (CAZymes) are two of the most recurrent traits associ-
ated with the mycorrhizal lifestyle (Kohler *et al.*, 2015). Although
complementary DNA-amplified fragment length polymorphism
analyses have been carried out to identify genes related to
*T. boudieri* × *Cistus incanus* symbiosis (Zaretsky *et al.*, 2006b),
and the genome of *T. boudieri* has been recently published
(Murat *et al.*, 2018), many key aspects of the biology of EEM
desert truffles remain to be elucidated. In addition, the impossi-
bility to mate these fungi under controlled conditions has pre-
vented mycologists shedding light into their reproductive modes
(Roth-Bejerano *et al.*, 2004). This problem could be overcome
by studying the structure and gene content of the mating type
(*MAT*) locus, the master regulator of sexual reproduction in
fungi. In the *MAT* locus of *Ascomycetes*, the *MAT1-1-1* and
*MAT1-2-1* genes, which regulate the expression of sex-specific
pheromone biosynthetic and signalling genes, are present
(Debuchy *et al.*, 2010). The two major mating systems in these
fungi are heterothallism and homothallism. In heterothallic
species, the isolates carry one of the two versions of the *MAT*
locus, with either *MAT1-1-1* or *MAT1-2-1* genes, and these two
versions are referred to as idiomorphs instead of alleles, because
overall DNA sequence homology between idiomorphs is absent

(Turgeon & Debuchy, 2007). In these fungi, only isolates that
differ at *MAT* are sexually compatible. Conversely, in homothal-
lic species, the isolates carry both *MAT1-1* and *MAT1-2* genes in
their genomes, and thus they are self-compatible (Bennett &
Turgeon, 2017). So far, *MAT* loci have not been identified in
any desert truffle, and therefore their mating mode (homothal-
lism vs heterothallism) is unknown (Murat *et al.*, 2018). As
shown for *Tuber melanosporum* (Rubini *et al.*, 2011, 2014; Linde
& Selmes, 2012; Zampieri *et al.*, 2012; Murat *et al.*, 2013; Le
Tacon *et al.*, 2014), *Tuber borchii* (Mello *et al.*, 2017; Leonardi
*et al.*, 2019), and *Tuber aestivum* (Molinier *et al.*, 2016; Splivallo
*et al.*, 2019), the characterization of the *MAT* locus has been criti-
cal for disclosing the reproduction mode of these fungi and, in
turn, for designing new cultivation strategies.

Here, we sequenced the genomes of *T. claveryi* strain T7 and *Tirmania nivea* strain G3 and compared them with those of other Pezizomycetes. We also assessed both plant and fungal transcriptomes of the symbiosis between *T. claveryi* and *H. almeriense* under well-watered (favouring intercellular colonization) and drought-stress (favouring intracellular colonization) conditions. We aimed to decipher the specific genomic and transcriptomic features of desert truffles compared with other Pezizaceae species. Our findings pave the way for a better understanding of sexual and vegetative propagation modes of the desert truffles and highlight singularities in their genomes and transcriptomes that may shape the ectendomycorrhizal interaction between these fungi and their host plants.

**Materials and methods**

Strains and fungal material for genome sequencing

*Terfezia claveryi* strain T7 was isolated from an ascocarp collected in an *H. almeriense × T. claveryi* experimental field located in Zarzadilla de Totana (Murcia, Spain) in April 2007. *Tirmania nivea* strain G3 was isolated from an ascocarp collected close to *Helianthemum lippi* plants in an arid natural site, located in Abu Dhabi (United Arab Emirates), in February 2014. Both strains have been maintained *in vitro* by subculturing every 3–6 months, using modified Melin–Norkrans (MMN) and/or optimal MMN (MMNo; Arenas *et al.*, 2018) liquid shaking culture (75 rpm) or solid media at 23°C in the dark.

Genome sequencing, assembly, and annotation

Genomic DNA and total RNA from axenic liquid culture of
*T. claveryi* and *T. nivea* mycelia were isolated by the cetyl
trimethylammonium bromide method (Chang *et al.*, 1993).
Genome sequencing details can be found in Supporting Informa-
tion Methods S1. They were annotated with the Joint Genome
Institute (JGI) Annotation Pipeline, which detects and masks
repeats and transposable elements (TEs), predicts genes, charac-
terizes each translated protein with sub-elements such as domains
and signal peptides, chooses a best gene model at each locus to
provide a filtered working set, clusters the filtered sets into draft
gene families, ascribes functional descriptions, and creates a JGI
MycoCosm genome portal with tools for public access and com-
munity-driven curation of the annotation (Grigoriev *et al.*, 2014;
Kuo *et al.*, 2014).

Phylogenetic analyses and comparative genomics

For phylogenetic tree construction and genome comparison,
3156 single-copy orthologues were selected among 13 published
fungal genomes belonging to the Pezizomycetes class (Table S1).
Single-copy orthologues are well suited for phylogenetic analyses
since they result from speciation events and go beyond the issues
of conventional markers (Debray et al., 2019). Orthologous
genes among the selected fungi were identified using FASTORTHO
program with the parameters set to 50% identity, 50% coverage,
inflation 3.0 (Wattam et al., 2014). Protein sequences were
downloaded from the JGI MycoCosm database (https://genome.
jgi.doe.gov/mycocosm/home). Clusters with single-copy genes
were identified and aligned with MAFFT 7.221 (Katoh & Stand-
ley, 2013), ambiguous regions eliminated, and single-gene align-
ments concatenated with GBLOCKS 0.91b (Cruickshank, 2000). A
phylogeny tree was constructed with RAxML 7.7.2 (Stamatakis,
2006) using the standard algorithm, the PROTGAMMAWAG
model of sequence evolution, 1000 bootstrap replicates, and ran-
dom seed number (12345).

Statistics of JGI genome assemblies (i.e. N50, number of
genes and scaffolds, genome size) were obtained from JGI
Mycocosm (https://genome.jgi.doe.gov/programs/fungi/index.jsf
). Genome completeness with single-copy orthologues was cal-
culated using BUSCO v.3.0.2 with default parameters (Sima
et al., 2015). The coverage of TEs in genomes was calculated
and visualized using our set of custom R scripts TINGO (Morin
et al., 2019). Secreted proteins were predicted as described pre-
viously by Pellegrin et al. (2015). The count and ratio of total
(occurrences in genomes) and secreted CAZymes, lipases, pro-
teases, and SSPs (< 300 amino acids) and the predicted plant
cell-wall-degrading enzymes (PCWDEs) and microbe cell-wall-
degrading enzymes (MCWDEs) were calculated, visualized,
and compared (note that some enzymes may belong to both
PCWDE and MCWDE categories). Global trends of ecologi-
cal groups were evaluated using nonmetric multidimensional
scaling with the count of total and secreted CAZymes using
metaMDS in R package VEGAN (Oksanen et al., 2016). The
output files generated were combined and visualized with the
set of custom R scripts PRINGO (Miyauchi et al., 2020a) incor-
porating R packages GGPLOT2, GGTREE and EGG (Wickham,
2009; Auguie, 2017; Yu et al., 2017).

Gene family evolution

ORTHOFINDER v.2.4.0 (Emms & Kelly, 2019) and CAFE v.4.2 (Han et al., 2013) were used for the computational analysis of gene family evolution. The proteomes from the 13 fungi (Table S1), plus *Cenoccum geophilum* and *Oidiodendron maius* as external group, were aligned to each other with BLASTP, and the orthogroups were inferred by the MCL algorithm (Enright et al., 2002) within the ORTHOFINDER pipeline.

The previously built phylogenetic tree was rooted and made
ultrametric with R package APE (Paradis & Schliep, 2019).
We estimated the species divergence times using the mcmctree
method implemented in PAML v.4.8. (Yang, 2007). The inde-
pendent-rates clock model, a WAG substitution model and
approximate likelihood calculation were used. The substitution
rates of each gene were estimated by CODEML under a global
clock model. We set the time unit to 100 Myr and applied
uniform priors on two calibration points with lower and
upper hard bounds: the Tuberaceae divergence time estima-
tion between 134 and 179 Ma (Bonito et al., 2013) and the
Pezizomycetes–Leotiomyceta divergence time estimation
between 400 and 487 Ma (Prieto & Wedin, 2013). The
orthogroups, together with the dated ultrametric tree, were
used by CAFE to estimate an error model accounting for
genome assembly errors (caerror command) and the birth–
death parameter λ. Such a model and λ value were finally
used to infer expansions/contractions for the protein families
analysed.

Pfam domain enrichment/depletion analysis

The top 100 Pfam domains (https://pfam.xfam.org/, v.32.0) identified in the gene repertoires of the 13 fungal species were selected. The occurrence frequency values were transformed into z-scores, which are a measure of relative enrichment and depletion. The hierarchical clustering was done with a Euclidian distance metric and average linkage clustering method. Data were visualized and clustered using MULTIEXPERIMENT VIEWER (http://mev tm4.org). The Mann–Whitney–Wilcoxon test was performed to compare each Pfam family between the four desert truffles (*T. claveryi*, *T. boudieri*, *T. nivea*, and *Kalaharituber pfeili*) and the rest of Pezizomycetes.

MAT and sexual-related genes identification

MAT 1-1-1 genes were identified by the JGI annotation pipeline
in the genomes of *T. nivea* and *K. pfeilii*. To search for all puta-
tive MAT genes in *T. claveryi*, *T. boudieri*, *T. nivea* and *K. pfeilii*
genomes, a local BLASTP was performed on the JGI website using
the *MAT 1-1-1* and *MAT 1-2-1* genes of filamentous
ascomycetes as queries and the default parameters. The genomic
regions surrounding the *MAT* locus in these species were then
scanned to search for conserved genes. Full coding sequences of
*MAT 1-1-1* and the HMG domain of *MAT 1-2-1* genes were
aligned using MAFFT 7.221 (Katoh & Standley, 2013), and maxi-
mum-likelihood analysis was done in IQ-TREE (Nguyen *et al.*,
2014) using an ultrafast bootstrap with 1000 replicates, random
seed number (567943) and the best-fit model (VT + I), chosen
according to Bayesian information criterion, determined with
MODELFINDER (Kalyaanamoorthy *et al.*, 2017). Phylogeny trees
were visualized using MEGA X (Kumar *et al.*, 2018). Genes related
to pheromone biosynthesis and signalling were also searched in
*T. claveryi* and *T. nivea* by BLASTP using orthologous genes
already identified in other Pezizomycetes as queries (Murat *et al.*,
2018).
Biological material and experimental conditions for
transcriptome analysis

Vegetative mycelium of strain T7 from *T. claverryi* was grown in Erlenmeyer flasks containing MMNo medium (Arenas *et al.*, 2018). The cultures were shaken at 75 rpm and maintained in the dark at 23°C for 3 months, then mycelium was collected, snap frozen in liquid N₂, and stored at -80°C until further analysis.

For host plant production, *H. almeriense* seeds were collected in Zarzadilla de Totana (Murcia, Spain) and germinated as described in Morte *et al.* (2008). Owing to the erratic growth and difficulties for *in vivo* mycorrhiza formation by *T. claverryi* mycelium (Arenas *et al.*, 2018), spore inoculation is the standard and preferred method for mycorrhizal plant production (Morte *et al.*, 2008, 2012, 2017). After 2 months of growth, 36 plants were transferred to 300 ml individual pots, where they were inoculated with approximately 10⁵ *T. claverryi* spore solution obtained from mature ascocarps from the same location stated earlier. Twelve plants were left uninoculated to act as nonmycorrhizal control plants (NMPs). Three months after spore inoculation, the mycorrhizal status of plants was assessed by using a stereomicroscope (Gutiérrez *et al.*, 2003), and 12 of the inoculated plants were subjected to drought stress. Soil moisture was maintained over 75% for the well-watered mycorrhizal plants (WWMPs) and NMPs, whereas it was lowered to 40% and maintained until shoot water potential was < -2 MPa in drought-stressed mycorrhizal plants (DSMPs). Shoot water potential Ψ was measured in the following way: 5 cm long plant apices were covered in dark during 1 h, cut, and immediately placed in a pressure chamber (Soil Moisture Equipment Corp., Santa Barbara, CA, USA) according to Scholander *et al.* (1965). Once DSMPs reached Ψ < -2 MPa, plant biomass was calculated and Chl content was estimated (Morte & Andrino, 2014). Secondary and tertiary roots containing apical tips from six plants per treatment were rinsed to remove soil and harvested. They were immediately stained with trypan blue and observed under an Olympus BH2 microscope (Tokyo, Japan), to calculate mycorrhizal colonization, mycorrhizal intensity, and intracellular colonization from six replicates per treatment. One hundred secondary and tertiary 2–5 mm root sections per sample were observed under the microscope and were classified as mycorrhizal or nonmycorrhizal depending on the presence/absence of *T. claverryi* mycorrhizal structures. The results are given as percentage of the total sections observed. To calculate the mycorrhiza formation intensity, each *T. claverryi* mycorrhizal structure observed was given a score of 1 (< 10% of mycorrhizal structures observed in one segment), 2 (< 50%), 3 (≥ 50%), or 4 (> 90%), according to the amount of mycorrhiza observed in a single section, similar to Derkowska *et al.* (2008). Then, the mycorrhiza formation intensity was calculated as the weighted average divided by 4. Each *T. claverryi* mycorrhizal colonization was also classified as intercellular or intracellular based on the presence of a well-differentiated intercellular Hartig net or intracellular hyphas detected inside root cells. To support these observations, another portion of roots was observed under a

stereomicroscope and some mycorrhizal root morphotypes were
selected and processed to make semi-thin sections of 0.5 µm,
using a Leica UC6 ultramicrotome (Vienna, Austria), then
stained with toluidine blue and observed under an Olympus
BH2 microscope. Secondary and tertiary roots, from three repli-
cates per treatment, were immediately snap frozen in liquid N₂
and stored at -80°C for further RNA extraction.

Transcriptome analysis of *Terfezia claverryi*

Total RNA was extracted from three biological replicates of *T. claverryi* free-living mycelium (FLM) and from mycorrhizal roots of WWMPs, DSMPs, and NMPs with the same methodology as reported earlier herein. RNA sequencing (RNA-seq) was performed using 2 × 75 bp Illumina HiSeq2500 sequencing after messenger RNA (mRNA) library construction by Lifesequencing SL (Valencia, Spain). Raw reads were analysed using FASTQC for quality (https://www.bioinformatics.babraham.ac.uk/projects/fa stqc/), paired, trimmed, and mapped against *T. claverryi* reference transcripts (https://mycocosm.jgi.doe.gov/Tercla1/Tercla1.home. html) using CLC Genomics WORKBENCH 11 (Qiagen, Courtaboeuf, France). For mapping, the minimum length fraction was 0.9, the minimum similarity fraction 0.9, and the maximum number of hits for a read was set to 10. Raw read counts were normalized (with the function *counts* in DESeq2), and differential transcription levels of the contigs (Bonferroni adjusted *P*<0.05) were calculated with DESeq2 (Love *et al.*, 2014). The RNA-seq data were consistent according to the quality assessment (Figs S1, S2), permitting our transcriptomic analyses. Condition-specific differentially expressed genes (log₂(fold-change) > 2 or < −2, false discovery rate FDR<0.05) were identified using a custom R script, SHIN + GO (Miyauchi *et al.*, 2016, 2017, 2018, 2020b). Symbiosis-induced and repressed genes were queried against the 13 genomes analysed using BLASTP, and they were classified according to eukaryotic orthologous groups and secretome categories (i.e. SSPs and secreted CAZymes) as described in Pellegrin *et al.* (2015).

De novo assembly of the plant transcriptome

A de novo transcriptome assembly for *H. almeriense* was per-
formed with sequencing reads from NMP, WWMP, and DSMP
samples using MEGAHIT v.1.1.3 (Li *et al.*, 2015). Transcripts from
each condition (NMP, WWMP, and DSMP) were mapped onto
the assembled contigs with BOWTIE2 v.2.3.0 (Langmead &
Salzberg, 2012), and reads mapped onto the contigs were
counted with SAMTOOLS v.1.7 (Li *et al.*, 2009). Only the contigs
with > 200 bases were selected for further steps, and then they
were annotated using BLAST on the nonredundant National Cen-
ter for Biotechnology Information (NCBI) sequence database
(March 2018 version) and JGI MycoCosm annotation data (July
2018 version) using DIAMOND v.0.9.19 (Buchfink *et al.*, 2014)
with the parameters more-sensitive max-target-seqs 1 max-hsps 1
evalue 0.00001. Taxonomic annotations for each contig were
selected based on the higher bit scores for the best hits obtained
from the NCBI or JGI databases.
The annotated contigs were further separated into three categories (plants, fungi, metazoa) according to their annotations, and unknown contigs were removed. Contigs showing sequence coverage > 80% with a known sequence on NCBI or JGI databases and with average reads > 5 combined from all conditions were selected for further analyses. A summary of the assembly and the contig selection process is shown in Table S2. Raw read counts were normalized (with the function *counts* in DESeq2) and differential transcription levels of the contigs (Bonferroni adjusted *P*<0.05) were calculated with DESeq2 (Love et al., 2014). Condition-specific differentially transcribed genes (log<sub>2</sub>(fold-change) > 2 or < -2, FDR < 0.05) were identified with the SHIN + GO pipeline (Miyauchi et al., 2016, 2017, 2018, 2020b).

## Results

### Desert truffles' genome features

As expected, the phylogenomic analyses placed both *T. claveryi* and *T. nivea* in the Pezizaceae family, within the Pezizales order, and close to the other sequenced desert truffles; that is, *T. boudieri* and *K. pfeilii* (Fig. S3). BUSCO score (*T. claveryi*, 83.1%; *T. nivea*, 86.5%) and RNA-seq read mapping (c. 96%, see Info page on JGI Genome portals) indicated that the assembly and annotation are of good quality. Genome size for *T. claveryi* (74.5 Mb) and *T. nivea* (63.1 Mb) is in the same range as other ectomycorrhizal Pezizaceae genomes analysed so far (70.0 ± 8.0 Mb), but significantly smaller than the genomes of ectomycorrhizal Tuberaceae (124.3 ± 19.8 Mb) (Table S3). *Terfezia claveryi* and *T. nivea* genome assemblies contain 12 765 and 10 591 genes, respectively, and the average gene content from Pezizaceae was not significantly different from Tuberaceae (Table S3). *Terfezia claveryi* and *T. nivea* assemblies contain 48.3% and 38.3% of TE sequences, respectively (Fig. S3), and this is significantly lower in Pezizaceae than in Tuberaceae (Table S3). The total and secreted CAZyme profiles of the four desert truffles are very similar, but they are different from Tuberaceae and non-mycorrhizal Pezizomycetes (Fig. S4). The number of total and secreted CAZymes is similar in the Pezizaceae desert truffles and the Tuberaceae truffles, but each family showed substantial differences in their CAZyme family distribution, such as fewer MCWDEs and a significantly lower number of genes targeting pectin and chitin in desert truffles (Fig. 1; Table S3).

### Desert truffle gene families

A total of 143 orthologous gene family expansions and 1474 contractions were inferred for the common ancestor of desert truffles. The number of orthologous gene families gained in desert truffle species relative to their most recent common ancestor ranges from 300 to 900 families. On the other hand, *K. pfeilii* presents a strikingly low number of orthologous gene family losses (nine) compared with the other desert truffles (Fig. 2a). Most of the expanded genes from desert truffles encode proteins with no known functions or functions related to protein-protein

interactions, signalling networks, stress, and/or cell death (Tables S4, S5).

In the hierarchical clustering of the abundance of the different Pfam protein domains annotated in the genomes of desert truffles and other Pezizomycetes, the desert truffles clustered together (Fig. S5). Only two gene families have an enriched number of genes, whereas 18 families have a lower gene content in desert truffles. Among the latter, it is worth noting the significant reduction in the gene copy number of several genes related to transport, like Major Facilitator Superfamily (MFS; 53.75 vs 77.56 on average, desert truffle vs nondesert truffles Pezizomycetes, respectively), ABC transporters (11.25 vs 14.88 on average, desert truffle vs nondesert truffles Pezizomycetes, respectively), or sugar transporters (20.50 vs 30.67 on average, desert truffle vs nondesert truffles Pezizomycetes, respectively) (Fig. 2b).

### Genes involved in sexual reproduction in desert truffles

Putative *MAT* genes were identified in *T. claveryi*, *T. nivea*, *T. boudieri* and *K. pfeilii*. The *MAT* 1-1-1 gene, coding for a transcription factor with an α-box domain, was located on scaffold 39 in *T. claveryi* assembly, scaffold 83 in *T. nivea*, and scaffold 10 in *K. pfeilii*, whereas the opposite mating type gene, *MAT* 1-2-1, was found on scaffold 37 in the *T. boudieri* assembly. Only in the *K. pfeilii* genome were both *MAT* genes found to also be present, on scaffold 14, the *MAT* 1-2-1 gene. This gene showed the typical HMG domain, maintained the intron position characteristic of *MAT* 1-2-1 genes (Fig. S6), and clustered together with Pezizaceae *MAT* 1-2-1 genes in a phylogenetic analysis (Fig. S7). The genomic regions flanking the *MAT* genes of all desert truffles showed synteny with other ascomycetes, highlighting a putative DNA lyase (APN) next to all desert truffle *MAT* loci except for *K. pfeilii* MAT 1-2 locus (Fig. 3, S8). Most of the key genes related to pheromone processing, perception, and signalling were also identified in *T. claveryi* and *T. nivea* as per BLASTP search (Table S6). Most of the genes related to the pheromone pathway were not differentially expressed between *T. claveryi* FLM and mycorrhizas, regardless of the fact that these structures were collected from WWMPs or DSMPs.

### Symbiosis-induced genes in *Terfezia claveryi*

When comparing the transcriptomic profiles of the mycorrhizal roots and FLM, 751 genes were found to be upregulated (5.89% of the total predicted genes, log<sub>2</sub>(fold-change) > 2, FDR < 0.05; Table S7), whereas 303 were found to be downregulated (2.37% of the total predicted genes, log<sub>2</sub>(fold-change) < -2, FDR < 0.05; Table S8). Among the highly upregulated genes, there were several genes that may be of importance for the establishment and development of functional mycorrhiza, such as SSPs, several transcription factors, N metabolism enzymes (e.g. nitrate assimilation cluster, cyanase), transporters (aquaporin, amino acid, or sugar transporters, among others), PCWDEs (expansins, cellulose-acting copper-dependent lytic polysaccharide monooxygenases (AA9), and carbohydrate esterase (CE8)),
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Secretomic profiles of 13 genomes. Left bubble plot shows the number of secreted genes for carbohydrate-active enzymes (CAZymes), lipases, proteases, and others (i.e. all secreted proteins not in these first three groups). SSP group refers to the number of small secreted proteins (SSPs; < 300 amino acids). The size of bubbles corresponds to the number of genes. Colours in the plots represent different taxonomic families. Left bar plots represent the ratio of CAZymes, lipases, and proteases to all secreted proteins (left) and the ratio of SSPs among the entire secretome (right). Right bubble plot shows the number of plant cell-wall-degrading enzymes (PCWDEs) and microbial cell-wall-degrading enzymes (MCWDEs; fungal cell walls and bacterial membranes), including lytic polysaccharide monooxygenase (LPMO), substrate-specific enzymes for cellulose, hemicellulose, lignin, and pectin (plant cell walls); chitin, glucan, mannan (fungal cell walls); peptidoglycan (bacterial membranes). Right bar plots show the total count of genes including PCWDEs and MCWDEs (left) and the ratio of PCWDEs to MCWDEs (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2** Gene family distribution in desert truffles. (a) Gene expansion and contraction along a time-calibrated phylogeny tree of 13 Pezizomycetes. Green and red numbers represent significantly (*P* < 0.01) expanded or contracted, respectively, gene family counts. (b) The heatmap depicts a double hierarchical clustering of the significant enriched or depleted (Mann-Whitney-Wilcoxon test, *P* < 0.05) Pfam families in the desert truffles Pezizaceae compared with the other nine Pezyzomicetes analysed. Numbers represent gene counts of each family. The top most frequent 100 domains analysed can be consulted in the Supporting Information Fig. S5.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Phylogenetic tree of *MAT* sequences from desert truffles and schematic organization of the mating-type loci. Phylogeny tree was constructed by maximum-likelihood method. Coloured arrows represent conserved genes among species. Blue, *MAT* 1-1-1; red, *MAT* 1-2-1; green, DNA lyase (APN); purple, Sam decarboxylase; orange, Cyclooxygenase 13 (COX13); yellow, mitochondrial ribosomal protein of the small subunit (RSM22); brown, homologue to GSTUMT2000009628001 from *Tuber melanosporum*; white, nonconserved coding regions. Full coding sequences from (a) *MAT* 1-1-1 and HMG domain from (b) *MAT* 1-2-1 proteins were aligned and phylogenetically analysed as explained in the Materials and Methods section.</figcaption></figure>

tryptophan metabolism, and auxin-related proteins (L-amino acid oxygenases, amidases, auxin efflux carrier; Table S7).

In *T. claveryi*, 33% of the symbiosis upregulated genes are conserved among the 13 fungi analysed (clusters VI and VII) (Fig. 4). In these clusters, CAZymes (e.g. AA9 and CE8), genes related to metabolism (e.g. oxidoreductases), cellular processes, and signalling (e.g. peptidases and heat shock proteins) are enriched in comparison with the whole gene repertoire (Fisher's exact test; Table S9). On the other hand, 20% of the upregulated genes are species specific and have no sequence similarity (cluster V). Clusters I to IV (47% of the upregulated genes) contain genes with sequence similarity with genes belonging to all or some of the desert truffles. SSPs are enriched in the whole set of upregulated genes and particularly in clusters I and II compared with the whole gene repertoire (Fisher's exact test; Table S9). Twelve out of the 15 SSPs shared by desert truffles (clusters I and II) also share a substantial level of genome synteny (Fig. S10).

## Fungal colonization and gene regulation under drought conditions

No significant differences were found between DSMPs and WWMPs in plant biomass or Chl content (Table S10). *Terfeza claveryi* colonization in DSMPs was greater and more intense (Table 1) than in WWMPs. A greater number of intracellular colonizations (up to 41.9%) was observed in DSMPs vs WWMPs (2.7%) (Fig. 5; Table 1). This morphological observation was corroborated by the higher number of *T. claveryi* RNA reads in DSMPs vs WWMPs *de novo* metatranscriptomes (see next subsection).

A total of 288 and 329 fungal genes (2.26% and 2.58%, respectively, of the entire gene set) were up or downregulated in DSMPs vs WWMPs, respectively (log₂(fold-change) > 2 or < -2,

respectively, FDR < 0.05; Tables S11, S12). Up to 26 fungal genes with a putative function in stress responses (catalase, thioredoxin, cupredoxins, superoxide dismutases, several heat shock proteins) were upregulated.

Under drought conditions, a significant downregulation of secreted CAZymes was found (Fisher's exact test, *P* = 0.0007167 for PCWDEs and *P* = 0.004683 for MCWDEs) compared with the whole gene repertoire. Seven PCWDEs targeting cellulose (AA9, GH131 and GH5_5), hemicellulose (GH43), and pectin (CE8, GH28 and PL4_1) and seven MCWDEs targeting chitin (AA11, CBM18, GH128, GH18 and GH132) and glucans (GH16 and GH17) were downregulated during drought stress conditions, whereas only a GH76 targeting mannan and a GH35 targeting hemicellulose were upregulated (Fig. 6).

## Plant co-assembly and gene regulation in mycorrhizal symbiosis

For the plant metatranscriptome analysis, 135 414 transcript contigs were assembled, and we selected 27 346 contigs based on stringent criteria for further transcriptomic analysis (see Materials and Methods section; Table S2). Within each treatment, > 80% of the reads were either from *T. claveryi* or the host plant with varied proportion of *T. claveryi* (Fig. S11): DSMPs was the treatment with more *T. claveryi* reads (21% vs 12% in WWMPs and 0.4% in NMPs; Table S13). The 0.4% of reads belonging to *T. claveryi* from NMPs can be attributed to slight nursery and/or laboratory contaminations, since no mycorrhizal structures were detected when observed under the optical microscope (Table 1). On average, 72% of the total contigs were from plant origin, as per BLAST results.

From the 19 696 plant contigs identified in *de novo* assembly, 163 were upregulated and 72 were downregulated in
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Sequence conservation and functional analysis of symbiosis-induced transcripts. The heatmap depicts a double hierarchical clustering of the symbiosis-upregulated *Terfezia claveryi* genes (rows, log₂(fold-change) > 2, FDR-corrected < 0.05; Supporting Information Table S7). Symbiosis-induced genes were queried against the genomes of the 13 genomes using BLASTP (see Fig. S9 for the heatmap of symbiosis-downregulated genes). Homologues are coloured from yellow to red depending on the percentage of similarity. Clusters of genes with coordinated expression are numbered and highlighted with different colours. Right of the heatmap, the percentages of putative functional categories are given for each cluster as bar plots. Clusters significantly enriched in small secreted proteins (SSPs) are marked with an asterisk (Fisher's exact test *P* < 0.05). CAZyme, carbohydrate-active enzymes; KOG, eukaryotic orthologous groups.</figcaption></figure>

<table><caption>**Table 1** Mycorrhizal colonization of *Terfezia claveryi* × *Helianthemum almeriense* plants under different treatments.</caption><thead><tr><th></th><th>Mycorrhizal colonization (%)</th><th>Mycorrhizal intensity (out of 4)</th><th>Intracellular colonization (%)</th></tr></thead><tbody><tr><td>Well-watered mycorrhizal plant</td><td>21.34 ± 4.83 a</td><td>1.52 ± 0.19 a</td><td>2.73 ± 1.88 a</td></tr><tr><td>Drought-stressed mycorrhizal plant</td><td>45.37 ± 6.62 b</td><td>2.20 ± 0.24 b</td><td>41.91 ± 9.98 b</td></tr><tr><td>Nonmycorrhizal plant</td><td>nd</td><td>nd</td><td>nd</td></tr></tbody></table>

Mean ± SE is presented for each variable. Data were submitted to Mann-Whitney-Wilcoxon test, and different letters represent significant differences (*P* < 0.05). nd, not detected.

WWMPs vs NMPs (log₂(fold-change) > 2 and < -2, respectively; Tables S14, S15), whereas 255 plant genes were upregulated and 611 downregulated (log₂(fold-change) > 2 and < -2, respectively; Tables S16, S17) in DSMPs vs WWMPs. Among the regulated genes during symbiosis in well-watered conditions, there were several ABC and MFS transporters, pathogen defence-related proteins (thaumatins, peroxidases), and oxidoreductases (e.g. multicopper oxidases, cytochrome P450). When comparing DSMPs with WWMPs, we found a majority of abiotic and biotic stress-response genes (e.g. heat shock proteins, spermidine synthases, thaumatin), transporters (e.g. sugar and amino acid transporters, MFS), carbohydrate-metabolism-related genes (mainly glycoside transferases and glycoside hydrolases), and several GH3 auxin-responsive genes were regulated.

In mycorrhizal plants vs NMPs, regardless of the host plant growing conditions, 20 and 30 genes were up and downregulated, respectively. They represent a core set of symbiosis-regulated *H. almeriense* genes (Fig. 7). Among the upregulated genes, those coding for proteins involved in carbohydrate and lipid metabolism and in transport were the most represented (30% of the total), whereas two genes coding for thaumatin-like pathogen-related proteins showed a further increase under stress. Conversely, among the set of downregulated genes, the majority were of unknown function, two coded for glycoside hydrolase, and a few showed a further decrease in DSMPs vs WWMPs. Among this subset of genes were one gene each for alternative oxidase (AOX) and pectinesterase and two for Pirins and protein kinases.

## Discussion

### Specific features of the desert truffle genomes

Although desert truffles share many of their genomic features with other Pezizomycetes species, they also exhibit several species-specific characteristics, as shown by our genomic analyses, coupled to the morpho-anatomical and transcriptomic analyses of *T. claveryi* × *H. almeriense* mycorrhizas. A large genome size and a high number of orphan SSPs, together with a reduced set of PCWDEs, are among the genomic traits shared by *T. claveryi* and *T. nivea* with other mycorrhizal fungi (Kohler *et al.*, 2015; Miyauchi *et al.*, 2020a). Conversely, the expansion of large sets of gene families with unknown Pfam domains are desert-truffles-specific traits, which may be linked to the specific processes that distinguish these fungi from other Pezizomycetes (e.g. different
**Fig. 5** *Terfezia claveryi* × *Helianthemum almeriense* mycorrhiza. (a) Club-shaped mycorrhizal root (CSMR), one of the four morphotypes described by Gutiérrez *et al.* (2003), with extraradical mycelium (EM). (b) Longitudinal section of a well-watered mycorrhiza (WWMP) with several hyphas forming a Hartig net (HN). (c) Cross-section of a WWMP with several hyphas forming an HN. (d) Cross-section of a drought-stressed mycorrhizal plant with majority of intracellular hyphas (IH) forming clumps.

<figure><img src="image_8.png" /></figure>

**Fig. 6** Regulation of *Terfezia claveryi* plant cell-wall-degrading enzymes and microbial cell-wall-degrading enzymes under drought stress in mycorrhizal symbiosis. Red, upregulation; blue, downregulation.

<figure><img src="image_9.png" /></figure>

host specificity, formation of EEMs, adaptation to arid climates). The significant reduction of the number of genes related to metabolite transport, such as MFS, ABC, or sugar transporters, is also a pattern that typifies both *Terfezia* spp. and *T. nivea* with respect to other Pezizomycetes. A similar trend was observed in the mutualistic endophyte *Piriformospora indica*, a well-adapted species to desert environments known to confer several benefits to its hosts under these conditions (Zuccaro *et al.*, 2011). These desert-truffle-specific signatures could be related to unknown mechanisms aimed to optimize water and nutrient resources, ultimately leading to adaptation to arid conditions.

## The reproductive modes in desert truffles

Genome sequencing is the only way to characterize the reproductive mode(s) of desert truffles, owing to our inability to master mating mechanisms under laboratory conditions. The occurrence of a single *MAT* gene (either *MAT1-1-1* or *MAT1-2-1*) in the genome of the sequenced strains of *T. claveryi*, *T. nivea* and *T. boudieri* is a signature that these fungi are heterothallic, whereas the presence of both *MAT* genes in *K. pfeilii* genome suggests that this species, differently from all the other symbiotic Pezizomycetes sequenced thus far (Murat *et al.*, 2018), might
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Core set of plant genes regulated in symbiosis. (a) Venn diagram representing the core upregulated (left) and downregulated (right) genes (log₂(fold-change) > 2 or < -2, respectively) and their eukaryotic orthologous groups classification. (b) The heatmap shows a double hierarchical clustering with the number of normalized transcripts by DESeq2 of core genes in three different conditions: nonmycorrhizal plant (NMP 1, 2 and 3), well-watered mycorrhiza (WWM 1, 4 and 6), and drought-stressed mycorrhiza (DSM 2, 3 and 9) replicates.</figcaption></figure>

self-fertilize. Future studies will be needed to verify whether both *K. pfeilii MAT* genes are functional and the conditions that promote their expression. The characterization of the structure and organization of the *MAT* locus and its surrounding genomic regions in the desert truffles provides crucial information to optimize the cultivation strategies for these fungi, model the genomic rearrangements within and around this locus, and trace the evolution of the breeding systems within Pezizomycetes.

## Genus and species-specific small secreted proteins

Mutual recognition of symbionts and development of mycorrhizal structures rely on a finely tuned exchange of signalling molecules between the partners, which are seen in the symbiosis of plants with Rhizobiaceae, arbuscular mycorrhiza, and ectomycorrhiza (Perret *et al.*, 2000; Bonfante & Genre, 2015; Martin, 2016). Recent studies showed that symbiosis-induced effector

proteins, such as *Laccaria bicolor* MiSSP7 and MiSSP7.6, interact with transcriptional regulators and dampen the local defence responses of *Populus* (Plett *et al.*, 2011, 2014; Kang *et al.*, 2020). Among the symbiosis upregulated genes of *T. claverryi*, we identified several mycorrhiza-induced SSPs (MiSSPs) showing sequence conservation and synteny within the desert truffles (*T. claverryi*, *T. boudieri*, *T. nivea* and *K. pfeilii*). Along with conserved MiSSPs, there are some proteins shared only between *T. claverryi* and its most closely related species (*T. boudieri*, *T. nivea*), and six MiSSPs are unique to *T. claverryi*. The three species *T. claverryi*, *T. boudieri* and *T. nivea* are known to form mycorrhizal symbiosis with the same host species (Moreno *et al.*, 2014). The presence of a common set of MiSSPs suggests the existence of conserved mechanism(s) among the desert truffles to interact with their hosts. At the same time, specific *T. claverryi* MiSSPs induced upon the symbiosis with *H. almeriense* suggest the possible existence of species-specific mechanisms. Further
functional analyses would be required to determine the involve-
ment of these effectors in EEM formation and their specificity
among different desert truffles and/or towards different host
species.

*Helianthemum almeriense* transcriptomic features related
to symbiosis

As the complete genome of *H. almeriense* is lacking, its tran-
scriptome was assembled *de novo* using our metatranscriptome
data set. This represents the first transcriptomics studies in Cis-
taceae, a family with a prominent role in the configuration of
Mediterranean xeric landscapes (Ellul *et al.*, 2002). We anno-
tated 19 696 transcript contigs, but a core set of only 50 genes
is strongly regulated in mycorrhizal conditions, regardless of
whether or not plants were under water stress. Upregulated
genes related to carbohydrate and lipid metabolism may reflect
the metabolic changes produced by the symbiotic state and the
exchange of nutrients between the plant and the fungus (Wang
et al., 2017). The upregulation of thaumatin-like genes, which
are known to play a role in pathogen defence (Liu et al.,
2010), could be either a direct defence response to limit
*T. claveryi* colonization or reflect an enhanced state of defence
against pathogens provoked by *T. claveryi* colonization, as it
happens in arbuscular mycorrhizal plants (Pozo et al., 2009).
Among the downregulated core genes, it is worth noting the
*alternative oxidase (AOX)* gene. The AOX pathway is an alter-
native to the cytochrome oxidase (COX) pathway for electron
transport in the respiration process, and it is uncoupled from
ATP production. AOX function is particularly important dur-
ing stress when it is involved in the optimization of respiratory
metabolism, which facilitates the continuous turnover of the
tricarboxylic acid and acts to prevent the excess generation of
reactive oxygen species (ROS) (Li et al., 2013; Vanlerberghe,
2013). The decrease of *AOX* mRNA levels in *T. claveryi* myc-
orrhizal roots concurs with previous studies that show suppres-
sion on AOX activity induced by arbuscular mycorrhizal
colonization, a condition in which the COX pathway seems to
be the predominant one (Liu et al., 2015; Del-Saz et al., 2018)
(see Notes S1 for additional discussion on this issue).

Intracellular colonization under dry conditions

*Terfezia claveryi* colonizes *H. almeriense* by differentiating symbi-
otic structures, which differ strikingly from ectomycorrhizal
structures. The extent of root colonization, but also the ability to
form either inter or intracellular structures, is driven by environ-
mental conditions. As previously observed (Navarro-Ródenas
et al., 2013; Marqués-Gálvez et al., 2019), drought stress
enhances the mycorrhizal colonization of *H. almeriense* roots by
*T. claveryi*. This morphological observation is supported by the
higher number of RNA reads corresponding to *T. claveryi* found
in dry conditions compared with well-watered conditions.
Intriguingly, we observed a shift in the mycorrhizal structures
formed under drought conditions. Approximately 40% of the
observed mycorrhizal roots entailed intracellular hyphal

structures, whereas *T. claveryi* does not generally penetrate into
host cells under well-watered conditions. Pectin-targeting
enzymes are thought to play a key role in the formation of the
intercellular Hartig net (Veneault-Fourrey *et al.*, 2014). We
observed that fungal genes coding for pectin targeting enzymes
(CE8, GH28 and PL4_1) were downregulated under the drought
condition, suggesting that the activity of these enzymes might be
dispensable when the intracellular colonization is favoured.
Whereas in ericoid mycorrhiza the intracellular penetration is
mainly attributed to PCWDEs because of their increased expres-
sion in mycorrhizal symbiosis (Martino *et al.*, 2018), the general-
ized downregulation of PCWDEs shown in this work suggests
that this may not be the case for EEM desert truffles.

The intracellular colonization of desert truffles under drought stress could depend more on the plant side than on the fungal symbiotic partner, as it has been suggested in arbuscular mycorrhizal associations (Tisserant *et al.*, 2013). Our transcriptomic data from *H. almeriense* show some CAZymes upregulated under drought stress. However, for most of them, the lack of full sequences makes it difficult to predict the signal peptides to infer their specific role in plant cell-wall remodelling. The role of plant hormones should also be taken into account. Auxin is already known to determine inter or intracellular colonization in *T. boudieri* × *C. incanus in vitro* (Zaretsky *et al.*, 2006a); and in the present study, certain auxin-related genes from both fungal and plant origins were upregulated under dry conditions. This hormone could act in enhancing lateral root development during stress (Korver *et al.*, 2018), thus increasing the surface available for fungal colonization and/or in loosening and remodelling plant cell walls and, possibly, facilitating the accommodation of the intracellular hyphas under drought stress, due to the promotion of ‘acid-growth’ (Majda & Robert, 2018).

Since intracellular fungal colonization appears to be driven by
water or nutritional stresses, it is difficult to infer whether differ-
entially expressed genes under drought stress are master regulators
of the switch in mycorrhizal types or are simply part of the myc-
orrhiza response to drought stress. Moreover, the lack of
drought-stressed nonmycorrhizal control plants in this experi-
ment also leaves open the possibility that differentially expressed
genes under drought might be a constitutive plant response to
drought. For instance, the massive upregulation of fungal stress-
response genes, in particular of several ROS scavengers, may have
a dual role. It could help in ROS alleviation and, hence, con-
tribute to mycorrhiza adaptation to dry conditions. At the same
time, it could serve as a mechanism to increase fungal coloniza-
tion by decreasing the oxidative burst that is normally produced
by plants to limit fungal expansion in the roots (Salzer *et al.*,
1999; Baptista *et al.*, 2007; Marqués-Gálvez *et al.*, 2019). Future
research aimed at decoupling the ecto/endo switch from environ-
mental stresses is required to shed more light on the specific
mechanisms regulating the host cell intracellular colonization by
ectendomycorrhizal fungi. The functional analyses of plant and
fungal genes highlighted in this work will also help to improve
our understanding of the interaction between desert truffles and
their host plants. Coupled to the knowledge of the mating sys-
tems of desert truffles gained through the present study, these
future analyses will be of pivotal importance to promote efficient
truffle farming in desiccated areas.

Acknowledgements

The project was funded by the JGI, a Department of Energy (DOE) Office of Science User Facility, and supported by the Office of Science of the US DOE under contract no. DE-AC02-05CH11231 within the framework of the Mycorrhizal Genomics Initiative (CSP no. 305), Metatranscriptomics of Forest Soil Ecosystems project (CSP no. 570), and the 1000 Fungal Genome projects (CSP no. 1974). This work has been partially financed by projects CGL2016-78946-R (AEI/FEDER, UE) and 20866/PI/18 (FEDER and Programa Regional de Fomento de la Investigación – Plan de Actuación 2019 – de la Fundación Séneca, Agencia de Ciencia y Tecnología of the Region of Murcia, Spain). This research was also supported by the Laboratory of Excellence ARBRE (ANR-11-LABX-0002-01), the Region Lorraine, the European Regional Development Fund, and the Plant-Microbe Interfaces Scientific Focus Area in the Genomic Science Program, and the Office of Biological and Environmental Research in the US DOE, Office of Science. JEMG was a beneficiary of a PhD grant (DI-14-06904) and ANR of a postdoctoral contract (IJCI-2016-28252), both from the Ministerio de Economía y Competitividad.

Author contributions

AM and FMM conceived, managed, and coordinated the project. AM, ANR, FMM and A Kohler planned and managed the research. IVG, KWB and A Kuo supervised the sequencing and annotation of the genomes. JEMG, SM, FP, FMM, AM, EM, LA, ANR, FA and A Kohler carried out the experimental set-up, bioinformatic analysis, and analysis and interpretation of the data. JEMG, SM and FP drafted and wrote the manuscript. MPG, ANR, FMM, A Kohler and AM contributed to sections of the manuscript. All the authors read and improved the manuscript.

ORCID

Francisco Arenas https://orcid.org/0000-0001-5850-1105

Kerrie W. Barry https://orcid.org/0000-0002-8999-6785

Igor V. Grigoriev https://orcid.org/0000-0002-3136-8903

Annegret Kohler https://orcid.org/0000-0002-9575-9567

José Eduardo Marqués-Gálvez https://orcid.org/0000-0003-1273-5082

Francis M. Martin https://orcid.org/0000-0002-4737-3715

Shingo Miyauchi https://orcid.org/0000-0002-0620-5547

Emmanuelle Morin https://orcid.org/0000-0002-7268-972X

Asunción Morte https://orcid.org/0000-0002-6426-0202

Alfonso Navarro-Ródenas https://orcid.org/0000-0001-9123-4547

Francesco Paolocci https://orcid.org/0000-0002-9394-876X

Manuela Pérez-Gilabert https://orcid.org/0000-0002-3706-5786

Data availability

Full genomes, predicted genes, and transcripts of *T. claveryi* T7 and *T. nivea* G3 can be found at the following JGI portals: https://genome.jgi.doe.gov/Tercla1/Tercla1.home.html and https://genome.jgi.doe.gov/Tirniv1/Tirniv1.home.html, respectively. RNA-seq data from *T. claveryi* × *H. almeriense* mycorrhiza and drought stress assay are available at Gene Expression Omnibus (GEO) with the accession nos. GSE154490 and GSE155042.

References

Arenas F, Navarro-Ródenas A, Chávez D, Gutiérrez A, Pérez-Gilabert M, Morte A. 2018. Mycelium of *Terfezia claveryi* as inoculum source to produce desert truffle mycorrhizal plants. *Mycorrhiza* **28**: 691–701.

Auguie B. 2017. *R package EGG*. [WWW document] URL https://cran.r-project.org/web/packages/egg/index.html [accessed 30 March 2019].

Baptista P, Martins A, Pais MS, Tavares RM, Lino-Nieto T. 2007. Involvement of reactive oxygen species during early stages of ectomycorrhiza establishment between *Castanea sativa* and *Pisolithus tinctorius*. *Mycorrhiza* **17**: 185–193.

Bennett RJ, Turgeon BG. 2017. Fungal sex: the Ascomycota. In: Heitman J, Howlett BJ, Crous PW, Stukenbrock EH, James TY, Gow NAR, eds. *The fungal kingdom*. Wiley Online Library: Hoboken, NJ, USA, 115–145.

Bonfante P, Genre A. 2015. Arbuscular mycorrhizal dialogues: do you speak “plantish” or “fungish”? *Trends in Plant Science* **20**: 150–154.

Bonifacio E, Morte A. 2014. Soil properties. In: Kagan-Zur V, Roth-Bejerano N, Sitrit Y, Morte A, eds. *Desert truffles. Phylogeny, physiology, distribution and domestication*. Soil Biology, vol. 38. Berlin, Germany: Springer-Verlag, 57–67.

Bonito G, Smith ME, Nowak M, Healy RA, Guevara G, Cazares E, Kinoshita A, Nouhra ER, Dominguez LS, Tedersoo L *et al*. 2013. Historical biogeography and diversification of truffles in the Tuberaceae and their newly identified Southern Hemisphere sister lineage. *PLoS ONE* **8**: e52765.

Buchfink B, Xie C, Huson DH. 2014. Fast and sensitive protein alignment using DIAMOND. *Nature Methods* **12**: 59–60.

Chang S, Puryear J, Cairney J. 1993. A simple and efficient method for isolating RNA from pine trees. *Plant Molecular Biology Reporter* **11**: 113–116.

Chen E, Morin E, Beaudet D, Noel J, Yildirim G, Ndikumana S, Charron P, St-Onge C, Giorgi J, Kruger M *et al*. 2018. High intraspecific genome diversity in the model arbuscular mycorrhizal symbiont *Rhizophagus irregularis*. *New Phytologist* **220**: 1161–1171.

Cruickshank R. 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. *Molecular Biology and Evolution* **17**: 540–552.

Debray K, Marie-Magdelaine J, Ruttink T, Clotault J, Foucher F, Malécot V. 2019. Identification and assessment of variable single-copy orthologous (SCO) nuclear loci for low-level phylogenomics: a case study in the genus *Rosa* (Rosaceae). *BMC Evolutionary Biology* **19**: e152.

Debuchy R, Berteaux-Lecellier V, Silar P. 2010. Mating systems and sexual morphogenesis in ascomycetes. In: Borkovich KA, Ebbole DJ, eds. *Cellular and molecular biology of filamentous fungi*. Washington, DC, USA: American Society of Microbiology, 501–535.

Del-Saz NF, Romero-Munar A, Cawthray GR, Palma F, Aroca R, Baraza E, Florez-Sarasa I, Lambers H, Ribas-Carbó M. 2018. Phosphorus concentration coordinates a respiratory bypass, synthesis and exudation of citrate, and the expression of high-affinity phosphorus transporters in *Solanum lycopersicum*. *Plant, Cell & Environment* **41**: 865–875.

Derkowska E, Sas-Paszt L, Sumorok B, Szwonek E, Gluszek S. 2008. The influence of mycorrhization and organic mulches on mycorrhizal frequency in apple and strawberry roots. *Journal of Fruit and Ornamental Plant Research* **16**: 227–242.

Ellul P, Boscaiu M, Vicente O, Moreno V, Rosello JA. 2002. Intra and interspecific variation in DNA content in *Cistus* (Cistaceae). *Annals of Botany* **90**: 345–351.
Emms DM, Kelly S. 2019. ORTHOFINDER: phylogenetic orthology inference for comparative genomics. *Genome Biology* 20: e238.

Enright AJ, Van Dongen S, Ouzounis CA. 2002. An efficient algorithm for large-scale detection of protein families. *Nucleic Acids Research* 30: 1575–1584.

Grigoriev IV, Nikitin R, Haridas S, Kuo A, Ohm R, Otillar R, Riley R, Salamov A, Zhao X, Korzeniewski F *et al*. 2014. MycoCosm portal: gearing up for 1000 fungal genomes. *Nucleic Acids Research* 42: 699–704.

Gutiérrez A, Morte A, Honrubia M. 2003. Morphological characterization of the mycorrhiza formed by *Helianthemum almeriense* Pau with *Terfezia claveryi* Chatin and *Picoa lefebvrei* (Pat.) Maire. *Mycorrhiza* 13: 299–307.

Han MV, Thomas GW, Lugo-Martinez J, Hahn MW. 2013. Estimating gene gain and loss rates in the presence of error in genome assembly and annotation using CAFE 3. *Molecular Biology and Evolution* 30: 1987–1997.

Honrubia M, Gutiérrez A, Morte A. 2002. Desert truffle plantation from south-east Spain. In: Hall IR, Wang Y, Danell E, Zambonelli A, eds. *Edible mycorrhizal mushrooms and their cultivation. Proceedings of the Second International Conference on Edible Mycorrhizal Mushrooms*. Christchurch, New Zealand, Crop & Food Research.

Huang J, Haipeng Y, Guan X, Wang G, Guo R. 2016. Accelerated dryland expansion under climate change. *Nature Climate Change* 6: 166–171.

Kagan-Zur V, Roth-Bejerano N, Sitrit Y, Morte A, editors. 2014. *Desert truffles. Phylogeny, physiology, distribution and domestication*. Soil Biology, vol. 38. Berlin/Heidelberg: Springer-Verlag. URL https://www.springer.com/gp/book/9783642400957.

Kalyaanamoorthy S, Minh BQ, Wong TK, von Haeseler A, Jermiin LS. 2017. MODELFINDER: fast model selection for accurate phylogenetic estimates. *Nature Methods* 14: 587–589.

Kang H, Chen X, Kempannien M, Pardo AG, Veneault-Fourrey C, Kohler A, Martin F. 2020. The small secreted effector protein MiSSP7.6 of *Laccaria bicolor* is required for the establishment of ectomycorrhizal symbiosis. *Environmental Microbiology* 22: 1435–1446.

Katoh K, Standley DM. 2013. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. *Molecular Biology and Evolution* 30: 772–780.

Kohler A, Kuo A, Nagy LG, Morin E, Barry KW, Buscot F, Canbäck B, Choi C, Cichoki N, Clum A *et al*. 2015. Convergent losses of decay mechanisms and rapid turnover of symbiosis genes in mycorrhizal mutualists. *Nature Genetics* 47: 410–415.

Korver RA, Koevoets IT, Testerink C. 2018. Out of shape during stress: a key role for auxin. *Trends in Plant Science* 23: 783–793.

Kovács G, Trappe J. 2014. Nomenclatural history and genealogies of desert truffles. In: Kagan-Zur V, Roth-Bejerano N, Sitrit Y, Morte A, eds. *Desert truffles. Phylogeny, physiology, distribution and domestication*. Soil Biology, vol. 38. Berlin, Germany: Springer-Verlag, 21–38.

Kumar S, Stecher G, Li M, Knyaz C, Tamura K. 2018. MEGA X: molecular evolutionary genetics analysis across computing platforms. *Molecular Biology and Evolution* 35: 1547–1549.

Kuo A, Bushnell B, Grigoriev IV. 2014. Fungal genomics: sequencing and annotation. In: Martin FM, ed. *Fungi. Advances in botanical research*, vol. 70. New York, NY, USA: Elsevier, 1–52.

Langmead B, Salzberg SL. 2012. Fast gapped-read alignment with BOWTIE 2. *Nature Methods* 9: 357–359.

Lavee H, Imeson AC, Sarah P. 1998. The impact of climate change on geomorphology and desertification along a Mediterranean–arid transect. *Land Degradation and Development* 9: 407–422.

Le Tacon F, Marçais B, Courvoisier M, Murat C, Montpied P, Becker M. 2014. Climatic variations explain annual fluctuations in French Périgord black truffle wholesale markets but do not explain the decrease in black truffle production over the last 48 years. *Mycorrhiza* 24: 115–125.

Leonardi P, Murat C, Puliga F, Iotti M, Zambonelli A. 2019. Ascoma genotyping and mating type analyses of mycorrhizas and soil mycelia of *Tuber borchii* in a truffle orchard established by mycelial inoculated plants. *Environmental Microbiology* 22: 964–975.

Li CR, Liang DD, Li J, Duan YB, Li H, Yang YC, Qin RY, Li L, Wei PC, Yang JB. 2013. Unravelling mitochondrial retrograde regulation in the abiotic stress induction of rice ALTERNATIVE OXIDASE 1 genes. *Plant, Cell & Environment* 36: 775–788.

Li D, Liu CM, Luo R, Sadakane K, Lam TW. 2015. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly by succinct de Bruijn graph. *Bioinformatics* 31: 1674–1676.

Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. 2009. The sequence alignment/map format and SAMTOOLS. *Bioinformatics* 25: 2078–2079.

Linde CC, Selmes H. 2012. Genetic diversity and mating type distribution of *Tuber melanosporum* and their significance to truffle cultivation in artificially planted truffières in Australia. *Application in Environmental Microbiology* 78: 6534–6539.

Liu JJ, Sturrock R, Ekramoddoullah AK. 2010. The superfamily of thaumatin-like proteins: its origin, evolution, and expression towards biological function. *Plant Cell Reports* 29: 419–436.

Liu Z, Li Y, Wang J, He X, Tian C. 2015. Different respiration metabolism between mycorrhizal and non-mycorrhizal rice under low-temperature stress: a cry for help from the host. *The Journal of Agricultural Science* 153: 602–614.

Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biology* 15: e950.

Majda M, Robert S. 2018. The role of auxin in cell wall expansion. *International Journal of Molecular Sciences* 19: e951.

Marqués-Gálvez JE, Morte A, Navarro-Ródenas A. 2020. Spring stomatal response to vapor pressure deficit as a marker for desert truffle fruiting. *Mycorrhiza* 30: 503–512.

Marqués-Gálvez JE, Morte A, Navarro-Ródenas A, García-Carmona F, Pérez-Gilabert M. 2019. Purification and characterization of *Terfezia claveryi* TcCAT-1, a desert truffle catalase upregulated in mycorrhizal symbiosis. *PLoS ONE* 14: e0219300.

Martin F, editor. 2016. *Molecular mycorrhizal symbiosis*. Hoboken, NJ, USA: John Wiley & Sons.

Martin F, Aerts A, Ahrén D, Brun A, Danchin EGJ, Duchaussoy F, Gibon J, Kohler A, Lindquist E, Pereda V *et al*. 2008. The genome of *Laccaria bicolor* provides insights into mycorrhizal symbiosis. *Nature* 452: 88–92.

Martin F, Kohler A, Murat C, Balestrini R, Coutinho PM, Jaillon O, Montanini B, Morin E, Noel B, Percudani R *et al*. 2010. Périgord black truffle genome uncovers evolutionary origins and mechanisms of symbiosis. *Nature* 464: 1033–1038.

Martino E, Morin E, Grelet GA, Kuo A, Kohler A, Daghino S, Barry KW, Cichoki N, Clum A, Dockter RB *et al*. 2018. Comparative genomics and transcriptomics depict ericoid mycorrhizal fungi as versatile saprotrophs and plant mutualists. *New Phytologist* 217: 1213–1229.

Mello A, Zampieri E, Zambonelli A. 2017. Truffle ecology: genetic diversity, soil interactions and functioning. In: Varma A, Prasad R, Tuteja N, eds. *Mycorrhiza – function, diversity, state of the art*. Cham, Switzerland: Springer, 231–252.

Miyauchi S, Hage H, Drula E, Lesage-Meessen L, Berrin J, Navarro D, Favel A, Chaduli D, Grisel S, Haon M *et al*. 2020b. Conserved white rot enzymatic mechanism for wood decay in the Basidiomycota genus *Pycnoporus*. *DNA Research* 27: dsaa011.

Miyauchi S, Kiss E, Kuo A, Drula E, Kohler A, Sánchez-García M, Andreopoulos B, Barry KW, Bonito G, Buée M *et al*. 2020a. Large-scale genome sequencing of mycorrhizal fungi provides insights into the early evolution of symbiotic traits. *Nature Communications* 11: e5125.

Miyauchi S, Navarro D, Grigoriev IV, Lipzen A, Riley R, Chevret D, Grisel S, Berrin JG, Henrissat B, Rosso MN. 2016. Visual comparative omics of fungi for plant biomass deconstruction. *Frontiers in Microbiology* 7: e1335.

Miyauchi S, Navarro D, Grisel S, Chevret D, Berrin JG, Rosso MN. 2017. The integrative omics of white-rot fungus *Pycnoporus coccineus* reveals co-regulated CAZymes for orchestrated lignocellulose breakdown. *PLoS ONE* 12: e0175528.

Miyauchi S, Rancon A, Drula E, Hage H, Chaduli D, Favel A, Grisel S, Henrissat B, Herpoël-Gimbert I, Ruiz-Dueñas FJ *et al*. 2018. Integrative visual omics of the white-rot fungus *Polyporus brumalis* exposes the biotechnological potential of its oxidative enzymes for delignifying raw plant biomass. *Biotechnology for Biofuels* 11: e201.

Molinier V, Murat C, Baltensweiler A, Büntgen U, Martin F, Meier B, Moser B, Sproll L, Stobbe U, Tegel W *et al*. 2016. Fine-scale genetic structure of natural *Tuber aestivum* sites in southern Germany. *Mycorrhiza* 26: 895–907.
Moreno G, Alvarado P, Manjón JL. 2014. Hypogeous desert fungi. In: Kagan-Zur V, Roth-Bejerano N, Sitrit Y, Morte A, eds. *Desert truffles. Phylogeny, physiology, distribution and domestication*. Soil Biology, vol. 38. Berlin Germany: Springer-Verlag, 3–20.

Morin E, Miyauchi S, Beaudet D, Hainaut M, Drula E *et al*. 2019. Comparative genomics of *Rhizophagus irregularis*, *R. cerebriforme*, *R. diaphanus* and *Gigaspora rosea* highlights specific genetic features in Glomeromycotina. *New Phytologist* 222: 1584–1598.

Morte A, Andrino A. 2014. Domestication: preparation of mycorrhizal seedlings. In: Kagan-Zur V, Roth-Bejerano N, Sitrit Y, Morte A, eds. *Desert truffles. Phylogeny, physiology, distribution and domestication*. Soil Biology, vol. 38. Berlin, Heidelberg: Springer-Verlag, 343–365.

Morte A, Andrino A, Honrubia M, Navarro-Ródenas A. 2012. *Terfezia* cultivation in arid and semiarid soils. In: Zambonelli A, Bonito GM, eds. *Edible ectomycorrhizal mushrooms*. Berlin, Germany: Springer-Verlag, 241–263.

Morte A, Cano A, Honrubia M, Torres P. 1994. *In vitro* mycorrhization of micropropagated *Helianthemum almeriense* plantlets with *Terfezia claveryi* (desert truffle). *Agricultural and Food Science* 3: 309–314.

Morte A, Honrubia M, Gutiérrez A. 2008. Biotechnology and cultivation of desert truffles. In: Varma A, ed. *Mycorrhiza*. Berlin, Germany: Springer-Verlag, 467–483.

Morte A, Lovisolo C, Schubert A. 2000. Effect of drought stress on growth and water relations of the mycorrhizal association *Helianthemum almeriense- Terfezia claveryi*. *Mycorrhiza* 10: 115–119.

Morte A, Navarro-Ródenas A, Nicolás E. 2010. Physiological parameters of desert truffle mycorrhizal *Helianthemum almeriense* plants cultivated in orchards under water deficit conditions. *Symbiosis* 52: 133–139.

Morte A, Pérez-Gilabert M, Gutiérrez A, Arenas F, Marqués-Gálvez JE, Bordallo JJ, Berná LM, Lozano-Carrillo C, Navarro-Ródenas A. 2017. Basic and applied research for desert truffle cultivation. In: Varma A, Prasad R. Tuteja N, eds. *Mycorrhiza – eco-physiology, secondary metabolites, nanomaterials*. Berlin, Germany: Springer-Verlag, 23–42.

Murat C, Payen T, Noel B, Kuo A, Morin E, Chen J, Kohler A, Krizsan K, Balestrini R, Da Silva C *et al*. 2018. Pezizomycetes genomes reveal the molecular basis of ectomycorrhizal truffle lifestyle. *Nature Ecology and Evolution* 2: 1956–1965.

Murat C, Rubini A, Riccioni C, De la Varga H, Akroume E, Belfiori B, Guarano M, Le Tacon F, Robin C, Halkett F *et al*. 2013. Fine-scale spatial genetic structure of the black truffle (*Tuber melanosporum*) investigated with neutral microsatellites and functional mating type genes. *New Phytologist* 199: 176–187.

Navarro-Ródenas A, Bárzana G, Nicolás E, Carra A, Schubert A, Morte A. 2013. Expression analysis of aquaporins from desert truffle mycorrhizal symbiosis reveals a fine-tuned regulation under drought. *Molecular Plant–Microbe Interactions* 26: 1068–1078.

Navarro-Ródenas A, Pérez-Gilabert M, Torrente P, Morte A. 2012. The role of phosphorus in the *ectendomycorrhiza continuum* of desert truffle mycorrhizal plants. *Mycorrhiza* 22: 565–575.

Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ. 2014. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. *Molecular Biology and Evolution* 32: 268–274.

Oksanen O, Blanchet FG, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P, Henry M *et al*. 2016. *VEGAN package v.2.5*. [WWW document] URL https://cran.r-project.org/web/packages/vegan/index.html.

Paradis E, Schliep K. 2019. *APE* 5.0: an environment for modern phylogenetics and evolutionary analyses in *R*. *Bioinformatics* 35: 526–528.

Pellegrin C, Morin E, Martin F, Veneault-Fourrey C. 2015. Comparative analysis of secretomes from ectomycorrhizal fungi with an emphasis on small-secreted proteins. *Frontiers in Microbiology* 6: e1278.

Perret X, Staehelin C, Broughton WJ. 2000. Molecular basis of symbiotic promiscuity. *Microbiology and Molecular Biology Reviews* 64: 180–201.

Plett JM, Daguerre Y, Wittulsky S, Vayssieres A, Deveau A, Melton SJ, Kohler A, Morrel-Falvey JL, Brun A, Veneault-Fourrey C *et al*. 2014. Effector MiSSP7 of the mutualistic fungus *Laccaria bicolor* stabilizes the *Populus* JAZ6 protein and represses jasmonic acid (JA) responsive genes. *Proceedings of the National Academy of Sciences, USA* 111: 8299–8304.

Plett JM, Kemppainen M, Kale SD, Kohler A, Legué V, Brun A, Tyler BM, Pardo AG, Martin F. 2011. A secreted effector protein of *Laccaria bicolor* is required for symbiosis development. *Current Biology* 21: 1197–1203.

Pozo MJ, Verhage A, García-Andrade J, García JM, Azcón-Aguilar C. 2009. Priming plant defence against pathogens by arbuscular mycorrhizal fungi. In: Azcón-Aguilar C, Barea JM, Gianinazzi S, Gianinazzi-Pearson V, eds. *Mycorrhizas: functional processes and ecological impact*. Berlin, Germany: Springer, 123–135.

Prieto M, Wedin M. 2013. Dating the diversification of the major lineages of Ascomycota (Fungi). *PLoS ONE* 8: e65576.

Roth-Bejerano N, Li YF, Kagan-Zur V. 2004. Homokaryotic and heterokaryotic hyphae in *Terfezia*. *Antonie van Leeuwenhoek* 85: 165–68.

Roth-Bejerano N, Navarro-Ródenas A, Gutiérrez A. 2014. Types of mycorrhizal association. In: Kagan-Zur V, Roth-Bejerano N, Sitrit Y, Morte A, eds. *Desert truffles. Phylogeny, physiology, distribution and domestication*. Soil Biology, vol. 38. Berlin, Germany: Springer-Verlag, 69–80.

Rubini A, Belfiori B, Riccioni C, Tisserant E, Arcioni S, Martin F, Paolocci F. 2011. Isolation and characterization of *MAT* genes in the symbiotic ascomycete *Tuber melanosporum*. *New Phytologist* 189: 710–722.

Rubini A, Riccioni C, Belfiori B, Paolocci F. 2014. Impact of the competition between mating types on the cultivation of *Tuber melanosporum*: Romeo and Juliet and the matter of space and time. *Mycorrhiza* 24: 19–27.

Salzer P, Corbière H, Boller T. 1999. Hydrogen peroxide accumulation in *Medicago truncatula* roots colonized by the arbuscular mycorrhiza-forming fungus *Glomus intraradices*. *Planta* 208: 319–325.

Schlesinger WH, Reynolds JF, Cunningham GL, Laura F, Jarrell WM, Virginia RA, Whitford WG. 1990. Biological feedbacks in global desertification conceptual models for desertification. *Science* 247: 1043–1048.

Scholander PF, Bradstreet ED, Hemmingsen EA, Hammel HT. 1965. Sap pressure in vascular plants. *Science* 148: 339–346.

Sima FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. 2015. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. *Bioinformatics* 31: 3210–3212.

Slama A, Fortas Z, Boudabous A, Neffati M. 2010. Cultivation of an edible desert truffle (*Terfezia boudieri Chatin*). *African Journal of Microbiology Research* 4: 2350–2356.

Splivallo R, Vahdatzadeh M, Maciá-Vicente JG, Molinier V, Peter M, Egli S, Uroz S, Paolocci F, Deveau A. 2019. Orchard conditions and fruiting body characteristics drive the microbiome of the black truffle *Tuber aestivum*. *Frontiers in Microbiology* 10: e1437.

Stamatakis A. 2006. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. *Bioinformatics* 22: 2688–2690.

Sun X, Chen W, Ivanov S, MacLean AM, Wight H, Rajmaraj T, Mudge J, Harrison MJ, Fei Z. 2019. Genome and evolution of the arbuscular mycorrhizal fungus *Diversispora epigaea* (formerly *Glomus versiforme*) and its bacterial endosymbionts. *New Phytologist* 221: 1556–1573.

Tisserant E, Malbreil M, Kuo A, Kohler A, Symeonidi A, Balestrini R, Charron P, Duensing N, Frei dit Frey N, Gianinazzi-Pearson V, *et al*. 2013. Genome of an arbuscular mycorrhizal fungus provides insight into the oldest plant symbiosis. *Proceedings of the National Academy of Sciences, USA* 110: 20117–20122.

Turgeon BG, Debuchy R. 2007. *Cochliobolus* and *Podospora*: mechanisms of sex determination and the evolution of reproductive lifestyle. In: Heitman J, Kronstad JW, Taylor JW, Casselton LA, eds. *Sex in fungi: molecular determination and evolutionary implications*. Hoboken, NJ, USA: John Wiley & Sons, 91–121.

United Nations Educational, Scientific and Cultural Organization. 1979. *Map of the world distribution of arid regions: map at scale 1:25,000,000 with explanatory note*. MAB Technical Notes 7. Paris, France: UNESCO.

Vanlerberghe GC. 2013. Alternative oxidase: a mitochondrial respiratory pathway to maintain metabolic and signaling homeostasis during abiotic and biotic stress in plants. *International Journal of Molecular Sciences* 14: 6805–6847.
Veneault-Fourrey C, Commun C, Kohler A, Morin E, Balestrini R, Plett J, Danchin E, Coutinho P, Wiebenga A, de Vries RP *et al.* 2014. Genomic and transcriptomic analysis of *Laccaria bicolor* CAZome reveals insights into polysaccharides remodelling during symbiosis establishment. *Fungal Genetics and Biology* 72: 168–181.

Venice F, Ghignone S, Salvioli di Fossalunga A, Amselem J, Novero M, Xianan X, Toro KS, Morin E, Lipzen A, Grigoriev IV *et al.* 2020. At the nexus of three kingdoms: the genome of the mycorrhizal fungus *Gigaspora margarita* provides insights into plant, endobacterial and fungal interactions. *Environmental Microbiology* 22: 122–141.

Wang W, Shi J, Xie Q, Jiang Y, Yu N, Wang E. 2017. Nutrient exchange and regulation in arbuscular mycorrhizal symbiosis. *Molecular Plant* 10: 1147–1158.

Wattam AR, Abraham D, Dalay O, Disz TL, Driscoll T, Gabbard JL, Gillespie JJ, Gough R, Hix D, Kenyon R *et al.* 2014. PATRIC, the bacterial bioinformatics database and analysis resource. *Nucleic Acids Research* 42: 581–591.

Wickham H. 2009. *GGPLOT2: elegant graphics for data analysis*. New York, NY, USA: Springer.

Yang Z. 2007. *PAML* 4: phylogenetic analysis by maximum likelihood. *Molecular Biology and Evolution* 24: 1586–1591.

Yu G, Smith DK, Zhu H, Guan Y, Lam TTY. 2017. *GGTREE: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data*. *Methods in Ecology and Evolution* 8: 28–36.

Yu TE, Egger KN, Peterson LR. 2001. Ectendomycorrhizal associations – characteristics and functions. *Mycorrhiza* 11: 167–177.

Zampieri E, Rizzello R, Bonfante P, Mello A. 2012. The detection of mating type genes of *Tuber melanosporum* in productive and non productive soils. *Applied Soil Ecology* 57: 9–15.

Zaretsky M, Kagan-Zur V, Mills D, Roth-Bejerano N. 2006a. Analysis of mycorrhizal associations formed by *Cistus incanus* transformed root clones with *Terfezia boudieri* isolates. *Plant Cell Reports* 25: 62–70.

Zaretsky M, Sitrit Y, Mills D, Roth-Bejerano N, Kagan-Zur V. 2006b. Differential expression of fungal genes at preinfection and mycorrhiza establishment between *Terfezia boudieri* isolates and *Cistus incanus* hairy root clones. *New Phytologist* 171: 837–846.

Zuccaro A, Lahrmann U, Güldener U, Langen G, Pfiffi S, Biedenkopf D, Wong P, Samans B, Grimm C, Basiewicz M *et al.* 2011. Endophytic life strategies decoded by genome and transcriptome analyses of the mutualistic root symbiont *Piriformospora indica*. *PLoS Pathogens* 7: e100229.

Supporting Information

Additional Supporting Information may be found online in the Supporting Information section at the end of the article.

**Fig. S1** *Terfezia claveryi* distribution and density of normalized log<sub>2</sub> transformed read counts of 9247 genes from nine replicates, three per treatment.

**Fig. S2** Correlation of *Terfezia claveryi* transcriptomes among nine replicates, three per treatment.

**Fig. S3** Phylogenomic tree with general features and statistics.

**Fig. S4** CAZyme-based clusters of 13 fungi.

**Fig. S5** Presence and abundance of different Pfam domain-containing proteins in 13 Pezizomycetes.

**Fig. S6** Alignment of the MAT 1-2-1 proteins of *K. pfeilii* with other Pezizomycotina species.

**Fig. S7** Phylogenetic tree of putative HMG domain containing proteins from *Kalaharituber pfeilii* and *Tuber melanosporum*, together with known MAT 1-2-1 proteins from *Tuberaceae* and *Pezizaceae*.

**Fig. S8** Schematic representation of the second *MAT* locus identified only in *Kalaharituber pfeilii* and the genome synteny with other desert truffles in the regions *framing* the *MAT* locus.

**Fig. S9** Sequence conservation and functional analysis of symbiosis-downregulated transcripts.

**Fig. S10** Genome synteny of SSPs shared by desert truffles.

**Fig. S11** Kingdom contig distribution from reads of the meta-transcriptome of *T. claveryi* × *H. almeriense* roots.

**Methods S1** *Terfezia claveryi* and *Tirmania nivea* genome sequencing.

**Notes S1** Additional discussion on *Helianthemum almeriense* alternative oxidase (AOX) genes downregulation in mycorrhizal plants under drought stress.

**Table S1** Genomic features of *Tefezia claveryi* and *Tirmania nivea* together with the rest of the fungi used for genome comparison.

**Table S2** Summary of the number of contigs during the selection process for de novo plant, fungi and metazoa assembly.

**Table S3** Differences in genomic and secretomic features between *Tuberaceae* and *Pezizaceae* mycorrhizal members used in this study.

**Table S4** Most significantly expanded orthogroup gene families according to CAFE analysis for *Terfezia claveryi* (Tercla, green), *Terfezia boudieri* (Terbo, red), *Tirmania nivea* (Tirniv, blue) and *Kalaharituber pfeilii* (Kalpfe, yellow).

**Table S5** Most significantly contracted orthogroup gene families according to CAFE analysis for *Terfezia claveryi* (Tercla, green), *Terfezia boudieri* (Terbo, red), *Tirmania nivea* (Tirniv, blue) and *Kalaharituber pfeilii* (Kalpfe, yellow).

**Table S6** Genes involved in pheromone biosynthesis and signalling in desert truffles, and its expression in *T. claveryi*.

**Table S7** Upregulated *T. claveryi* genes (log<sub>2</sub> fold-change > −2) in mycorrhizal roots vs FLM.

**Table S8** Downregulated *T. claveryi* genes (log<sub>2</sub> fold-change > −2) in mycorrhizal roots vs FLM.
**Table S9** Fisher's exact test summary results for *Terfezia claveryi* mycorrhiza differentially expressed genes clusters.

**Table S10** Plant growth under different water conditions.

**Table S11** Upregulated (log₂ fold-change > 2) *T. claveryi* genes comparing mycorrhizal roots drought stressed with control well-watered mycorrhizal roots.

**Table S12** Downregulated (log₂ fold-change < -2) *T. claveryi* genes comparing mycorrhizal roots drought stressed with control well-watered mycorrhizal roots.

**Table S13** Fungal reads present in plant de novo co-assembly.

**Table S14** Upregulated (log₂ fold-change > 2) *H. almeriense* genes comparing mycorrhizal roots with non mycorrhizal roots.

**Table S15** Downregulated (log₂ fold-change < -2) *H. almeriense* genes comparing mycorrhizal roots with non mycorrhizal roots.

**Table S16** Upregulated (log₂ fold-change > 2) *H. almeriense* genes comparing mycorrhizal roots drought stressed with control well-watered mycorrhizal.

**Table S17** Downregulated (log₂ fold-change < -2) *H. almeriense* genes comparing mycorrhizal roots drought stressed with control well-watered mycorrhizal roots.

Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the *New Phytologist* Central Office.