# Identification of quantitative trait loci (QTLs) for key cheese making phenotypes in the blue-cheese mold Penicillium roqueforti

Thibault Caron, Ewen Crequer, Mélanie Le Piver, Stéphanie Le Prieur, Sammy Brunel, Alodie Snirc, Gwennina Cueff, Daniel Roueyre, Michel Place, Christophe Chassard, et al.

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

## To cite this version:

Thibault Caron, Ewen Crequer, Mélanie Le Piver, Stéphanie Le Prieur, Sammy Brunel, et al. Identification of quantitative trait loci (QTLs) for key cheese making phenotypes in the blue-cheese mold Penicillium roqueforti. PLoS Genetics, 2025, 21 (4), pp.e1011669. 10.1371/journal.pgen.1011669. hal-05072919

HAL Id: hal-05072919

https://hal.inrae.fr/hal-05072919v1

Submitted on 19 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
RESEARCH ARTICLE

# Identification of quantitative trait loci (QTLs) for key cheese making phenotypes in the blue-cheese mold Penicillium roqueforti

Thibault Caron<sup>1,2,*</sup>, Ewen Crequer<sup>1,3*</sup>, Mélanie Le Piver<sup>2</sup>, Stéphanie Le Prieur<sup>1</sup>, Sammy Brunel<sup>2</sup>, Alodie Snirc<sup>1</sup>, Gwennina Cueff<sup>3</sup>, Daniel Roueyre<sup>2</sup>, Michel Place<sup>2</sup>, Christophe Chassard<sup>4</sup>, Adeline Simon<sup>1,5</sup>, Ricardo C. Rodríguez de la Vega<sup>1,5</sup>, Monika Coton<sup>1,3</sup>, Emmanuel Coton<sup>1,3</sup>, Marie Foulongne-Oriol<sup>1,6‡</sup>, Antoine Branca<sup>1,7‡</sup>, Tatiana Giraud<sup>1,1*</sup>

1. Ecologie Systématique Evolution, IDEEV, Gif-sur-Yvette, France, 2. Laboratoire Interprofessionnel de Production—SAS L.I.P., Aurillac, France, 3. Univ Brest, Laboratoire Universitaire de Biodiversité et Ecologie Microbienne, Plouzané, France, 4. Université Clermont Auvergne, INRAE, Vetagro Sup, Aurillac, France, 5. Université Paris-Saclay, INRAE, Palaiseau, France, 6. INRAE, MycSA, Mycologie et Sécurité des Aliments, Villenave d'Ornon, France, 7. Université Paris-Saclay, CNRS, IRD, UMR Évolution, Génomes, Comportement et Écologie, Gif-sur-Yvette, France

* These authors contributed equally to this work.

‡ These authors jointly supervised the study.

* thibaultcaron864@gmail.com (TC); tatiana.giraud@universite-paris-saclay.fr (TG)

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

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

OPEN ACCESS

**Citation:** Caron T, Crequer E, Le Piver M, Le Prieur S, Brunel S, Snirc A, et al. (2025) Identification of quantitative trait loci (QTLs) for key cheese making phenotypes in the blue-cheese mold *Penicillium roqueforti*. PLoS Genet 21(4): e1011669. https://doi.org/10.1371/journal.pgen.1011669

**Editor:** Paula Gonçalves, FCT NOVA: Universidade Nova de Lisboa Faculdade de Ciencias e Tecnologia, PORTUGAL

**Received:** March 25, 2024

**Accepted:** March 28, 2025

**Published:** April 15, 2025

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pgen.1011669

**Copyright:** © 2025 Caron et al. This is an open access article distributed under the terms of the Creative Commons Attribution License,

## Abstract

Elucidating the genomic architecture of quantitative traits is essential for our understanding of adaptation and for breeding in domesticated organisms. *Penicillium roqueforti* is the mold used worldwide for the blue cheese maturation, contributing to flavors through proteolytic and lipolytic activities. The two domesticated cheese populations display very little genetic diversity, but are differentiated and carry opposite mating types. We produced haploid F1 progenies from five crosses, using parents belonging to cheese and non-cheese populations. Analyses of high-quality genome assemblies of the parental strains revealed five large translocations, two having occurred via a circular intermediate, one with footprints of *Starship* giant mobile elements. Offspring genotyping with genotype-by-sequencing (GBS) revealed several genomic regions with segregation distortion, possibly linked to degeneration in cheese lineages. We found transgressions for several traits relevant for cheese making, with offspring having more extreme trait values than parental strains. We identified quantitative trait loci (QTLs) for colony color, lipolysis, proteolysis, extrolite production, including mycotoxins, but not for growth rates. Some genomic regions appeared rich in QTLs for both lipid and protein metabolism, and other regions for the production of multiple extrolites, indicating that QTLs have pleiotropic effects. Some QTLs corresponded to known biosynthetic gene clusters, e.g., for the production of melanin or extrolites. F1 hybrids constitute valuable strains for cheese producers, with new traits and new allelic combinations, and allowed identifying target genomic
which permits unrestricted use, distribution,
and reproduction in any medium, provided the
original author and source are credited.

**Data availability statement:** Reads and
assemblies were deposited at GenBank under
BioProject PRJNA1079247 for LCP06136,
PRJNA1211835 for the eleven RxN offspring
and PRJNA1216917 for the GBS data. Other
available data are presented in S5 Table.

**Funding:** This study has been funded by the ANR-19-CE20-0002-02 Fungadapt ANR (to ECo), ERC Genomesfun 309403 Stg (To TG), ERC Blue Proof of Concept (To TG), Fondation Louis D (French Academy of Sciences) grants (To TG), the LIP SAS and the ANRT (Association Nationale Recherche Technologie) (To DR and MP). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** I have read the journal's policy and the authors of this manuscript have the following competing interests: TC, MLP, SB, DR and MP were employed by SAS LIP (laboratoire Interprofessionnel de Production), which produces starters for fermented food products, during the course of the study and therefore declare a competing financial interest. None of the other authors has any conflict of interest to declare, with the exception of TC employed by the LIP during his PhD, but not afterwards during his post-doc representing most of the present work.

regions for traits important in cheese making, paving the way for strain improvement.
The findings further contribute to our understanding of the genetic mechanisms
underlying rapid adaptation, revealing convergent adaptation targeting major gene
regulators.

Author summary

Understanding the genetic determinants underlying quantitative traits is crucial for our understanding of adaptation and for improving varieties in domesticated organisms. *Penicillium roqueforti* is the mold used worldwide for blue cheese making. The two domesticated cheese populations in *P. roqueforti* each display very little genetic diversity, but are highly differentiated. We therefore generated progenies between these populations despite reduced fertility. As expected, offspring displayed high variation in multiple traits important for cheese making, involved in aspect, flavor or food security, such as color, lipolysis, proteolysis and extrolite production. We detected genetic determinants for most of these traits, paving the way for generating new varieties. Several genomic regions impacted multiple traits, as often found in domesticated organisms. These findings further contribute to our understanding of the genetic mechanisms underlying rapid adaptation to new environments, likely involving master gene regulation impacting multiple genes.

**Introduction**

Many traits of agricultural or evolutionary significance are quantitative and genet-
ically complex, involving multiple genes that interact with one another and the
environment [1]. Genomic loci showing diversity associated with quantitative trait
variation are called QTLs, for “quantitative trait loci”. The QTL approach has largely
advanced our understanding of the genomic architecture of important traits, for
example in crops and cattle [2–4]. Studies of the genomic architecture of adaptation
in crops have revealed that traits involved in domestication were often controlled by
only a few QTLs with pleiotropic effects, often being major regulators [5–16].

In the domesticated fungi used for food production, QTL mapping proved to be
an efficient tool for identifying the genes affecting phenotypes impacting technolog-
ical performances of industrial yeast strains (Saccharomyces cerevisiae), such as
thermotolerance, chemical resistance, and phenotypes associated with the fermen-
tation process, e.g., volatile compound and ethanol production [17–23]. In addition,
traits have also been reported to be impacted in yeasts by aneuploidy [24] and by
genomic rearrangements [25]. QTL mapping has also been successfully used in the
button mushroom *Agaricus bisporus*, for identifying the genomic regions controlling
color and spore number [26,27].
<figure><img src="image_5.png" /></figure>

*Penicillium roqueforti* is the mold used for the maturation of all types of blue cheeses worldwide, responsible for the typical blue-veined aspect, and contributing to their specific flavor and aroma, in particular through high levels of proteolytic and lipolytic activities [28–30]. When this study was performed, four populations had been identified in *P. roqueforti*, two populations being used to produce cheeses (named Roquefort and non-Roquefort) and two populations being found in molded silage, lumber or spoiled food [31]. The two populations used for cheese maturation each show footprints of bottlenecks and a domestication syndrome, with phenotypes contrasting with those of the non-cheese populations and beneficial for several important aspects of cheese safety, appearance and flavor [31–33].

The Roquefort population, found primarily in Roquefort protected-designation-of-origin cheeses, exhibits traits beneficial for pre-industrial cheese production, e.g., relatively slow growth in cheese medium and high spore production on bread medium, and produce high quantity and diversity of positive aromatic compounds in cheeses [31,32,34]. The other cheese population, named non-Roquefort, corresponds to a single clonal lineage, and is found in most types of blue cheese worldwide but rarely in Roquefort PDO cheeses. The non-Roquefort population displays phenotypes more suited for industrial cheese production, such as a faster growth in cheese, higher tolerance to salt and to lactic acid compared to other populations [31,33,35]. The non-Roquefort clonal lineage acquired large genomic regions by horizontal transfers, in particular two very large regions called *Wallaby* and *CheesyTer* [31,36], mediated by giant *Starship* mobile elements [37]. *Starships* are involved in adaptation in fungi domesticated for food fermentation [36,38], and also in plant pathogenic fungi [39,40]. The nucleotide diversity has been found very low in the non-Roquefort cheese population (π = 0.000115), about three-fold higher in the Roquefort population (π = 0.000403), and an order of magnitude higher in the lumber/spoiled food and silage populations (π values of 0.00106 and 0.00160, respectively; [31,33]).

*Penicillium roqueforti* can secrete metabolites, also known as extrolites, including some considered as toxins (e.g., Roquefortine C and PR toxin). However, these compounds are either absent in blue cheeses, or present at such low concentration that they do not pose acute health hazards [41–43]. The non-Roquefort lineage has lost the ability to secrete mycophenolic acid (MPA), an immunosuppressant used for preventing transplant organ rejection [44], due to a deletion in a key gene (*mpaC*) of its biosynthesis pathway [33,45].

Identifying the genomic regions controlling the phenotypes contrasting between the four *P. roqueforti* populations thriving in different niches has important fundamental implications for our understanding of the genomic architecture of adaptation and also applied consequences for strain improvement by marker-assisted selection in progenies. Due to the clonal structure and low diversity of the domesticated cheese lineages [31], identifying the genetic determinants for cheese-making traits through phenotype/genotype association requires producing a recombinant population, which would also generate new allelic combinations and phenotypic diversity. Fortunately, sexual reproduction can be induced in *P. roqueforti* [46]. The two cheese clonal lineages carry opposite mating types and some cheese strains can be crossed, despite the general degeneration for sexual reproduction in the cheese populations [47]. Furthermore, the haploid life cycle of *P. roqueforti* enables carrying out association genetics in first-generation offspring.

Here, we aimed at generating new allelic combinations in cheese *P. roqueforti* populations and at identifying the genomic regions controlling the phenotypic differences between *P. roqueforti* populations. For this goal, we produced five sexual F1 haploid progenies by crossing six fertile haploid *P. roqueforti* strains from the four originally identified populations. We studied in the progenies several phenotypes important for cheese making, i.e., lipolysis, proteolysis, growth, color and the production of the main known extrolites: roquefortine C, PR toxin, with its intermediates eremofortins A and B, as well as andrastin A, (iso)-fumigaclavine A and MPA [48,49]. We analyzed 1,073 offspring based on single nucleotide polymorphisms (SNPs) obtained using genotype-by-sequencing (GBS). Because genotyping revealed the presence of the two parental alleles in some genomic regions in the otherwise haploid offspring, we also generated high-quality genome assemblies for the parental strains and eleven offspring, to compare the parental genomes and identify genomic translocations. We built genetic maps and ran QTL analyses to identify the genetic architecture of the phenotypes important for cheese making.
## Results

### Five fertile crosses involving at least one cheese strain

We chose the five most fertile crosses, based on qualitative assessment of the number of cleistothecia, among 17 trials involving in each case at least one cheese strain, chosen based on fertility data from a previous study [47]. Pairwise nucleotide identity estimates between parental strains ranged from 99.48% to 99.63% across crosses (Table 1). We identified the isolated spores that were actually recombinant ascospores and not asexual conidia by using 11 microsatellite markers in the Roquefort x non-Roquefort cross (386 spores were excluded out of 1620 tested spores). We focused on the cross between the two cheese populations, i.e., Roquefort (R; LCP06136 strain) and non-Roquefort (N; LCP06133 strain), as it may be the most interesting for strain improvement in the cheese making context, and we isolated 387 recombinant offspring (Table 1). The other crosses each involved a cheese strain (either R or N) and a strain from a non-cheese population, either silage (S) or lumber/food spoiler populations (L). We isolated between 157 and 185 recombinant offspring for each of these four crosses (Table 1).

### Translocations between parental genomes

The comparison of high-quality genome assemblies revealed five large genomic chromosomal rearrangements between the parents in four crosses (LxN, RxN, RxB, RxL; Fig 1), which had impacts on genetic maps. These translocations ranged from 60 to 530 kb in size. The two silage parents and the non-Roquefort parent showed syntenic genomes between each other, suggesting that they kept the ancestral arrangement order. We could therefore infer that four of the

<table><caption>Table 1. The five crosses performed using six strains of Penicillium roqueforti from different populations (R for Roquefort, N for non-Roquefort, L for lumber/food spoiler, S for silage/food spoiler) with opposite mating types (MAT1-1 and MAT1-2), with the number of isolated and phenotyped offspring, the number of marker obtained from GBS (genotyping-by-sequencing) and the average nucleotide identity (ANI) between parental strains. LCP: « Laboratoire de Cryptogamie, Paris ».</caption><thead><tr><th>Cross</th><th>MAT1-1</th><th>MAT1-2</th><th>Number of ascospores isolated</th><th>Number of offspring GBS-genotyped</th><th>Number of offspring phenotyped</th><th>Number of off-spring color-phenotyped</th><th>Number of off-spring lip-olysis phenotyped</th><th>Number of off-spring prote-olysis phenotyped</th><th>Number of off-spring (myco-toxins) used for map construction</th><th>Number of marker</th><th>ANI (%)</th></tr><tr><th></th><th>Population</th><th>Strain</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>RxN</td><td>Roquefort</td><td>LCP06136</td><td>non-Roquefort</td><td>LCP06133</td><td>389</td><td>384<sup>a</sup></td><td>307</td><td>311</td><td>387</td><td>382</td><td>274</td><td>358</td><td>2462</td><td>99.578</td></tr><tr><td>RxB</td><td>Roquefort</td><td>LCP06136</td><td>silage/food spoiler</td><td>LCP06043</td><td>157</td><td>157</td><td>157</td><td>157</td><td>157</td><td>156</td><td>0</td><td>148</td><td>2902</td><td>99.593</td></tr><tr><td>RxL</td><td>Roquefort</td><td>LCP06136</td><td>lumber/food spoiler</td><td>LCP06037</td><td>185</td><td>185</td><td>129</td><td>131</td><td>150</td><td>141</td><td>120</td><td>159</td><td>2943</td><td>99.630</td></tr><tr><td>SxN</td><td>silage/food spoiler</td><td>LCP06059</td><td>non-Roquefort</td><td>LCP06133</td><td>176</td><td>176</td><td>92</td><td>92</td><td>150</td><td>150</td><td>75</td><td>160</td><td>2860</td><td>99.484</td></tr><tr><td>LxN</td><td>lumber/food spoiler</td><td>LCP06039<sup>a</sup></td><td>non-Roquefort</td><td>LCP06133</td><td>171</td><td>171</td><td>91</td><td>91</td><td>150</td><td>151</td><td>0</td><td>165</td><td>2618</td><td>99.511</td></tr></tbody></table>

aThe strain LCP06039 was found to be a mixture of two morphotypes ("dark" and "light"), which was detected after the cross had been performed.

bThese offspring have been genotyped with indels, in addition to GBS, as all other offspring.

https://doi.org/10.1371/journal.pgen.1011669.t001
<figure><img src="image_6.png" /></figure>

translocations (L1, L2, L3 and L4) occurred in the lumber population while the fifth one, R1, occurred in the Roquefort population (Fig 1).

The L1 and L4 regions displayed internal rearrangements between their copies in different genomic locations: the last 3' portion of each of the L1 and L4 regions in the parents with the ancestral arrangement was located in 5' of the L1 and L4 regions in the translocated parent (S1 Fig). This internal rearrangement suggests that the translocation occurred via a circular intermediate mechanism, with an insertion cutting site in each of the L1 and L4 regions different from their excision site (S1 Fig), as suggested to be the case for *Starship* mobile element movements [50]. As a matter of fact, we detected in the L4 translocation tyrosine recombinases (i.e., *Starship* "captains") and a DUF3723 domain, being all typical of *Starship* elements, using the dedicated *Starfish* pipeline [51].

Both parental alleles for the R1 region insertion were found together in 38, 31 and 24% of the RxN, RxL and Rxs off-spring in GBS data, respectively. Only the R1 region was found in duplicate copies in these progenies, not the whole chromosome(s), indicating a lack of whole chromosome aneuploidy. The Oxford Nanopore long-read mapping depths confirmed the presence of two R1 region copies in the six analyzed offspring from the RxN cross with both R1 region parental

<figure><img src="image_7.png" /><figcaption>**Fig 1. Genome comparisons between parental strains for each cross.** Chromosomes are represented with colors corresponding to the population-of-origin of the strains: light orange, blue, purple and green for silage (S), non-Roquefort, (N) Roquefort (R) and lumber (L) populations, respectively. From outside to inside, the tracks represent the i) GC content, ii) translocated regions of more than 50kb between parental genomes in red (with their IDs, the first letter indicating the population-of-origin of the parent in which the translocation occurred), specific regions that are lacking in other genomes (i.e., acquired horizontally transferred regions) with their color indicating the populations in which they are present, putative centromeres in black, and large recombination cold spots in light yellow, iii) predicted genes in gray and iv) transposable elements in green. Orthologs between genomes are linked one to each other, with blue, orange and red links, indicating synteny, inversions or >50kb translocations, respectively.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1011669.g001
<figure><img src="image_8.png" /></figure>

alleles detected in GBS data and the lack of whole chromosome aneuploidy. Similarly, the long-read mapping confirmed that the five offspring with a single parental allele in GBS data in the R1 region carried a single copy of the translocated region in their genomes. The absence of this translocated region is likely lethal as we did not detect any offspring with no R1 copy, while some recombination events should have generated offspring with none of the two copies. We detected no offspring with both parental alleles nor an absence of the region for the translocations L1 to L4, indicating that both the duplication or absence of these regions lead to inviability.

We checked that the presence of the two R1 copies in haploid genomes was stable after culture and replication of 24 offspring of the RxN cross carrying the two parental alleles. After cultivating these 24 offspring from the RxN cross for one week on plates and transferring conidia from the edge of the colony to new plates 19 times, PCRs showed that the cultivated lineages still carried the two R1 parental alleles. We did not detect any RIP footprints (C-to-T repeat-induced point mutations) in the R1 regions of the six offspring with a duplicated R1 region and with high-quality genome assemblies. However, RIP footprints would only be expected in the following generation (F2), after a sex event involving genomes with the duplicated region, as RIP is known to occur during the short dikaryotic phase of sexual reproduction in ascomycetes [52].

## Genetic maps, recombination pattern and segregation biases

We identified between 2,462 and 2,943 reliable markers in the progenies across the five crosses, evenly distributed along the non-Roquefort reference genome (LCP06133 strain). For constructing genetic maps, we excluded redundant markers, i.e., those physically close and with identical segregation patterns in a given cross. We thus constructed the genetic maps with 1,448 markers for the RxN cross, and 676–785 markers for the other four crosses with fewer analyzed offspring (Table 2). The genetic maps for the five crosses displayed lengths between 775 and 913 cM, with a mean of 857 cM (Table 2). We detected four linkage groups in each genetic map, corresponding to the four chromosomes in the parental genome assemblies (Fig 2). The mean number of crossing-overs per chromosome was 2.15, with a maximum of 2.58 for chromosome 1, the longest chromosome, and a minimum of 1.25 for chromosome 4, the smallest chromosome. The mean recombination rate was estimated between 27.0 cM/Mb and 31.8 cM/Mb, with a mean of 29.8 ± 1.8 cM/Mb (Table 2). The LxN cross had a particularly low mean recombination rate, mainly due to the presence of two cold spots of recombination on chromosome 1 (indicated by a plateau in the plot of genetic map positions against genomic positions; empty rectangles in Fig 2), with local recombination rates of 6.7 cM/Mb and 3.1 cM/Mb, respectively (Table 2). The edges of the two cold spots corresponded to the translocated regions in the chromosome 1 of the LCP06039 strain from the lumber/food spoiler population (regions L1, L2, L3, Figs 1 and 2). This is likely due to the inviability of offspring carrying either zero or two copies of the translocated region because of recombination events between their insertion sites (Fig 1). Other regions with low or zero recombination rates corresponded to horizontally transferred regions only present in the reference genome, putative centromeres, other translocated regions and TE-rich regions (Figs 1 and 2).

In the five crosses, large regions presented segregation biases, with under-representation of alleles from the cheese populations (i.e., either from the Roquefort or non-Roquefort parents; Fig 3). Such segregation biases could allow the purge of deleterious mutations in the cheese clonal lineages, but could also render the selection of valuable alleles that are in linkage with deleterious alleles in these regions more challenging. The segregation bias against the non-Roquefort parent allele in the second half of chromosome 4 (between 3.5 Mb and 4.2 Mb) was present in all crosses involving the non-Roquefort parent (Fig 3) and was very strong (>90% of the alternative allele in all three crosses). Such strong under-representation of alleles from the cheese population in all of these three crosses can be due to the lower fertility of domesticated populations [47], with likely deleterious alleles accumulated through clonal replication that lead to low viability of the ascospores carrying them. We observed another segregation bias with strong over-representation (> 90%) of silage parental alleles, in the Rxs and SxN crosses, in chromosome 1 from 4.2 to 4.5 Mb (Fig 3). The under-represented alleles, that of the cheese parents, displayed no segregation bias in the other crosses involving the cheese parents (Fig 3),
suggesting the presence of either a selfish element (e.g., a spore killer) or highly beneficial allele(s) for early growth in the corresponding region in the silage/spoiled food parents rather than a deleterious one in the cheese parent. This interpretation is supported by the over-representation of the silage alleles in two different crosses, and the lack of bias in the Rxl and LxN crosses in this region, while similar biases would be expected in all crosses if the cheese parents had deleterious alleles in this region. The segregation bias with over-representation of Roquefort alleles at the end of chromosome 1 is likely due to the lack of the R1 translocated region that leads to offspring inviability. Symmetrically, the alleles from the other parent tended to be over-represented at the R1 locus on chromosome 2, but not significantly so. This suggests that the location of R1 may be more advantageous on chromosome 1 than on chromosome 2.

<table><caption>Table 2. The five crosses performed using six strains of Penicillium roqueforti from different populations (R for Roquefort, N for non-Roquefort, L for lumber/food spoiler, S for silage/food spoiler), the linkage group and the total, the number of markers used for genetic map construction, the total length of the linkage group, the average and maximum spacing in cM between markers, the mean recombination rate in cM/Mb, and the reference chromosome and genome size in Mb.</caption><thead><tr><th>Cross</th><th>Linkage group</th><th>Number of markers</th><th>Mean number of crossing over</th><th>Length (cM)</th><th>Average spacing (cM)</th><th>Maximum spacing (cM)</th><th>Mean recombination rate (cM/Mb)</th><th>Reference size (Mb)</th></tr></thead><tbody><tr><td rowspan="5">SxN</td><td>L.1</td><td>213</td><td>2.74</td><td>274</td><td>1.3</td><td>8.8</td><td>32.21</td><td>8.51</td></tr><tr><td>L.2</td><td>242</td><td>2.49</td><td>249</td><td>1</td><td>12.8</td><td>29.97</td><td>8.31</td></tr><tr><td>L.3</td><td>189</td><td>1.97</td><td>197</td><td>1</td><td>8.8</td><td>25.51</td><td>7.72</td></tr><tr><td>L.4</td><td>94</td><td>1.19</td><td>120</td><td>1.3</td><td>8.8</td><td>28.47</td><td>4.22</td></tr><tr><td>overall</td><td>738</td><td>8.39</td><td>841</td><td>1.1</td><td>12.8</td><td>29.25</td><td>28.75</td></tr><tr><td rowspan="5">LxN</td><td>L.1</td><td>161</td><td>1.98</td><td>199</td><td>1.2</td><td>9.2</td><td>23.39</td><td>8.51</td></tr><tr><td>L.2</td><td>248</td><td>2.55</td><td>255</td><td>1</td><td>9.2</td><td>30.69</td><td>8.31</td></tr><tr><td>L.3</td><td>164</td><td>1.95</td><td>195</td><td>1.2</td><td>6.7</td><td>25.25</td><td>7.72</td></tr><tr><td>L.4</td><td>103</td><td>1.28</td><td>128</td><td>1.3</td><td>8.6</td><td>30.36</td><td>4.22</td></tr><tr><td>overall</td><td>676</td><td>7.75</td><td>776</td><td>1.2</td><td>9.2</td><td>26.99</td><td>28.75</td></tr><tr><td rowspan="5">RxN</td><td>L.1</td><td>466</td><td>2.78</td><td>279</td><td>0.6</td><td>7.3</td><td>32.79</td><td>8.51</td></tr><tr><td>L.2</td><td>435</td><td>2.48</td><td>248</td><td>0.6</td><td>8.7</td><td>29.85</td><td>8.31</td></tr><tr><td>L.3</td><td>378</td><td>2.22</td><td>223</td><td>0.6</td><td>9.6</td><td>28.88</td><td>7.72</td></tr><tr><td>L.4</td><td>179</td><td>1.25</td><td>126</td><td>0.7</td><td>8.7</td><td>29.89</td><td>4.22</td></tr><tr><td>overall</td><td>1458</td><td>8.74</td><td>875</td><td>0.6</td><td>9.6</td><td>30.43</td><td>28.75</td></tr><tr><td rowspan="5">RxA</td><td>L.1</td><td>191</td><td>2.74</td><td>275</td><td>1.4</td><td>11.7</td><td>32.32</td><td>8.51</td></tr><tr><td>L.2</td><td>242</td><td>2.47</td><td>248</td><td>1</td><td>13.9</td><td>29.85</td><td>8.31</td></tr><tr><td>L.3</td><td>227</td><td>2.55</td><td>256</td><td>1.1</td><td>12.4</td><td>33.15</td><td>7.72</td></tr><tr><td>L.4</td><td>115</td><td>1.33</td><td>133</td><td>1.2</td><td>7.5</td><td>31.55</td><td>4.22</td></tr><tr><td>overall</td><td>775</td><td>9.10</td><td>913</td><td>1.2</td><td>13.9</td><td>31.75</td><td>28.75</td></tr><tr><td rowspan="5">RxL</td><td>L.1</td><td>186</td><td>2.64</td><td>264</td><td>1.4</td><td>14.9</td><td>31.03</td><td>8.51</td></tr><tr><td>L.2</td><td>295</td><td>2.66</td><td>267</td><td>0.9</td><td>13.5</td><td>32.14</td><td>8.31</td></tr><tr><td>L.3</td><td>181</td><td>2.30</td><td>231</td><td>1.3</td><td>9.5</td><td>29.91</td><td>7.72</td></tr><tr><td>L.4</td><td>129</td><td>1.19</td><td>119</td><td>0.9</td><td>12.9</td><td>28.23</td><td>4.22</td></tr><tr><td>overall</td><td>791</td><td>8.79</td><td>881</td><td>1.1</td><td>14.9</td><td>30.64</td><td>28.75</td></tr><tr><td rowspan="5">Crosses mean</td><td>L.1</td><td>243.4</td><td>2.58</td><td>258.2</td><td>1.18</td><td>10.38</td><td>30.348</td><td>8.51</td></tr><tr><td>L.2</td><td>292.4</td><td>2.53</td><td>253.4</td><td>0.9</td><td>11.62</td><td>30.5</td><td>8.31</td></tr><tr><td>L.3</td><td>227.8</td><td>2.20</td><td>220.4</td><td>1.04</td><td>9.4</td><td>28.54</td><td>7.72</td></tr><tr><td>L.4</td><td>124</td><td>1.25</td><td>125.2</td><td>1.08</td><td>9.3</td><td>29.7</td><td>4.22</td></tr><tr><td>overall</td><td>887.6</td><td>8.55</td><td>857.2</td><td>1.04</td><td>12.08</td><td>29.812</td><td>28.75</td></tr></tbody></table>

https://doi.org/10.1371/journal.pgen.1011669.t002
<figure><img src="image_9.png" /></figure>

Transgression and heterosis in phenotype distributions

We measured several phenotypes important for blue cheese production, i.e., mycelium growth, colony color, lipolysis, proteolysis and extrolite production. We analyzed nine variables summarizing these traits in the five progenies, plus seven traits related to extrolite production in three crosses given time and cost limitations; the three crosses were Rxl, RxN, SxN, choosing crosses involving at least one cheese population. Overall, these measures resulted in 66 trait distributions (9x5+7x3). Among these 66 distributions, 28 did not significantly deviate from normality, 10 best fitted a unimodal distribution and 30 a bimodal distribution (S2 and S3 Figs and S1 Table). None of the phenotype distributions were trimodal or amodal (S1 Table).

<figure><img src="image_10.png" /><figcaption>**Fig 2.** Genetic and physical maps for the four chromosomes of *Penicillium roqueforti*, obtained by analyzing five inter-population crosses. Each panel corresponds to an individual chromosome, line graphs on the left represent linkage groups (one color per cross, with cross IDs as in Table 1; note that different parents were used for the silage and lumber populations between crosses); the lines link markers to chromosomes of the reference genome (LCP06133 non-Roquefort parent). In each panel, plots on the right represent the relationship between genetic and physical distance, with one plot per cross and the same color code as in left panels. At the bottom of each panel, the reference chromosome is represented, with filled black, red and blue rectangles representing the putative centromeres, the translocated regions and the horizontally transferred regions of more than 50kb when present, respectively. The red translocated regions are labeled as in Fig 2. The orange bar in chromosome 3 represents the mating-type locus position. The empty black rectangles indicate large recombination cold spots situated between translocated or parent-specific regions.</figcaption></figure>
The distribution of the five progenies on the principal component analysis based on trait values (S4 Fig) shows that they represent overall a higher phenotypic diversity than each cross considered separately. A few traits appeared positively or negatively correlated (Fig 4D), some associations being expected, such as between color traits (e.g., blue and red) or between the production of different targeted extolites [34,53–55]. Other associations were more surprising, such as the negative correlation between proteolysis rate on the one hand, and green and hue level on the other hand (Fig 4D).

We found heterosis in the progenies, computed as the difference in mean trait values between the parents and their progenies (S2 Table). Out of the 66 phenotype distributions, 53 displayed significant positive or negative heterosis. We found negative heterosis in all crosses for the production of PR toxin (the most toxic *P. roqueforti* mycotoxin) and andrastin A, indicating that progenies produced on average less of these two extolites than their parent mean. In contrast, we found positive heterosis in all crosses for MPA production, a less toxic extolite. All tested phenotypes

<figure><img src="image_11.png" /><figcaption>**Fig 3.** Representation of identified quantitative trait loci (QTLs) and segregation distortion along the four chromosomes of *Penicillium roqueforti* in five progenies. Only QTLs explaining more than 5% of the total variance are represented. The color of the bars indicates the considered trait while the letter indicates the parent having the alleles increasing the trait value. Dotted lines link the QTLs with significant interactions and color indicates the phenotype. The cross ID is indicated on the right, the first parent carrying the MAT1-1 mating type. (A) QTLs identified for phenotypes linked to lipolysis, proteolysis, growth and colony color for the five crosses (in lines, indicated on the right with the same code as in Table 1; note that different parents were used for the silage and lumber populations between crosses). Black curves represent the proportion in progenies of the allele from the parent carrying the MAT1-1 mating type. Transparent colored rectangles on curves indicate regions with significant segregation distortion, their color corresponding to the under-represented parental allele: blue, purple, green and yellow for non-Roquefort, Roquefort, lumber/spoiled food and silage/ spoiled food parents, respectively. The x-axis represents genomic physical positions. At the bottom, the four chromosomes of the reference genome (LCP06133) are represented, with rectangles representing genomic features: large horizontally transferred regions specific to the reference genome in blue, and translocated regions in the reference genome in red (with same IDs as in Fig 2). The location of the horizontally transferred regions (W for Wallaby and C for CheesytTer) are indicated. The location of the dihydroxynaphthalene (DHN) melanin production cluster is indicated by a black arrow. The large vertical empty black rectangles indicate pleiotropic QTL regions, with effects on multiple traits. (B) QTLs identified for extolite production (with colors corresponding to the different toxins), analyzed in three crosses (SxN, RxN and LxN). The locations of the gene clusters of the toxin biosynthesis pathways are indicated by arrows with the same color code (the PR toxin cluster being in gray, eremofortin A and B being intermediates).</figcaption></figure>
<figure><img src="image_12.png" /></figure>

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

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

Cross

RgL
RoS
RxN
LxN
SxN

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

D

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

Parental population

Lumber

non-Roquefort

Roquefort

Silage

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

Fig 4. Distributions of phenotypes in the five progenies in Penicillium roqueforti, for (A) the lipolysis rate (mm/day), (B) the proteolysis rate (mm/day) and (C) the relative green color channel of the progeny (RxL in red, RxS in beige, RxN in yellow, LxN in turquoise and SxN in pink). The black point and line represent the mean and standard deviation, respectively, for each progeny. The colored points represent the parental values (lumber in green, non-Roquefort in blue, Roquefort in purple, silage in orange). (D) Pearson correlation coefficients, as a color gradient, between trait values across offspring and crosses.

https://doi.org/10.1371/journal.pgen.1011669.g004

exhibited both positive and negative transgressions across all crosses, i.e., with some offspring displaying higher
values and others lower values than either parent. We found a large proportion of offspring with low or no PR toxin
production in the RxN and SxN crosses. We detected in all crosses negative transgression in roquefortine C production
level (i.e., the offspring mean was lower than the parental mean) and many offspring with no MPA production, which is
promising for the selection of hypotoxinogenic strains.
<figure><img src="image_18.png" /></figure>

## High values of phenotype heritability

Heritability can be estimated either as broad-sense heritability, accounting for total genetic variance (e.g., additive and epistatic), or narrow-sense heritability, only accounting for additive genetic variance. Narrow-sense heritability ($h^2$) was estimated for all phenotypic traits across the five crosses (RxD, Rxl, RxD, SxD and LxD) using a linear mixed model approach, based on the genetic relatedness between individuals and their phenotypes (S3 Table). Most traits analyzed exhibited relatively high narrow-sense heritability ($h^2$>50%) in at least one cross. Lipolysis and proteolysis rates, green proportion and extrolite production showed heritability estimates even above 70%, suggesting that these traits are primarily influenced by additive genetic effects, thus allowing selection in breeding programs. In contrast, the growth rate displayed consistently low heritability values, ranging from 9.6% (RxD) to 26.9% (SxD).

The broad sense heritability could be estimated for extrolite production. Estimates were high, ranging between 91.3% for roquefortine C and 99.6% for isofumigaclavine A, which is also promising for selection (S3 Table).

## Detection of QTLs with pleiotropic effects

We investigated QTLs for the nine traits associated with lipolysis, proteolysis, growth and color in all descendants, as well as for the seven traits associated with extrolite production specifically in the RxD, Rxl and SxD crosses (S4 Table). We found QTLs for all tested phenotypes except for growth rate (Fig 3 and S4 Table), indicating that marker-assisted selection for these traits may be performed for strain improvement and diversification. We detected 123 QTLs across all phenotypes and crosses, 109 of which explained 5% or more of the phenotypic variance in a given cross (hereafter called "major QTLs"). We identified 58 major QTLs for the eight parameters related to lipolysis, proteolysis and color, across the five crosses, with 8–17 QTLs per cross, a mean of 1.3 and maximum of 3 major QTLs per trait and cross. For extrolite production, we identified 13–23 major QTLs across the three progenies analyzed for these traits, with a mean of 2.4 and a maximum of 5 major QTLs per extrolite and cross (S4 Table and Fig 3).

The identified QTLs and the direction of their effects were in agreement with the past occurrence of selection for color, proteolysis, lipolysis and extrolite production in cheese populations. We indeed identified QTLs impacting colony color, bluer colonies being associated with cheese parent alleles, and in particular in the genomic region containing the dihydroxynaphthalene-melanin biosynthesis cluster, known to be involved in melanin production (Fig 3). The cheese parent alleles were also both associated with slower proteolysis and faster lipolysis, in connection with firmer texture and longer cheese storage and the production of typical blue cheese flavors, respectively. We identified QTLs for the production of three extrolites, MPA, PR toxin and roquefortine C, the non-Roquefort alleles being associated with lower production levels (Fig 3). This is consistent with a selection for bluer color, more efficient lipolysis and less efficient proteolysis in domesticated lineages, as well as lower extrolite production levels.

Most of the major QTL regions presented pleiotropic effects. The 109 major QTLs were indeed not homogeneously distributed across the genome, with regions instead carrying clustered QTLs for different phenotypes (Figs 3 and S5). We considered as pleiotropic the regions presenting overlaps in QTL intervals of at least two different trait classes in a given cross (proteolysis, lipolysis, color and extrolite production). We identified five such regions, among which three regions displayed pleiotropic impact on the four trait classes, across the different progenies, and also within some progenies for multiple traits (empty vertical rectangles in Figs 3 and S5). For the pleiotropic region on chromosome 2, the Roquefort allele was associated with faster lipolysis and lower proteolysis. This QTL region displayed significant interactions for lipolysis with two other QTLs, in chromosomes 3 and 4; such interactions between loci in determining a trait value constitute positive epistasis. The same region also impacted color. For the pleiotropic regions detected in chromosomes 1 and 3, the non-Roquefort alleles were associated with slower proteolysis, faster lipolysis and bluer colony color (S4 Table and Fig 3). We further found epistasis for other types of phenotypic traits, particularly between different QTLs controlling proteolysis (S4 Table and Fig 3). QTL interactions had the same sign (positive or negative) across the
<figure><img src="image_19.png" /></figure>

various crosses in some cases and different signs in other cases, indicating that epistatic interactions are dependent on the genetic background. Such epistasis could impair selection, as recombination breaks up allelic combinations, while specific allelic combinations are beneficial under positive epistasis.

We also identified QTLs for extrolite level production, including some outside of their known biosynthesis gene clusters, suggesting that these QTLs correspond to *trans*-acting regulators. For three extrolites (MPA, PR toxin and roquefortine C), we identified QTLs with intervals included in their biosynthesis cluster, with non-Roquefort alleles associated with lower production levels, consistent with selection for an hypotoxinogenicity in this population. The lower MPA production level associated with the non-Roquefort allele is likely due to the deletion in the *mpaC* gene that was identified in the non-Roquefort population [34,45]. We also found a QTL in a region including the PR toxin biosynthesis gene cluster, with lower production level of the PR toxin but accumulation of its production intermediates, eremofortins A and B, associated with the non-Roquefort allele at the PR toxin biosynthesis cluster. This suggests a disruption in the cluster, driving low production levels of the toxin and accumulation of its intermediates in offspring presenting the non-Roquefort allele, likely due to a premature stop codon in the ORF11 gene of the cluster [45]. Most of the QTLs associated with extrolite production were not located within their biosynthesis gene cluster, which suggests the implication of *trans*-acting regulators. These QTLs were often shared between different extrolites, suggesting a co-regulation of their pathways, which could facilitate the selection for low extrolite production in strains bred for cheese production. In particular, isofumigaclavine A production shared five of its eight major QTLs with andrastin A production, suggesting co-regulation of the production of the two extrolites. We detected QTLs for all tested extrolites, except roquefortine C, in the R1 region for its two insertion sites (Fig 3), with lower production associated with the presence of the R1 region. Such a QTL clustering suggests the presence of a master regulator of extrolite production in the R1 region. We, in fact, found an ortholog of a master regulator gene in the R1 region, namely the *srk1* gene, regulating osmotic and oxidative responses in fungi [56].

We detected an andrastin A QTL next to the *Wallaby* horizontally transferred region (HTR) in the RxN cross and a proteolysis QTL in the region containing the *CheesyTer* HTR in each of the SxN and LxN crosses (Fig 3). Multiple other phenotypes have previously been suggested to be controlled by these HTRs, such as milk sugar metabolism and interactions with other microorganisms [36], but the indel nature of these HTRs may render hard to detect their associations with phenotypes, as markers are lacking in one of the parents.

We also checked whether QTL regions encompassed other regulators known to be involved in extrolite production (*prlaeA*, [55], *pcz1*, [57], and *sfk1*, [53]), as well as other essential phenotypes in *P. roqueforti*, such as conidiation or growth rate. We located *pcz1* in the pleiotropic region of chromosome 2, encompassing QTLs for andrastin A production and color, melanin being involved in both color and conidiation [57]; *pcz1* could therefore be a candidate for the QTLs controlling andrastin A production and color via expression regulation. While no QTL interval overlapped with *prlaeA*, a QTL region controlling MPA production in one cross encompassed *sfk1*. However, this QTL region overlaps with a shorter MPA QTL present in another cross without *sfk1* in the confidence interval, which suggests no implication of *sfk1* in MPA production regulation.

The QTL intervals were too large to identify particular genes or functions beyond the horizontal gene transfers, biosynthesis gene clusters of the targeted extrolites or other a priori candidate genes. Indeed, with on average one gene every 3kb in the annotated Roquefort parental strain and a median QTL length of 760 kb, half of the QTL regions are expected to encompass more than 250 genes. We therefore only tried to identify candidate genes in the smallest QTL intervals explaining more than 20% of the phenotypic variance. We found an enrichment in genes involved in carbohydrate catabolism for a QTL associated with color in the pleiotropic region of chromosome 2 (g3621 to g3634), and a pepsin gene (endopeptidase, g5778) for a QTL associated with proteolysis in chromosome 3. Carbohydrate catabolism and pepsin are of paramount importance for secondary starters, on one hand, to use the remaining lactose residues to establish themselves early in the cheese matrix [28], and, on the other, to degrade caseins and therefore to participate in proteolysis [58].
<figure><img src="image_20.png" /></figure>

Discussion

Variability in phenotypes in progenies, with high heritability, heterosis and positive transgression for most traits

We found a large phenotypic variation, transgression and heterosis in progenies for multiple relevant phenotypes for blue cheese production, i.e., mycelium growth, lipolysis, proteolysis and extrolite production. The generation of new allelic combinations is particularly important in *P. roqueforti*, as the two main cheese populations, Roquefort and non-Roquefort, have lost most of their genetic diversities due to recent bottlenecks driven by strong selection in industrialization times [31]. These clonal lineages have degenerated in terms of sexual fertility and probably for other traits of interest for cheese making [47], and the use of progenies could purge accumulated deleterious mutations [59].

We found high values of heritability for most traits, in particular for toxin production, as previously reported for aflatoxin production in the fungus *Aspergillus flavus* [60], an opportunistic pathogen causing aspergillosis disease in humans. High values of heritability were also found for many traits in the yeast *S. cerevisiae*, in particular for traits important for beverage production [61,62]. High heritability was also reported for traits related to fitness in wild fungal species, such as the plant pathogens *Monilinia vaccinii-corymbosi* and *Rhizoctonia solani* [63,64].

We observed transgressive segregation in *P. roqueforti* progenies. Some offspring indeed displayed more extreme trait values than both parents in key traits such as color, metabolite production, as well as proteolytic and lipolytic activities. Transgression, where offspring exhibit extreme phenotypes, with values beyond those of their parents, has been documented in a wide range of organisms. For example, transgression has been leveraged in the yeast *S. cerevisiae* to enhance fermentation ability and improve morphological characteristics [65–67]. Transgressive effects have also been reported in crops like rice, sorghum, and lettuce for agronomic features [68–71] and in the button mushroom *A. bisporus* for yield-related traits and color [26].

Translocations and segregation biases

The high-quality genome assemblies led to the identification of five translocations between the parental strains, having occurred independently in three lineages. These translocations and the resulting copy-number variation in some offspring could be one of the reasons why the fertility of crosses is low and variable among crosses in *P. roqueforti* [47]. In *S. cerevisiae*, non-reciprocal translocation can result in segmental duplication in offsprings, leading to genomic imbalances and reduced viability [72]. A similar reduction in offspring viability due to translocations has also been reported in *Neurospora crassa* [73]. We detected lower production levels for the targeted metabolites in offspring carrying two copies of the R1 region. These offspring may therefore be highly valuable for cheese making, especially as the presence and locations of the translocations were stable through multiple replications and culture events.

The L4 translocation appeared to have occurred via a circular form and encompassed domains typical of mobile *Starship* elements, indicating that additional *Starships* are present in *P. roqueforti* compared to the previously detected ones, shown to be involved in adaptation [36,37]. The L1 translocation, having also occurred by a circular intermediate, may also be a *Starship*, even if the *Starfish* pipeline did not detect the typical domains, but it can sometimes miss them. The other translocations may also be *Starships* or have occurred without being mobile elements.

Chromosomal translocations between domesticated lineages and their wild relatives, possibly affecting crosses and thereby breeding, have been reported in other fungi, for example in the yeast *Lachancea cidri* [74], but also in plants: soybean [75], wheat [76], oat [77], quinoa [78] and chili pepper [79]. These translocations sometimes provide beneficial phenotypes, such as fermentation and ethanol tolerance in yeast [74]. Structural variation has however been suggested to be mostly deleterious in domesticated lineages and to represent a cost of domestication, for example in rice [80] and grapevine [81].
<figure><img src="image_21.png" /></figure>

Crossing-over numbers per chromosome and recombination rate

In eukaryotes, the mean number of crossing-overs is often around two per chromosome, with variability across crosses [82], which is in agreement with our findings in *P. roqueforti*. Among fungi, the number of crossing-overs per chromosome is between 4 and 10 in ascomycetous yeasts, such as *Saccharomyces cerevisiae* [83] and *Komagataella phaffii* [84], above 15 in the ascomycete *Aspergillus fumigatus* [85], between 0.25 and 1.8 in the basidiomycete *Lentinula edodes* [86], and between 0.19 and 1.33 in another basidiomycete, *Hericium erinaceus* [87]. Thus, our estimates in *P. roqueforti* of 2.15 crossing-over per chromosome, on average, is among the lowest numbers of crossing-overs per chromosome described in ascomycete fungi. Our estimated recombination rate in *P. roqueforti* (29.81 cM/Mb) was also lower than the average reported in fungi of 48.68 cM/Mb [88], but closer to the average for eukaryotes (22.93 cM/Mb; [88]). Fungi indeed generally display high recombination rates compared to other eukaryotes, mainly because of their small genome size, given the negative correlation between genome size and recombination rate [88].

QTL identification, strain breeding and maintenance of genetic diversity in domesticated populations

We identified QTLs for all phenotypes but growth. There were only a few major QTLs per trait and cross, which should facilitate strain breeding. The QTLs for color, lipolysis and proteolysis on the one hand, and for extrolite production on the other hand, clustered in the same genomic regions, suggesting pleiotropy, or clustering of genes with effects on these different pathways. The control of traits involved in domestication by only a few QTLs and their pleiotropy have been reported in animals and plants, often mainly due to major transcription regulators, human selection having unconsciously targeted "masterminds" [10], but also to the genetic determinant of important crop traits being clustered in genomes [9,15]. Here, we show that such pleiotropic regions are present in two independently domesticated lineages of *P. roqueforti*, with the same type of effect on the phenotype, suggesting evolutionary convergence [89] and selection of major regulators with pleiotropic effects. The pleiotropy of many QTLs can either facilitate or impair the selection for desired trait combinations. For example, antagonistic pleiotropy impaired industrial strain optimization in *S. cerevisiae* and *Schizosaccharomyces pombe*, because of a trade-off between stress resistance and ethanol production [90].

The discovery of additional, non-pleiotropic QTLs in *P. roqueforti* indicates the potential for a finer selection on individual traits without negatively impacting others. QTLs can serve as valuable markers for breeding programs, as previously done for improving fermentation traits in yeast [91]. Marker-assisted selection has been successfully applied to mitigate phenolic off-flavor or hydrogen sulfide production in *S. cerevisiae* [92] and for the rapid selection of the sporeless trait in the basidiomycete *Pleurotus pulmonarius* [93].

For strain breeding in *P. roqueforti*, further efforts could include integrating the recently discovered Termignon cheese population. The Termignon population is indeed genetically and phenotypically distinct from other cheese populations, likely representing remnants of an ancestral-like population of early domesticates [33]. This population thus offers a unique genetic reservoir for strain breeding and for enhancing diversity in cheese production. In general, more efforts are needed in agro-ecosystems for sustainable use of crops and domesticated micro-organisms, to maintain extant genetic diversity, and reintroduce diversity from wild populations or even closely related species [94–96].

Materials and methods

Progeny production

The production of progenies was obtained following the protocol in [46] with some modifications. Six *P. roqueforti* strains (Table 1) were selected for their fertility as observed in a previous study [46] and based on their phenotypic differences estimated from proteolysis and lipolysis measurements [31]. They were crossed in all possible pairs given their mating types, either MAT1-1 (three strains) or MAT1-2 (three strains) (Table 1).
<figure><img src="image_22.png" /></figure>

For each cross, three Petri dishes were inoculated on an oat medium supplemented with biotin after sterilization (6.4 µg/L; [97]). On each Petri dish, 5 µL of an uncalibrated spore suspension was used to inoculate agar sections at opposite sites on the Petri dish, located perpendicular to the sections inoculated with aliquots of conidia of the opposite sexual type (Fig 5A; [98]). Petri dishes were sealed with food-grade plastic film and incubated upside down at 19°C in the dark. The contact zone between strains (Fig 5A) was examined with a binocular loupe and an optical microscope regularly for four to six weeks to check for the presence of cleistothecia (sexual structures). When cleistothecia were found, three to five of them were sampled, crushed with tweezers, and observed under an optical microscope with lactophenol cotton blue dye solution for ascospore detection (Fig 5B and 5C). We took the pictures with a Nikon DS-Fi2 camera, on an inverted optical microscope with apodized phase contrast (Nikon France, Champigny-sur-Marne, France). When ascospores were observed, all the cleistothecia of the corresponding Petri dish were collected in a 0.05% Tween-20 water suspension (P1379-100mL Sigma Aldrich, Saint Quentin Fallavier, France). Because cleistothecia are covered with conidia derived from asexual reproduction, we developed a treatment to kill them before releasing ascospores from cleistothecia. Approximately 500 µL of the cleistothecia suspension was pooled with 500 µL of a commercial 2.6% sodium hypochlorite solution and gently agitated for 20 seconds; then 1 mL of water was added to stop the sodium hypochlorite action. Cleistothecia were washed twice with one-minute centrifugation at 7,000 rpm, removal of the supernatant and resuspension in Tween 20. Then, we crushed cleistothecia with a sterile piston (Piston Pellet Eppendorf, Montesson, France) and filtered residues with a 40 µm filter (cell screen, Greiner, Les Ulis, France). Single-ascospore isolation was carried out with a dilution method on malt extract agar Petri dishes (MEA, [99].

<figure><img src="image_23.png" /><figcaption>**Fig 5.** Observations of *Penicillium roqueforti* crosses and examples of phenotypic trait determination. (A) Petri dish with a cross between the LCP06136 (MAT1-1; Roquefort; at the top and bottom) and LCP06173 (MAT1-2; non-Roquefort; at left and right) strains. Sexual structures are formed at the contact zones between the two strains. (B) Cleistothecia (sexual structures of *Penicillium roqueforti*) shown by black arrow. (C) Asci containing ascospores shown by the black arrow. (D) Petri dish used for estimating colony growth rate, by counting the number of pixels occupied by the fungal colony (one offspring grown on malt medium for 120 hours). (E) Two lipolysis test tubes showing the lysis dynamics, for two different strains, with marks on the tubes indicating the limit of lipid medium degradation at different times (blue, red and blue marks, drawn at 0, 7 and 14 days, respectively). (F) Petri dish on which an offspring grew on raw sheep milk medium during 13 days for color measure.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1011669.g005
<figure><img src="image_24.png" /></figure>

## Recombinant offspring detection

We chose five fertile crosses for further analyses (Table 1). Genomic DNA was extracted from fresh mycelium and conidia suspension after single-ascospore isolation and growth for five days on MEA media using the Nucleospin soil kit (Macherey-Nagel, Düren, Germany). For checking that isolated spores were recombinant offspring and not asexual conidia, we used 11 microsatellite markers [46] labeled with fluorescence, out of which between five and nine were polymorphic depending on the cross. These markers were amplified with the Multiplex PCR kit (Qiagen, Les Ulis, France) using a touchdown program with an initial denaturation of 15 min at 95°C, 35 cycles of 30 s at 94°C, a decrease of 1°C every 90 s from 60 to 50°C, and 60 s at 72°C. The PCR program ended with a final 30 min extension step at 60°C. Genotyping by capillary fractionation electrophoresis was performed at INRAe Clermont-Ferrand (INRAe Platform GENTYANE, Clermont-Ferrand, FRANCE). The profiles obtained were analyzed with the GENEMAPPER v4.0 software (Applied Biosystem, Villebon-sur-Yvette, France) to detect recombination between parental genotypes. Only strains with recombinant genotypes were retained, i.e., carrying alleles from one parent at some markers and alleles for the other parent at other markers.

## Genome sequencing, assembly and analysis

We generated long-read-based genome assemblies by sequencing the genomes of the LCP06136 parental strain, as well as of 11 F1 offspring, with Oxford Nanopore MinION technology using an R9 flow cell, a high-quality genome assembly of the other parental strain (LCP06133) being available [33]. We also generated genome assemblies for the other parents, based on Illumina sequencing and assemblies guided by the Oxford Nanopore assemblies of the LCP06136 and LCP06133 parents (S5 Table).

DNA was extracted from mycelium and conidia with the Nucleospin soil kit for the progeny and with the NucleoBond High Molecular Weight DNA kit (Macherey-Nagel, Düren, Germany) for the parental strain LCP06136, with mechanical disruption of about 30 mg of lyophilized mycelium with two tungsten beads (3 mm diameter) for 5 min at 30 Hz. The Nanopore library was prepared with the SQK-LSK109 ligation sequencing kit, starting with 1.5 µg DNA, and sequencing was performed in-house with a MinION combined with a MinIT (version 19.06.9), a 72 hour run, and with the FastCALLing baseCALLing algorithm. The genome of the LCP06136 strain was sequenced alone in a R9 flow cell, whereas the 11 F1 offspring genomes of the LCP06133 x LCP06136 cross were sequenced in a 12-plex with the Rapid Barcoding Sequencing kit SQK-RBK004. We assessed run quality using Porechop v0.2.3_sequan2.1.1 [100], and when necessary, demultiplexed the Nanopore raw reads with the default parameters.

De novo assemblies of the genomes were constructed from both Illumina and Nanopore reads. For the LCP06136 parental genome, the raw Nanopore reads were trimmed and assembled with Canu v1.8 [101] with the option genome-Size=28m. The assembly obtained with Canu was polished twice with Illumina reads [31] using Pilon v1.24 [102] with the default settings. For each round of polishing, Bowtie2 [103] was used to align the Illumina trimmed reads with the assembly for polishing with a maximum length (-X) of 1000 bp. The Illumina reads were first trimmed with Trimmomatic v0.36 [104] with the following options: ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:10 TRAILING:10 SLIDINGWINDOW:5:10 MINLEN:50. Redundant contigs were removed based on a self-alignment using NUCmer v3.1 [105]. For other parental strains (LCP06037, LCP06039, LCP06043, LCP06059), available Illumina reads [31] were trimmed with Trimmomatic v0.36 [104] the same way as previously described, and assembled with SOAPdenovo2 v2.04 [106] using the corresponding default size 23 k-mer length. The assemblies obtained were then scaffolded using Ragout v2.3 with both the LCP06136 and the LCP06133 genomes as references [107]. For the 11 offspring genomes, the raw Nanopore reads were trimmed and assembled both with Canu v1.8 [101] and with Flye v2.9 [108], and then merged using Quickmerge [109]. These obtained assemblies were polished once with Illumina reads of the parental strain LCP06136 and LCP06133 [31] in the same way as for the LCP06136 parental genome. The obtained assemblies were scaffolded using Ragtag v2.1.0
[110]. We evaluated the quality of the final assemblies with Quast v5.0.2 [111] and their completeness with BUSCO v5.3.2 [112] using the eurotiales_odb10 lineage dataset. We checked the assemblies for contaminations with the GX cross-species aligner [113], part of NCBI's Foreign Contamination Screen tool suite, using the database shipped with version 0.5.0 (64,164 genome assemblies).

Transposable elements (TEs) were annotated using the REPET package (https://urgi.versailles.inra.fr/Tools/REPET). Briefly, the TEdenovo pipeline [114] was used to detect repeated elements in LCP06136 and LCP06133 (https://doi.org/10.57745/SIP7CH) genomes, and to provide consensus sequences. These consensus sequences were classified with PASTEC v1.3 [115], based on the Wicker hierarchical TE classification system [116] and manually curated. The resulting bulk library of consensus sequences was then used to annotate TE copies in the six whole genomes (namely, LCP06136, LCP06133, LCP06037, LCP06039, LCP06043 and LCP06059) using the TEannot pipeline [117]. Repeat-induced point mutations (RIP) in assemblies were detected using the web-based tool The RIPper [118].

All final assemblies were formatted (contigs being renamed and ordered in descending order of size), and repeats of the LCP06133 and LCP06136 parental genomes were masked with RepeatMasker v4.1.2 [119] after de novo repeat detection using RepeatModeler v2.0.2 [120]. We performed gene annotation on the masked assemblies of the parental strain LCP06136 using the Funnannotate v1.8.9 pipeline [121], with Braker v2.1.6 [122–125], which uses a combination of the *ab initio* gene predictors Augustus and GeneMark-ET [126] with NCBI blast [127] and blast+ [128]. We ran BRAKER twice, first using the BUSCO dataset of proteins eurotiales_odb10 [112] with the ProtHint pipeline [129–133], and then using the RNA-seq dataset from [134] mapped using STAR v2.5.4b [135] with the default settings. We combined the results of the two BRAKER runs with TSEBRA [136].

We studied synteny between assemblies by investigating collinearity of one-to-one whole genome alignments, using nucmer v3.1 [105], with 500 pb as the minimum match length. We visualized synteny between genomes by plotting one-to-one ortholog links using Circos v0.69-6 [137].

To determine the genetic divergence between parental strains, ANI (average nucleotide identity) values were estimated using FastANI v1.34 [138], based on the assembled genomes, for each pair of parental strains. The resulting ANI values are provided in Table 1.

We looked for potential *Starship* elements in LCP06133 (non-Roquefort) and LCP06136 (Roquefort) with the “annotate” and “sketch” modules of the Starfish pipeline [51], using the alignment profiles and protein sequences of tyrosine recombinases (referred as “captain”) and auxiliary proteins included in the pipeline.

## Phenotyping

We selected four traits for their applied interest in cheesemaking that are easily measurable with precision and at high throughput, for which we assessed the phenotypes in all the isolated offspring obtained for the five analyzed crosses. The growth rate was measured using a ScanStation (Intersciences, Mourjou, France), which is a temperature-controlled incubator. We set the ScanStation at 22.5°C and it took pictures regularly, in batches of 100 Petri dishes, i.e., the full capacity of the ScanStation. We spread 5 µL of standardized spore suspension (250 spores/inoculation) on each Petri dish, containing malt medium (Biokar BK045 HA) at 15g.L⁻¹ in source water (Cristalline, Chelles, France), sterilized at 120°C during 15 minutes. The ScanStation took pictures of the Petri dishes right side up on a black background, for five days every 30 minutes, i.e., 240 measures per strain (Fig 5D). The ScanStation estimated the colony area and diameter for each picture, by counting the number of pixels corresponding to the fungal colony. Because the colony diametral growth was near linear when plotted against days (S6 Fig), we used, as a measure of growth rate, the slope of the least square linear regression computed with the *lm* function of the base package R software v4.3.2 (http://www.r-project.org/). The coefficients of determination (r²), measuring the fit of the data to a line, are presented in the S6 Table and the linear regression of an offspring with a median r² in the S6 Fig. To ensure the robustness of the linear fit, strains with coefficients of determination below 0.98 were not considered in the QTL analysis of this trait (21 out of 781 offspring phenotyped, i.e., 2.6% of strains filtered out).
The colony color was measured after growth on raw sheep milk powder medium (Biocoop, Paris, France) at 21% in source water, sterilized at 105°C for five minutes, and mixed with 1.7% agar in source water, sterilized at 121°C for 15 minutes. After 13 days of growth, a picture was taken by the ScanStation under standardized conditions of light and on a white background (Fig 5F). We recorded for each strain the colony color decomposition using RGB (red, green and blue) and HSB (hue, saturation, and brightness) with ImageJ v1.52n [139] (S7 Table). In the RGB system, the levels of red, green and blue are represented by a range of integers from 0 to 255 (256 levels for each color). In the HSB system, hue ranges from 0 to 360 (0 is red, 120 is green, 240 is blue), saturation and brightness from 0 to 100% (0% saturated is neutral grey, 100% saturated is the full color; 0% brightness is black, 50% brightness is normal, 100% brightness is white). In order to obtain uncorrelated components from the RGB data, we divided each value by the sum of the three components. For each of these parameters, we recorded average values across all the colony pixels for each strain. For QTL detection, we chose a *posteriori*, among these color parameters, those whose distributions and heritabilities were optimal for QTL analysis (i.e., largest and uncorrelated variance); we thus retained the RGB proportion, as well as HSB.

The lipolytic and proteolytic activities of *P. roqueforti* strains were measured *in vitro* following [31]: for each strain, 50 µl of standardized spore suspensions (ca. 2,500 spores/inoculation) were inoculated at the top of a test tube containing agar and tributyrin for lipolytic activity measure (10mL·L⁻¹, ACROS Organics) or semi-skimmed cow milk for the proteolytic activity measure (40g·L⁻¹, Casino). The lipolytic and proteolytic activities were estimated by the degree of compound degradation, which changes the media from opaque to translucent. We measured the distance between the initial mark and the hydrolyzed, translucent front, after 7, 14, 21 and 28 days of growth at 20°C in the dark (Fig 5E). As lipolysis curves were nearly linear when plotted against days, we estimated the lipolysis rate with a linear regression computed with the R base package function *lm*. Proteolysis profiles appeared less linear but with no clear sigmoid nor other regression patterns, so we also chose to fit linear models for more simplicity. The coefficients of determination (r²), measuring the fit of the data to a line, are presented in the S6 Table and an example of a linear regression of an offspring with a median r² in the S6 Fig.

For extrolite production, we grew fungal cultures in 24-well sterile microplates containing 2mL of yeast extract sucrose (YES) agar medium buffered at pH 4.5 with phosphate-citrate buffer and characterized by a high C/N ratio to favor extrolite production as previously described [140]. For each strain, 1µL of a calibrated spore suspension (10⁶ spores.mL⁻¹) prepared from a 7-day culture was inoculated in the center of the well. Two replicates per strain were performed for extrolite analyses. The plates were incubated at 25°C in the dark for seven days and then stored at -20°C until extrolite analysis. For extrolite extractions, we used an optimized "high-throughput" extraction method [38,45]. Briefly, 2 g-aliquots (the entire YES culture obtained from a well) were homogenized after thawing samples with a sterile flat spatula then 12.5mL of acetonitrile (ACN) supplemented with 0.1% formic acid (v/v) was added, samples were vortexed for 30 sec followed by 15 min sonication. Extracts were again vortexed before being centrifuged for 10 min at 5000g at 4°C. The supernatants were directly collected and filtered through 0.45 µm PTFE membrane filters (GE Healthcare Life Sciences, UK) into amber vials and stored at -20°C until analyses. Extrolite detection and quantification were performed using an Agilent 6530 Accurate-Mass Quadropole Time-of-Flight (Q-TOF) LC/MS system equipped with a Binary pump 1260 and degasser, well plate autosampler set to 10°C and a thermostatted column compartment. Filtered 2µL aliquots were injected into a ZOR-BAX Extend C-18 column (2.1x50 mm and 1.8 µm, 600 bar) maintained at 35°C with a flow rate set to 0.3mL.min⁻¹. The mobile phase A contained milli-Q water + 0.1% formic acid (v/v) and 0.1% ammonium formate (v/v) while mobile phase B was ACN + 0.1% formic acid. Mobile phase B was maintained at 10% for 4 min, followed by a gradient from 10 to 100% for 18 min. Then, mobile phase B was maintained at 100% for 5 min before a 5-min post-time. Samples were ionized in positive (ESI+) electrospray ionization mode in the mass spectrometer with the following parameters: capillary voltage 4 kV, source temperature 325°C, nebulizer pressure 50 psig, drying gas 12 L.min⁻¹, ion range 100–1000 m/z. Target extrolite characteristics used for quantifications are given in S8 Table and included commercially available extrolites produced by Penicillium species, namely andrastin A, eremofortins A & B, (iso)-fumigaclavin A, mycophenolic acid and roquefortine C.
<figure><img src="image_25.png" /></figure>

Andrastin A, eremofortins A & B and (iso)-fumigaclavin A standards were obtained from Bioviotica (Goettingen, Germany), while all others were from Sigma-Aldrich (St Louis, MO, USA). All stock solutions were prepared in dimethyl sulfoxide (DMSO) at 1 mg.mL⁻¹ in amber vials. As the PR toxin was not commercially available, a previously produced purified stock solution [34] without known concentration was used to ensure MPA [mycophenolic acid) and PR toxin were separated (as they have the same mass) as well as linearity of PR toxin quantification. For these analyses, metabolite identification was performed using both the mean retention time ± 1 min and the corresponding ions listed in S8 Table. We used a matrix-matched calibration curve (r²>0.99 for all extrolites except 2>0.96) to confirm linearity the relation between signal area and extrolite concentration with final concentrations ranging from 1 to 10000 ng.mL⁻¹ according to the target metabolite and method performance was carried out as previously described [34].

## Statistical analysis on phenotypes

For phenotype distributions, we performed tests of normality with the *shapiro.test* function of the base package of R software v4.3.2 (http://www.r-project.org/). We ran tests for unimodality with the *dip.test* function of the Diptomtest package v0.75-7 [141] and tests for bimodality with the *bimodalityCoefficient* function of the Mousetrap package v3.2.1 [142] with a threshold of 0.555 according to [143]. We performed the trimodality and amodality tests with the *is.trimodal* and *is.amodal* functions, respectively, of the LaplaceS demon v16.1.6 package [144]. We computed the average heterosis as ((F1 - MP)/MP) * 100, F1 and MP being the trait mean values of the progeny and the parents, respectively.

We computed the matrix of Pearson and Spearman correlation coefficients with the R stats package v4.3.2 (http://www.r-project.org/). We produced the violin plots using vioplot v0.4.0 [145], ggplot2 v3.4.3 [146] and reshape2 v1.4.4 [147] R packages. We performed the principal component analysis (PCA) using the function *PCA* from the R base package, assigning trait means to missing values. Plots were drawn using FactoMineR v2.8 [148], factoextra v1.0.7 [149] and corrplot v0.92 [150] packages.

## Phenotype heritability

Narrow-sense heritability (h²) for each cross x phenotype combination was computed using a linear mixed model approach implemented in the "rrBLUP" R package [151,152]. Phenotypic values were modeled in a mixed model, using "mixedsolve" as a function of a genetic relatedness matrix derived from the genotypes matrix with the function "A.mat" of the "rrBLUP" package, to estimate additive genetic variance (σ<sub>A</sub>) and residual variance (σ<sub>e</sub>). The narrow-sense heritability was calculated as the proportion of the total phenotypic variance explained by additive genetic effects (h² = σ<sub>A</sub> / (σ<sub>A</sub> + σ<sub>e</sub>)). Narrow-sense heritability estimates for each cross x phenotype combination were further refined using the delete-one jackknife method to assess the stability of the estimates and compute standard errors. For each trait within each cross, individual segregants were iteratively excluded, and heritability was recalculated to derive standard errors.

Extrolite production was determined for three crosses (SxN, RxN and Rxl), with two replicates for each strain; we therefore estimated broad-sense heritability (H²) for these phenotypes for each cross based on variability among replicates, in addition to narrow-sense heritability as for the other phenotypes. Broad sense heritability was determined by performing a one-way ANOVA with the *lm* function of the stats v4.2.1 package in R.

## Genotyping based on indels

We genotyped 389 offspring from the cross RxN between LCP06136 (Roquefort; MAT1-1; https://www.ebi.ac.uk/SAMEA103939766) and LCP06133 (non-Roquefort; MAT1-2; SAMEA103939763) with 200 markers with a polymorphism of amplification size between the two parental strains, due to the presence of indels (S9 Table). Primers were designed (i) to yield an amplicon size between 150-410 bp and a difference ranging from six to 20 bp between alleles, (ii) to be compatible in multiplex and (iii) to be evenly distributed along the assembly of the non-Roquefort reference genome available at the beginning of the study (FM164; EMBL accession numbers HG792015 to HG792062, [153]. Indels were detected
<figure><img src="image_26.png" /></figure>

after a whole assembly alignment (blast v2.9.0; [127] of the two parental strains with the gap and extension penalties set to 1. After adding 250 base pairs upstream and downstream from the hit position, we performed a multi-alignment using MAFFT [154] with the whole genomes of the parental strains of the other crosses. We only kept a consensus sequence generated using consambig [96] and designed primers using Primer3 [155].

To reduce the cost of fluorescence primers, we only used four universal fluorescent primers (M13(-21), D8S1132f, D12S1090f, DYS437f), labeled with different fluorescent dyes (FAM, ATTO550, ATTO565, YAKIMA Yellow; Eurofins, Ebersberg, Germany). The specific primers attached during PCR to universal fluorescent primers thanks to a tail added to each specific primer matching an edge of the universal primers [156,157]. For each locus, the PCR amplification thus required three primers: the forward fluorescent universal primer, the forward specific primer with a 5' tail corresponding to the chosen universal primer and the regular specific reverse primer. The 200 markers were genotyped in 25 multiplex groups of 8 locus organized according to their amplicon size and their color fluorescent dye (S9 Table). For each multiplex, two separate sets of quadruplex PCR reactions were performed to minimize the interaction between primers and then pooled into octoplexes for the genotyping step. The PCR mix contained 0.16 µM of each of the fluorescent universal primer and of the reverse specified primer and 0.04 µM of the 5' tail forward primer in a final 15 µl reaction volume (2x QIAGEN Multiplex PCR Master Mix with 3 mM Mg²⁺, 10x primer mix with 1.6 or 0.4 µM of each primer, 3 µl of DNA diluted 50 fold) using the same PCR cycle conditions as for the microsatellites loci previously described in the recombinant offspring detection paragraph, except that it started with a hybridization temperature of 62°C instead of 60°C. Genotyping and analysis were also done as previously described in the recombinant offspring detection paragraph.

## Genotyping by sequencing

For all five crosses, we genotyped the obtained progeny (1073 offspring in total) at the AGAP CIRAD platform with genotyping-by-sequencing using the ApeKI enzyme. DNAs were purified after digestion on QIAquick columns (Qiagen, Les Ulis, France) and their qualities were controlled using a TapeStation instrument (Agilent Technologies, Les Ulis, France). Sequencing was performed on an Illumina Hiseq 4000 (2x150 bp). The following numbers of recombinant offspring were genotyped by GBS: 384 for the cross RxN, 157 for the cross RxB, 185 for the cross RxA, 176 for the cross SxN and 171 for the cross LxN.

Demultiplexed reads were mapped to the LCP06136 high-quality genome for the crosses RxA and RxB and to the LCP06133 high-quality genome for the crosses RxN, LxN and SxN using Bowtie2 v2.3.4.1 [103,158]. In Bowtie2, we set the maximum fragment length (-X) to 1000 and used the preset "very-sensitive-local". We used SAMtools v1.7 [159] to sort and filter out duplicate reads and reads with a mapping quality score above ten for SNP calling. Single nucleotide polymorphisms (SNPs) were called using the GATK v4.1.2.0 Haplotype Caller [160], generating a gVCF file per strain with option ploidy 2 to detect potential duplicated regions in progenies, as we had found indels with two alleles in some offspring. GVCFs were combined using GATK CombineGVCFs, genotypes with GATK GenotypeGVCFs, and biallelic SNPs were selected after filtration using GATK SelectVariants. We filtered SNPs using GATK VariantFiltration and options QUAL <30, DP<5, QD<20.0, FS>60.0, MQ<40, 0.05<AF<0.95, SOR>3.0, MQRankSum<-12.5, ReadPosRankSum<-8.0, DP>5.0 and GQ<10. We generated the genotype matrix based on the VCF file with Bcftools v1.11 [159,161]. Markers with no data or with no differences between the two parental genomes were filtered out. To reconstruct offspring genotypes, we performed SNP calling in the diploid mode in order to detect double copies of translocated regions. Markers corresponding to the translocated regions were in fact called as heterozygous in some offspring; they were therefore placed in offspring at the two parental locations, which were each assigned the corresponding parental allele.

After performing crosses and isolating ascospores, we detected sectors of two slightly different morphotypes in subsequent cultures of the parental strain LCP06039 aiming at phenotype measurements. However, genetic data indicated that the cross was performed with a pure strain and not a mixture: (i) the indel genotypes of the two morphotypes were assessed and were identical, and (ii) the GBS genotypes of the two morphotypes were not typed directly, but we did not
<figure><img src="image_27.png" /></figure>

detect genome-wide segregation distortion in the progeny in GBS data for this cross, while a cross involving a mixture of strains with different genotypes should result in over-representation of the alleles by the pure parent; on the contrary, the segregation distortions detected in some genomic regions corresponded to over-representations of the strain showing two morphotypes rather than its under-representation. We therefore considered that a single strain contributed to the cross and that the two morphotypes were the result of plasticity or a few somatic mutations in a single strain. We therefore built the genetic map and ran a QTL analysis for the LxN cross. For analyses involving parental traits, and therefore potentially affected by the presence of two morphotypes (i.e., for heritability estimates and transgression detection) we computed estimates using all three possible values for the LCP06039 parent strain phenotype: the means of the trait values for the two morphotypes or the value of one of the two morphotypes. This did not affect conclusions as the two morphotypes were actually very close for all trait measures.

## Experiment on the stability of genotypes across replication

As we identified a genomic region with two alleles present in multiple offspring in the RxN cross with associated QTLs, we wanted to check whether this particular genotype was stable during strain culture and across multiple replication steps. We therefore cultivated, on malt agar, the two parental strains from the RxN cross and the 11 offspring for which genomes were sequenced with Nanopore, six with alleles of both parents for this genomic region, four with the LCP06136 parental allele and one with the LCP06133 allele. We cultivated the strains for one week on malt agar and replated them on a new plate by transplanting some conidia from the edge of the colonies. We performed 19 replication steps. In the end, we re-genotyped the lines with the six indel markers present in the rearranged region.

## Linkage map and QTL detection

Obtained offspring phenotypes are presented in S10 Table and offspring genotypes are presented in S11-S15 Tables for SxN, LxN, RxN, RxS and RgL crosses respectively. The genetic map and statistical tests for QTL association were performed using the R software v3.5.1 (http://www.r-project.org/) with the ASMap v1.0-5 package [162] for the genetic map and qtl v1.60 package [163] for statistical tests. For the genetic map, we used unique segregating markers with less than 10% missing data and markers with segregation distortion with significant Chi Square p-values adjusted with Bonferroni corrections. Individuals with more than 30% of missing data were filtered out (Table 1), concerning up to 14% of the offspring of a cross. Linkage groups and marker ordering were determined with the *mstmap* function of AsMap package, with parameter "bychr" set as FALSE and a p.value of 10⁻¹². Markers with less than 20% missing data and less than 95% of one parental allele were then reintegrated in linkage groups with the *pull-cross* function with "max.rf" and "min.lod" parameters set at 0.1 and 20, respectively. We reconstructed marker order in linkage groups with the *mstmap* function, the parameter "bychr" set as TRUE, and a p.value < 10⁻¹². We estimated genetic distances between markers with the AsMap *quickEst* function and the Kosambi's mapping function. We compared the obtained order to the reference genome order using the *compareorder* function of the qtl package v1.6 with error probability set at 0.005. We changed the order only if it improved the LOD score and minimized the linkage group length. We represented the relationships between chromosome physical positions and genetic maps using ALLMAPS in JCVI utility libraries v1.3.6 [164], adapted to include genetic features.

We performed single QTL detection using the *scanone* function of the qtl package with the threshold at 5% being estimated with 1000 permutations. We used the *addqtl* function of the qtl package to mask the first set of QTLs and thus detect additional QTLs that may have been hidden by the main ones. We estimated confidence intervals for QTL using the *lodint* function of the qtl package with default parameters. We identified QTL interactions using the *scantwo* function with the Haley-Knott method using a threshold at 5% estimated by 1,000 permutations. We performed analyses of variance (ANOVAs) using a linear model with the *fitqtl* function of the qtl package, keeping interactions between loci only when they were significant in both *scantwo* and the full model AVOVAs. We extended QTL confidence intervals including the first
or last marker of a linkage group to the start or end of the corresponding chromosome. For comparisons of QTL interval positions between crosses, we transposed the intervals determined with the positions on the Roquefort parent genomes (for RxL and RxB crosses) onto the non-Roquefort genome. To do this, we extracted genome sequences using bedtools getfasta v2.30.0 [165] and mapped the extracted sequences onto the non-Roquefort parent genome using minimap2 v2.17-r941 [166]. The representation of the aligned QTL intervals was constructed using the karyoplotR v1.24.0 package [167].

## Supporting information

**S1 Fig.** Focus on the L4 translocated region in the *Penicillium roqueforti* LCP06037 strain, the same pattern being observed for the L1 region in the LCP06039 strain. The L4 region is subdivided into two subregions, L4.1 in red and L4.2 in blue, with different relative positions, suggesting a circular intermediate for the translocation. The cut sites for the proposed circular intermediate are indicated by purple and green arrows, respectively, and numbered 1 and 2. The figure shows the synteny of the L4 translocated region present in chromosome 1 in the LCP06037 lumber parental strain and in chromosome 4 in other parents and the proposed circular intermediate.

(TIF)

**S2 Fig.** Distributions of phenotypes measured in the five progenies: (A) diametral growth in pixel per hour, (B) lipolysis rate in mm per day, (C) proteolysis rate in mm per day, (D) relative red, (E) green, (F) blue, (G) hue, (H) saturation and (I) brightness; progenies are represented along the x axis, RxL in red, RxB in beige, RxN in yellow, LxN in turquoise and SxN in pink. The black point and line represent the mean and standard deviation, respectively, in progenies. The coloured points represent the parental values (lumber in green, non-Roquefort in blue, Roquefort in purple, silage in orange).

(TIF)

**S3 Fig.** Distributions of extrolite production in three progenies: (A) PR toxin, (B) mycophenolic acid, (C) andras tin A, (D) eremofortin A, (E) eremofortin B, (F) (iso)-fumigaclavine A and (G) roquefortine C; progenies are repre- sented along the x axis: RxL in red, RxN in yellow and SxN in pink. The coloured points represent the parental values (lumber in green, non-Roquefort in blue, Roquefort in purple, silage in orange). All units are in mL.L⁻¹ except the PR toxin, which is expressed in arbitrary units.

(TIF)

**S4 Fig.** Biplot of a principal component analysis with the two first components (A) and the third and fourth (B). For every cross the progeny scores are displayed in a different color (LxN in red, RxL in brown, RxN in yellow, RxB in turquoise, SxN in pink) and with two density levels (from 0.05 to 0.10 in light color and above 0.10 in dark color). Large dots represent the barycentre of each progeny's scores using the same color code (the red one is under the yellow one). The 4-fold loadings of the phenotypes are represented by the black arrows. Small dots show the scores of the parental strains from different populations (lumber in green, non-Roquefort in blue, Roquefort in purple and silage in orange) in the labels show the cross they are involved in. The percentage in the axis labels shows the total variance explained by the component.

(TIF)

**S5 Fig.** Number of phenotype classes (i.e., lipolysis, proteolysis, color and extrolite production) presenting a QTL along the genome positions. The crosses are represented one per line, with the cross ID indicated on the right, the first parent carrying the MAT1-1 mating type; the bottom line represents a pool of all crosses. The x-axis represents genomic physical positions (Mb). At the bottom, the four chromosomes of the reference genome (LCP06133) are represented by rectangles depicting genomic features: large horizontally transferred regions specific to the reference genome,
in blue, and translocated regions in the reference genome, in red. The large vertical empty red rectangles indicate pleio-
tropic QTL regions, with effects on multiple traits in at least one cross.

(TIF)

S6 Fig. Example of phenotype measures and parameter estimation, for an offspring having a median r² value for the least squares regression fit: (A) growth dynamic, i.e., colony diameter (in pixel numbers on pictures) as a function of time (hours), (B) lipolysis dynamics and (C) proteolysis dynamics, B and C representing the distance of the lysis front (i.e., characterized by medium fading) in mm per day. The linear regression lines are shown in blue, and their equations and coefficients of determination are given at the top.

(TIF)

S1 Table. Result of normality (Wilk), unimodality, bimodality, trimodality and amodality test for each phenotype x
cross combination.

(CSV)

S2 Table. Heterosis and heterobeltosis values and significance for each phenotype x cross combination.

(CSV)

S3 Table. Narrow sense and broad sense heritability for each tested phenotype.

(CSV)

S4 Table. List of identified QTLs for each cross and phenotype.

(CSV)

S5 Table. List of genome and GBS data used in this article.

(CSV)

S6 Table. Determination coefficient (Rsq) of growth, lipolysis and proteolysis linear regression for parents and offsprings.

(CSV)

S7 Table. Script used for color phenotype determination.

(CSV)

S8 Table. Method performance characteristics for metabolite quantification in YES medium.

(CSV)

S9 Table. INDEL list, primer sequence and amplicon size.

(CSV)

S10 Table. Phenotypes of parents and offsprings.

(CSV)

S11 Table. Genotypes of offspring from the SxN cross. S: LCP06059 allele; N: LCP06133 allele.

(CSV)

S12 Table. Genotypes of offspring from the LxN cross. W: LCP06039 allele (L); N: LCP06133 allele.

(CSV)

S13 Table. Genotype of offspring from the RxN cross. R: LCP06136 allele; N: LCP06133 allele.

(CSV)
**S14 Table. Genotypes of offspring from the RxB cross.** R: LCP06136 allele; S: LCP06043 allele.
(CSV)

**S15 Table. Genotypes of offspring from the RzL cross.** R: LCP06136 allele; W: LCP06037 allele (L).
(CSV)

**Acknowledgments**

TC and ECr acknowledge the E3GP3 network for facilitating collaborations and skills transfer from MF laboratory. We acknowledge the staff at the GENTYANE genotyping platform (INRAe GDEC, Clermont-Ferrand, France) for microsatellite and indel genotyping and the AGAP CIRAD platform for GBS sequencing. We gratefully acknowledge Intersciences for the loan and set-up of the ScanStation (https://www.interscience.com/en/products/real-time-incubator-and-colyony-counter/). We thank Jeanne Ropars for her help in obtaining the progeny and Jean-Philippe Vernadet for help with bioinformatic analyses.

**Author contributions**

**Conceptualization:** Thibault Caron, Daniel Roueyre, Michel Place, Christophe Chassard, Antoine Branca, Tatiana Giraud.

**Data curation:** Thibault Caron, Ewen Crequer.

**Formal analysis:** Thibault Caron, Ewen Crequer, Ricardo C. Rodríguez de la Vega.

**Funding acquisition:** Daniel Roueyre, Michel Place, Emmanuel Coton, Tatiana Giraud.

**Investigation:** Thibault Caron, Ewen Crequer, Mélanie Le Piver, Stéphanie Le Prieur, Sammy Brunel, Alodie Snirc, Gwennina Cueff, Adeline Simon.

**Methodology:** Thibault Caron, Ewen Crequer.

**Project administration:** Tatiana Giraud.

**Software:** Ricardo C. Rodríguez de la Vega.

**Supervision:** Daniel Roueyre, Michel Place, Christophe Chassard, Monika Coton, Emmanuel Coton, Marie Foulongne-Oriol, Antoine Branca, Tatiana Giraud.

**Writing – original draft:** Thibault Caron, Ewen Crequer, Tatiana Giraud.

**Writing – review & editing:** Thibault Caron, Ewen Crequer, Mélanie Le Piver, Stéphanie Le Prieur, Alodie Snirc, Gwennina Cueff, Daniel Roueyre, Michel Place, Christophe Chassard, Adeline Simon, Ricardo C. Rodríguez de la Vega, Monika Coton, Emmanuel Coton, Marie Foulongne-Oriol, Antoine Branca, Tatiana Giraud.

**References**

1. Falconer DS, Mackay TFC. Introduction to Quantitative Genetics. Pearson Education Limited; 1995.

2. Adamczyk K, Pokorska J, Makulska J, Earley B, Mazurek M. Genetic analysis and evaluation of behavioural traits in cattle. Livest Sci. 2013;154:1–12.

3. Kumar J, Gupta DS, Gupta S, Dubey S, Gupta P, Kumar S. Quantitative trait loci from identification to exploitation for crop improvement. Plant Cell Rep. 2017;36(8):1187–213. https://doi.org/10.1007/s00299-017-2127-y PMID: 28352970

4. Silva MVB, dos Santos DJA, Boison SA, Utsunomiya ATH, Carmo AS, Sonstegard TS, et al. The development of genomics applied to dairy breeding. Livest Sci. 2014;166(1):66–75.

5. Andargie M, Pasquet RS, Gowda BS, Muluvi GM, Timko MP. Molecular mapping of QTLs for domestication-related traits in cowpea (V. unguiculata (L.) Walp.). Euphytica. 2014;200(3):401–12.

6. Bachlava E, Tang S, Pizarro G, Schuppert GF, Brunick RK, Draeger D, et al. Pleiotropy of the branching locus (B) masks linked and unlinked quantitative trait loci affecting seed traits in sunflower. Theor Appl Genet. 2010;120(4):829–42. https://doi.org/10.1007/s00122-009-1212-1 PMID: 19921140
7. Bomblies K, Doebley JF. Pleiotropic effects of the duplicate maize FLORICAULA/LEAFY genes zfl1 and zfl2 on traits under selection during maize domestication. Genetics. 2006;172(1):519–31. https://doi.org/10.1534/genetics.105.048595 PMID: 16204211

8. Johnsson M, Rubin C-J, Höglund A, Sahlqvist A-S, Jonsson KB, Kerje S, et al. The role of pleiotropy and linkage in genes affecting a sexual orna-
ment and bone allocation in the chicken. Mol Ecol. 2014 May;23(9):2275–86.

9. Kantar MB, Nashoba AR, Anderson JE, Blackman BK, Riesenberg LH. The Genetics and Genomics of Plant Domestication. BioScience. 2017;67(11):971–82.

10. Martínez-Ainsworth NE, Tenaillon MI. Superheroes and masterminds of plant domestication. C R Biol. 2016;339(7–8):268–73. https://doi.org/10.1016/j.crvi.2016.05.005 PMID: 27317057

11. Somta P, Chen J, Yimram T, Yundaeng C, Yuan X, Tomooka N, et al. QTL Mapping for Agronomic and Adaptive Traits Confirmed Pleiotropic Effect of mog Gene in Black Gram (Vigna mungo (L.) Hepper). Front Genet. 2020; 11:635.

12. Sweeney M, McCouch S. The Complex History of the Domestication of Rice. Ann Bot. 2007;100(5):951-7.

13. Wang J, Liao X, Li Y, Zhou R, Yang X, Gao L, et al. Fine mapping a domestication-related QTL for spike-related traits in a synthetic wheat. Genome. 2010;53(10):798–804. https://doi.org/10.1139/g10-066 PMID: 20962886

14. Wirén A, Jensen P. A Growth QTL on Chicken Chromosome 1 Affects Emotionality and Sociality. Behav Genet. 2011;41(2):303-11.

15. Wright D. The Genetic Architecture of Domestication in Animals. Bioinform Biol Insights. 2015;9(Suppl 4):BBI.S28902. https://doi.org/10.4137/BBI.S28902 PMID: 26512200

16. Wright D, Rubin CJ, Martinez Barrio A, Schütz K, Kerje S, Brändström H, et al. The genetic architecture of domestication in the chicken: effects of pleiotropy and linkage. Mol Ecol. 2010;19(23):5140–56.

17. Eder M, Sanchez I, Brice C, Camarasa C, Legras JL, Dequin S. Qtl mapping of volatile compound production in *Saccharomyces cerevisiae* during alcoholic fermentation. BMC Genomics. 2018;19:166.

18. Kessi-Pérez EI, Molinet J, Martínez C. Disentangling the genetic bases of Saccharomyces cerevisiae nitrogen consumption and adaptation to low nitrogen environments in wine fermentation. Biol Res. 2020;53(1):1–10.

19. Liti G, Louis EJ. Advances in quantitative trait analysis in yeast. PLoS Genet. 2012;8(8):e1002912. https://doi.org/10.1371/journal.pgen.1002912 PMID: 22916041

20. Nguyen BAN, Lawrence KR, Rego-Costa A, Gopalakrishnan S, Temko D, Michor F, et al. Barcoded bulk QTL mapping reveals highly polygenic and epistatic architecture of complex traits in yeast. eLife. 2022;11:e73983.

21. Steinmetz LM, Sinha H, Richards DR, Spiegelman JI, Oefner PJ, McCusker JH, et al. Dissecting the architecture of a quantitative trait locus in yeast. Nature. 2002;416(6878):326-30.

22. Swinnen S, Thevelein JM, Nevoigt E. Genetic mapping of quantitative phenotypic traits in Saccharomyces cerevisiae. FEMS Yeast Res. 2012;12(2):215-27. https://doi.org/10.1111/j.1567-1364.2011.00777.x PMID: 22150948

23. Wilkening S, Lin G, Fritsch ES, Tekkedil MM, Anders S, Kuehn R, et al. An evaluation of high-throughput approaches to QTL mapping in Saccharomyces cerevisiae. Genetics. 2014;196(3):853–65. https://doi.org/10.1534/genetics.113.160291 PMID: 24374355

24. Todd RT, Forche A, Selmeki A. Ploidy Variation in Fungi - Polyploidy, Aneuploidy, and Genome Evolution. Microbiol Spectr. 2017;5(4):1-31.

25. Zimmer A, Durand C, Loira N, Durrens P, Sherman DJ, Marullo P. QTL dissection of lag phase in wine fermentation reveals a new translocation responsible for Saccharomyces cerevisiae adaptation to sulfite. PLoS ONE. 2014;9(1):e86298. https://doi.org/10.1371/journal.pone.0086298 PMID: 24489712

26. Foulongne-Oriol M, Rodier A, Rousseau T, Savoie J-M. Quantitative trait locus mapping of yield-related components and oligogenic control of the cap color of the button mushroom, Agaricus bisporus. Appl Environ Microbiol. 2012;78(7):2422-34. https://doi.org/10.1128/AEM.07516-11 PMID: 22267676

27. Imbernon M, Callac P, Gasqui P, Kerrigan RW, Velcko AJ. BSN, the primary determinant of basidial spore number and reproductive mode in Agaricus bisporus, maps to chromosome I. Mycologia. 1996;88(5):749–61. https://doi.org/10.1080/00275514.1996.12026713

28. Cerning J, Gripon JC, Lamberet G, Lenoir J. Les activités biochimiques des Penicillium utilisés en fromagerie. Le Lait. 1987;67(1):3-39.

29. Collins YF, McSweeney PLH, Wilkinson MG. Lipolysis and free fatty acid catabolism in cheese: a review of current knowledge. Int Dairy J. 2003;13(11):841–66.

30. Moreau C. Le Penicillium roqueforti, morphologie, physiologie, intérêt en industrie fromagère, mycotoxines. Le Lait. 1980;60(595-596):254-71.

31. Dumas É, Feurtey A, Rodríguez de la Vega RC, Le Prieur S, Snirc A, Coton M, et al. Independent domestication events in the blue-cheese fungus Penicillium roqueforti. Mol Ecol. 2020;33913:451773.

32. Caron T, Piver ML, Péron A-C, Lieben P, Lavigne R, Brunel S, et al. Strong effect of Penicillium roqueforti populations on volatile and metabolic compounds responsible for aromas, flavor and texture in blue cheeses. Int J Food Microbiol. 2021;354:109174. https://doi.org/10.1016/j.ijfoodmicro.2021.109174 PMID: 34103155

33. Crequer E, Ropars J, Jany L, Caron T, Coton M, Snirc A, et al. A new cheese population in Penicillium roqueforti and adaptation of the five popula-
tions to their ecological niche. Evol Appl. 2023;1–20.
34. Gillot G, Jany J-L, Poirier E, Maillard M, Debaets S, Thierry A, et al. Functional diversity within the Penicillium roqueforti species. Int J Food Microbiol. 2017;241:141–50. https://doi.org/10.1016/j.ijfoodmicro.2016.10.001 PMID: 27771579

35. Ropars J, Caron T, Lo YC, Bennetot B, Giraud T. The domestication of Penicillium cheese fungi. Comptes Rendus - Biol. 2020;1.

36. Ropars J, Rodríguez de la Vega RC, López-Villavicencio M, Gouzy J, Sallet E, Dumas É, et al. Adaptive Horizontal Gene Transfers Between Multiple Cheese-Associated Fungi. Curr Biol. 2015;25(19):2562–9. https://doi.org/10.1016/j.cub.2015.08.025 PMID: 26412136

37. Gluck-Thaler E, Ralston T, Konkel Z, Ocampos CG, Ganeshan VD, Dorrance AE, et al. Giant Starship Elements Mobilize Accessory Genes in Fungal Genomes. Mol Biol Evol. 2022;39(5):msac109. https://doi.org/10.1093/molbev/msac109 PMID: 35588244

38. Lo Y, Bruxaux J, Rodríguez de la Vega RC, O'Donnell S, Snirc A, Coton M, et al. Domestication in dry-cured meat Penicillium fungi: Convergent specific phenotypes and horizontal gene transfers without strong genetic subdivision. Evol Appl. 2023;16(9):1637–60. https://doi.org/10.1111/eva.13591 PMID: 37752962

39. Hartmann FE, Rodríguez De La Vega RC, Demené A, Badet T, Vernadet JP, Rougemont Q, et al. An inversion polymorphism under balancing selection, involving giant mobile elements, in an invasive fungal pathogen. Wilson M, editor. Mol Biol Evol. 2025;msaf026.

40. Peck LD, Llewellyn T, Bennetot B, O'Donnell S, Nowell RW, Ryan MJ, et al. Horizontal transfers between fungal Fusarium species contributed to successive outbreaks of coffee wilt disease. PLOS Biol. 2024;22(12):e3002480. https://doi.org/10.1371/journal.pbio.3002480 PMID: 39637834

41. Scott PM. Toxins of Penicillium Species Used in Cheese Manufacture. J Food Prot. 1981;44(9):702–10. https://doi.org/10.4315/0362-028X-44.9.702 PMID: 30856725

42. Fontaine K, Hymery N, Lacroix MZ, Puel S, Puel O, Rigalma K, et al. Influence of intraspecific variability and abiotic factors on mycotoxin production in Penicillium roqueforti. Int J Food Microbiol. 2015;215:187–93. https://doi.org/10.1016/j.ijfoodmicro.2015.07.021 PMID: 26320771

43. Hymery N, Puel O, Tadrist S, Canlet C, Le Scouarnec H, Coton E, et al. Effect of PR toxin on THP1 and Caco-2 cells: an in vitro study. World Mycotoxin J. 2017;10(4):375–86. https://doi.org/10.3920/wmj2017.2196

44. Matas AJ, Smith JM, Skeans MA, Thompson B, Gustafson SK, Stewart DE, et al. OPTN/SRTR 2013 Annual Data Report: Kidney. Am J Transplant. 2013;15:1–34.

45. Crequer E, Coton E, Cueff G, Cristiansen JV, Frisvad JC, de la Vega RR, et al. Different metabolite profiles across Penicillium roqueforti populations associated with ecological niche specialisation and domestication. Microbiology. 2024 [cited 2024 Jan 24]. Available from: http://biorxiv.org/lookup/doi/10.1101/2024.01.12.575369

46. Ropars J, López-Villavicencio M, Dupont J, Snirc A, Gillot G, Coton M, et al. Induction of sexual reproduction and genetic diversity in the cheese fungus Penicillium roqueforti. Evol Appl. 2014;7(4):433–41. https://doi.org/10.1111/eva.12140 PMID: 24822078

47. Ropars J, Lo Y-C, Dumas É, Snirc A, Begerow D, Rollnik T, et al. Fertility depression among cheese-making Penicillium roqueforti strains suggests degeneration during domestication. Evolution. 2016;70(9):2099–109. https://doi.org/10.1111/evo.13015 PMID: 27470007

48. Coton E, Coton M, Hymery N, Jany J-L, Mounier J. Penicillium roqueforti: an overview of its genetics, physiology, metabolism and biotechnological applications. Fungal Biology Reviews. 2020;34(2):59–73. https://doi.org/10.1016/j.fbr.2020.03.001

49. Chávez R, Vaca I, García-Estrada C. Secondary Metabolites Produced by the Blue-Cheese Ripening Mold Penicillium roqueforti; Biosynthesis and Regulation Mechanisms. J Fungi. 2023;9(4):459.

50. Bucknell AH, McDonald MC. That's no moon, it's a Starship: Giant transposons driving fungal horizontal gene transfer. Mol Microbiol. 2023;mmi.15118.

51. Gluck-Thaler E, Vogan AA. Systematic identification of cargo-mobilizing genetic elements reveals new dimensions of eukaryotic diversity. Nucleic Acids Res. 2024;52(10):5496–513.

52. Galagan JE, Selker EU. RIP: the evolutionary cost of genome defense. Trends Genet. 2004;20(9):417–23. https://doi.org/10.1016/j.tig.2004.07.007 PMID: 15313550

53. Torrent C, Gil-Durán C, Rojas-Aedo JF, Medina E, Vaca I, Castro P, et al. Role of sfk1 gene in the filamentous fungus Penicillium roqueforti. Front Microbiol. 2017;8(DEC):1–11.

54. Rojas-Aedo JF, Gil-Durán C, Goity A, Vaca I, Levicán G, Larrondo LF. The developmental regulator Pcz1 affects the production of secondary metabolites in the filamentous fungus Penicillium roqueforti. Microbiol Res. 2018;212–213:67–74.

55. Marcano Y, Montanares M, Gil-Durán C, González K, Levicán G, Vaca I, et al. PrlæA Affects the Production of Roquefortine C, Mycophenolic Acid, and Andrastin A in Penicillium roqueforti, but It Has Little Impact on Asexual Development. J Fungi. 2023;9(10):954.

56. Marquina M, Lambea E, Carmona M, Sánchez-Marinas M, López-Aviles S, Ayte J, et al. A new negative feedback mechanism for MAPK pathway inactivation through Srk1 MAPKAP kinase. Sci Rep. 2022;12(1):19501. https://doi.org/10.1038/s41598-022-23970-8 PMID: 36376357

57. Gil-Durán C, Rojas-Aedo JF, Medina E, Vaca I, García-Rico RO, Villagrán S, et al. The *pcz1* Gene, which Encodes a Zn(II)2Cys6 Protein, Is Involved in the Control of Growth, Conidiation, and Conidial Germination in the Filamentous Fungus *Penicillium roqueforti*. Pöggeler S, editor. PLOS ONE. 2015;10(3):e0120740.

58. Gripon JC, Desmazeaud MJ, Le Bars D, Bergere JL. Role of Proteolytic Enzymes of Streptococcus lactis, Penicillium roqueforti, and Penicillium caseicolum during Cheese Ripening. J Dairy Sci. 1977;60(10):1532–8.

59. Otto SP, Lenormand T. Resolving the paradox of sex and recombination. Nat Rev Genet. 2002;3(4):252–61. https://doi.org/10.1038/nrg761 PMID: 11967550
60. Olarte RA, Horn BW, Dorner JW, Monacell JT, Singh R, Stone EA, et al. Effect of sexual recombination on population diversity in aflatoxin production by Aspergillus flavus and evidence for cryptic heterokaryosis. Mol Ecol. 2012;21(6):1453-76.

61. Peter J, De Chiara M, Friedrich A, Yue JX, Pflieger D, Bergström A, et al. Genome evolution across 1,011 Saccharomyces cerevisiae isolates. Nature. 2018;556(7701):339-44.

62. Jakobson CM, Jarosz DF. Molecular Origins of Complex Heritability in Natural Genotype-to-Phenotype Relationships. Cell Syst. 2019; 8(5):363-79. e3.

63. Lin Y-C, Liu H-H, Tseng M-N, Chang H-X. Heritability and gene functions associated with sclerotia formation of Rhizoctonia solani AG-7 using whole genome sequencing and genome-wide association study. Microb Genom. 2023; 9(3):mgen000948. https://doi.org/10.1099/mgen.0.000948 PMID: 36867092

64. Lehman JS, Oudemans PV. Variation and Heritability of Phenology in the Fungus Monilinia vaccinii-corymbosi on Blueberry. Phytopathology. 2000;90(4):390-5. https://doi.org/10.1094/PHYTO.2000.90.4.390 PMID: 18944589

65. Brem RB, Kruglyak L. The landscape of genetic complexity across 5,700 gene expression traits in yeast. Proc Natl Acad Sci U S A. 2005;102(5):1572-7. https://doi.org/10.1073/pnas.0408709102 PMID: 15659551

66. Nogami S, Ohya Y, Yvert G. Genetic Complexity and Quantitative Trait Loci Mapping of Yeast Morphological Traits. Snyder M, editor. PLoS Genet. 2007;3(2):e31.

67. Greetham D, Wimalasena TT, Leung K, Marvin ME, Chelandia Y, Hart AJ, et al. The Genetic Basis of Variation in Clean Lineages of Saccharomyces cerevisiae in Response to Stresses Encountered during Bioethanol Fermentations. Du C, editor. PLoS ONE. 2014; 9(8):e103233.

68. Thomson MJ, Edwards JD, Septiningsih EM, Harrington SE, McCouch SR. Substitution Mapping of dth1.1, a Flowering-Time Quantitative Trait Locus (QTL) Associated With Transgressive Variation in Rice, Reveals Multiple Sub-QTL. Genetics. 2006;172(4):2501-14.

69. Li J, Thomson M, McCouch SR. Fine mapping of a grain-weight quantitative trait locus in the pericentromeric region of rice chromosome 3. Genetics. 2004;168(4):2187-95. https://doi.org/10.1534/genetics.104.034165 PMID: 15611185

70. Hartman Y, Uwimana B, Hooftman DAP, Schranz ME, van de Wiel CCM, Smulders MJM, et al. Genomic and environmental selection patterns in two distinct lettuce crop-wild hybrid crosses. Evol Appl. 2013;6(4):569-84. https://doi.org/10.1111/eva.12043 PMID: 23789025

71. Liu H, Liu H, Zhou L, Lin Z. Genetic Architecture of domestication- and improvement-related traits using a population derived from Sorghum virgatum and Sorghum bicolor. Plant Sci. 2019;283:135-46. https://doi.org/10.1016/j.plantsci.2019.02.013 PMID: 31128683

72. Koszul R, Dujon B, Fischer G. Stability of large segmental duplications in the yeast genome. Genetics. 2006;172(4):2211-22. https://doi.org/10.1534/genetics.105.048058 PMID: 16489235

73. Perkins DD. Chromosome Rearrangements in and Other Filamentous Fungi. In: Advances in Genetics [Internet]. Elsevier; 1997 [cited 2025 Feb 17]. p. 239-398. Available from: https://linkinghub.elsevier.com/retrieve/pii/S0065266008603119

74. Villarreal P, O'Donnell S, Agier N, Muñoz-Guzman F, Benavides-Parra J, Urbina K, et al. Domestication signatures in the non-conventional yeast Lachancea cidri. Gibbons JG, editor. mSystems. 2024; 9(1):e01058-23.

75. Wang W, Chen L, Wang X, Duan J, Flynn RD, Wang Y, et al. A transposon-mediated reciprocal translocation promotes environmental adaptation but compromises domesticability of wild soybeans. New Phytol. 2021;232(4):1765-77. PMID: 34363228

76. Akpinar BA, Biyiklioglu S, Alptekin B, Havránková M, Vrána J, Doležel J, et al. Chromosome-based survey sequencing reveals the genome organization of wild wheat progenitor Triticum dicoccoides. Plant Biotechnol J. 2018;16(12):2077-87. https://doi.org/10.1111/pbi.12940 PMID: 29729062

77. Zhou X, Jellen EN, Murphy JP. Progenitor Germplasm of Domesticated Hexaploid Oat. Crop Sci. 1999;39(4):1208-14.

78. Maughan PJ, Jarvis DE, De La Cruz-Torres E, Jaggi KE, Warner HC, Marcheschi AK, et al. North American pitseed goosefoot (Chenopodium berlandieri) is a genetic resource to improve Andean quinoa (C. quinoa). Sci Rep. 2024;14(1):12345.

79. Lopez-Moreno H, Basurto-Garduño AC, Torres-Meraz MA, Diaz-Valenzuela E, Arellano-Arciniega S, Zalapa J, et al. Genetic analysis and QTL mapping of domestication-related traits in chili pepper (Capsicum annuum L.). Front Genet. 2023;14:1101401.

80. Kou Y, Liao Y, Toivainen T, Lv Y, Tian X, Emerson JJ, et al. Evolutionary Genomics of Structural Variation in Asian Rice (Oryza sativa) Domestication. M Purugganan, editor. Mol Biol Evol. 2020; 37(12):3507-24.

81. Zhou Y, Minio A, Massonnet M, Solares E, Lv Y, Beridze T, et al. The population genetics of structural variants in grapevine domestication. Nat Plants. 2019; 5(9):965-79.

82. Pelé A, Falque M, Trotoux G, Eber F, Nègre S, Gilet M, et al. Amplifying recombination genome-wide and reshaping crossover landscapes in Bras-sicas. PLOS Genet. 2017;13(5):e1006794. https://doi.org/10.1371/journal.pgen.1006794 PMID: 28493942

83. Chakraborty P, Pankajam AV, Lin G, Dutta A, Krishnaprasad GN, Tekkedil MM, et al. Modulating Crossover Frequency and Interference for Obligate Crossovers in Saccharomyces cerevisiae Meiosis. G3 GenesGenomesGenetics. 2017;7(5):1511-24.

84. Braun-Galleani S, Dias JA, Coughlan AY, Ryan AP, Byrne KP, Wolfe KH. Genomic diversity and meiotic recombination among isolates of the biotech yeast Komagataella phaffii (Pichia pastoris). Microb Cell Fact. 2019;18(1):211. https://doi.org/10.1186/s12934-019-1260-4 PMID: 31801527

85. Auxier B, Debets AJM, Stanford FA, Rhodes J, Becker FM, Reyes Marquez F, et al. The human fungal pathogen Aspergillus fumigatus can produce the highest known number of meiotic crossovers. PLOS Biol. 2023;21(9):e3002278. https://doi.org/10.1371/journal.pbio.3002278 PMID: 37708139
86. Gong W, Shen N, Zhang L, Bian Y, Xiao Y. Chromosome-Wide Characterization of Intragenic Crossover in Shiitake Mushroom, Lentinula edodes. *J Fungi*. 2021; 7(12):1076.

87. Gong W, Song X, Xie C, Zhou Y, Zhu Z, Xu C, et al. Landscape of meiotic crossovers in Hericium erinaceus. *Microbiol Res*. 2021;245:126692.

88. Stapley J, Feulner PGD, Johnston SE, Santure AW, Smadja CM. Variation in recombination frequency and distribution across eukaryotes: patterns and processes. *Philos Trans R Soc B Biol Sci*. 2017;372(1736):20160455.

89. Giraud TRJ. La domestication des champignons du fromage. La Recherche. 2022;569:46-9.

90. Ho P-W, Piamplantsant S, Gallone B, Del Cortona A, Peeters P-J, Reijbroek F, et al. Massive QTL analysis identifies pleiotropic genetic determinants for stress resistance, aroma formation, and ethanol, glycerol and isobutanol production in Saccharomyces cerevisiae. *Biotechnol Biofuels*. 2021;14(1):211. https://doi.org/10.1186/s13068-021-02059-w PMID: 34727964

91. Vion C, Peltier E, Bernard M, Muro M, Marullo P. Marker Assisted Selection of Malic-Consuming Saccharomyces cerevisiae Strains for Winemaking. Efficiency and Limits of a QTL's Driven Breeding Program. *J Fungi*. 2021;7(4):304.

92. Marullo P, Yvert G, Bely M, Aigle M, Dubourdieu D. Efficient use of DNA molecular markers to construct industrial yeast strains. *FEMS Yeast Res*. 2007;7(8):1295-306.

93. Okuda Y, Murakami S, Matsumoto T. Development of STS markers suitable for marker-assisted selection of sporeless trait in oyster mushroom, *Pleurotus pulmonarius*. *Breed Sci*. 2009;59(3):315-9.

94. Molinet J, Cubillos FA. Wild Yeast for the Future: Exploring the Use of Wild Strains for Wine and Beer Fermentation. *Front Genet*. 2020;11:589350.

95. Zhang F, Batley J. Exploring the application of wild species for crop improvement in a changing climate. *Curr Opin Plant Biol*. 2020;56:218-22. https://doi.org/10.1016/j.pbi.2019.12.013 PMID: 32029361

96. Renzi JP, Coyne CJ, Berger J, Von Wettberg E, Nelson M, Ureta S. How Could the Use of Crop Wild Relatives in Breeding Increase the Adaptation of Crops to Marginal Environments? *Front Plant Sci*. 2022;13:886162.

97. Böhm J, Hoff B, O'Gorman CM, Wolfers S, Klix V, Binger D, et al. Sexual reproduction and mating-type-mediated strain development in the penicillin-producing fungus Penicillium chrysogenum. *Proc Natl Acad Sci U S A*. 2012;110(4):1476-81. https://doi.org/10.1073/pnas.1217943110 PMID: 23307807

98. O'Gorman CM, Fuller HT, Dyer PS. Discovery of a sexual cycle in the opportunistic fungal pathogen Aspergillus fumigatus. *Nature*. 2009;457:471-5. https://doi.org/10.1038/nature07528 PMID: 19043401

99. Galloway LD, Burgess R. Applied Mycology and Bacteriology 3rd Edition Leonard Hill, London. pp54 and 57. *Appl Mycol Bacteriol*. 1952;3:54-7.

100. Wick RR, Judd LM, Gorrie CL, Holt KE. Completing bacterial genome assemblies with multiplex MinION sequencing. *Microb Genom*. 2017;3(10):e000132. https://doi.org/10.1099/mgen.0.000132 PMID: 29177090

101. Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. *Genome Res*. 2017;27(5):722-36. https://doi.org/10.1101/gr.215087.116 PMID: 28298431

102. Walker BJ, Abeel T, Shea T, Priest M, Abouelliel A, Sakthikumar S, et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. *PLoS ONE*. 2014;9(11):e112963. https://doi.org/10.1371/journal.pone.0112963 PMID: 25409509

103. Langmead B, Wilks C, Antonescu V, Charles R, Hancock J, editor. Scaling read aligners to hundreds of threads on general-purpose processors. *Bioinformatics*. 2019; 35(3):421-32.

104. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics*. 2014; 30(15):2114-20.

105. Kurtz S, Phillippy A, Delcher AL, Smoot M, Shumway M, Antonescu C, et al. Versatile and open software for comparing large genomes. *Genome Biol*. 2004;5:R12.

106. Luo R, Liu B, Xie Y, Li Z, Huang W, Yuan J, et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. *GigaScience*. 2012;1(1):18. https://doi.org/10.1186/2047-217X-1-18 PMID: 23587118

107. Kolmogorov M, Armstrong J, Raney BJ, Streeter I, Dunn M, Yang F, et al. Chromosome assembly of large and complex genomes using multiple references. *Genome Res*. 2018;28(11):1720-32. https://doi.org/10.1101/gr.236273.118 PMID: 30341161

108. Kolmogorov M, Yuan J, Lin Y, Pevzner PA. Assembly of long, error-prone reads using repeat graphs. *Nat Biotechnol*. 2019;37(5):540-6. https://doi.org/10.1038/s41587-019-0072-8 PMID: 30936562

109. Solares EA, Chakraborty M, Miller DE, Kalsow S, Hall K, Perera AG, et al. Rapid Low-Cost Assembly of the Drosophila melanogaster Reference Genome Using Low-Coverage, Long-Read Sequencing. *G3 GenesGenomesGenetics*. 2018;8(10):3143-54.

110. Alonge M, Lebeigle L, Kirsche M, Jenike K, Ou S, Aganezov S, et al. Automated assembly scaffolding using RagTag elevates a new tomato system for high-throughput genome editing. *Genome Biol*. 2022;23(1):258.

111. Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. *Bioinformatics*. 2013; 29(8):1072-5.

112. Manni M, Berkeley MR, Seppey M, Simão FA, Zdobnov EM. BUSCO Update: Novel and Streamlined Workflows along with Broader and Deeper Phylogenetic Coverage for Scoring of Eukaryotic, Prokaryotic, and Viral Genomes. *Kelley J*, editor. *Mol Biol Evol*. 2021;38(10):4647-54.

113. Astashyn A, Tvedte ES, Sweeney D, Sapojnikov V, Bouk N, Joukov V, et al. Rapid and sensitive detection of genome contamination at scale with FCS-GX. *Genome Biol*. 2024;25(1):60.
114. Flutre T, Duprat E, Feuillet C, Quesneville H. Considering transposable element diversification in de novo annotation approaches. PLoS ONE. 2011;6(1):e16526. https://doi.org/10.1371/journal.pone.0016526 PMID: 21304975

115. Hoede C, Arnoux S, Moisset M, Chaumier T, Inizan O, Jamilloux V, et al. PASTEC: An Automatic Transposable Element Classification Tool. Cordaux R, editor. PLoS ONE. 2014; 9(5):e91929.

116. Wicker T, Sabot F, Hua-Van A, Bennetzen JL, Capy P, Chalhoub B, et al. A unified classification system for eukaryotic transposable elements. Nat Rev Genet. 2007;8(12):973-82. https://doi.org/10.1038/nrg2165 PMID: 17984973

117. Quesneville H, Bergman CM, Andrieu O, Autard D, Nouaud D, Ashburner M, et al. Combined evidence annotation of transposable elements in genome sequences. PLoS Comput Biol. 2005;1(2):e22. https://doi.org/10.1371/journal.pcbi.0010022 PMID: 16110336

118. Van Wyk S, Harrison CH, Wingfield BD, De Vos L, Van Der Merwe NA, Steenkamp ET. The RIPper, a web-based tool for genome-wide quantification of Repeat-Induced Point (RIP) mutations. PeerJ. 2019;7:e7447.

119. Smit A, Hubley R, Green P. RepeatMasker Open-4.0. [Internet]. 2013. Available from: http://www.repeatmasker.org

120. Smit A, Hubley R. RepeatModeler Open-1.0 [Internet]. 2008. Available from: http://www.repeatmasker.org

121. Palmer J, Stajich JE. Funannotate v1.8.1: Eukaryotic genome annotation [Internet]. 2020. Available from: https://zenodo.org/record/4054262

122. Brüna T, Hoff KJ, Lomsadze A, Stanke M, Borodovsky M. BRAKER2: automatic eukaryotic genome annotation with GeneMark-EP+ and AUGUSTUS supported by a protein database. NAR Genomics Bioinforma. 2021;3(1):lqaa108. https://doi.org/10.1093/nargab/lqaa108 PMID: 33575650

123. Hoff K, Lomsadze A, Borodovsky M, Stanke M. Whole-Genome Annotation with BRAKER. Methods Mol Biol. 2019;162:65-95. https://doi.org/10.1007/978-1-4939-9173-0_5 PMID: 31020555

124. Stanke M, Diekhans M, Baertsch R, Haussler D. Using native and syntenically mapped cDNA alignments to improve de novo gene finding. Bioinformatics. 2008; 24(5):637-44.

125. Hoff KJ, Lange S, Lomsadze A, Borodovsky M, Stanke M. BRAKER1: Unsupervised RNA-Seq-Based Genome Annotation with GeneMark-ET and AUGUSTUS. Bioinformatics. 2016;32(5):767-9.

126. Lomsadze A, Burns PD, Borodovsky M. Integration of mapped RNA-Seq reads into automatic training of eukaryotic gene finding algorithm. Nucleic Acids Res. 2014;42(15):e119. https://doi.org/10.1093/nar/gku557 PMID: 24990371

127. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403-10. https://doi.org/10.1016/S0022-2836(05)80360-2 PMID: 2231712

128. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10(1):421. https://doi.org/10.1186/1471-2105-10-421 PMID: 20003500

129. Brüna T, Lomsadze A, Borodovsky M. GeneMark-EP+: eukaryotic gene prediction with self-training in the space of genes and proteins. NAR Genomics Bioinforma. 2020;2(2):lqaa026.

130. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12(1):59-60. https://doi.org/10.1038/nmeth.3176 PMID: 25402007

131. Gotoh O, Morita M, Nelson DR. Assessment and refinement of eukaryotic gene structure prediction with gene-structure-aware multiple protein sequence alignment. BMC Bioinformatics. 2014 Dec;15(1):189. https://doi.org/10.1186/1471-2105-15-189 PMID: 24927652

132. Iwata H, Gotoh O. Benchmarking spliced alignment programs including Spaln2, an extended version of Spaln that incorporates additional species-specific features. Nucleic Acids Res. 2012;40(20):e161.

133. Lomsadze A, Ter-Hovhannisyan V, Chernoff YO, Borodovsky M. Gene identification in novel eukaryotic genomes by self-training algorithm. Nucleic Acids Res. 2005;33(20):6494-506. https://doi.org/10.1093/nar/gki937 PMID: 16314312

134. Punt M, van den Brule T, Teertstra WR, Dijksterhuis J, den Besten HMW, Ohm RA, et al. Impact of maturation and growth temperature on cell-size distribution, heat-resistance, compatible solute composition and transcription profiles of Penicillium roqueforti conidia. Food Res Int. 2020;136:109287. https://doi.org/10.1016/j.foodres.2020.109287 PMID: 32846509

135. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15-21.

136. Gabriel L, Hoff KJ, Brüna T, Borodovsky M, Stanke M. TSEBRA: transcript selector for BRAKER. BMC Bioinformatics. 2021;22(1):566. https://doi.org/10.1186/s12859-021-04482-0 PMID: 34823473

137. Krzywinski M, Schein J, Birol İ, Connors J, Gascoyne R, Horsman D, et al. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19(9):1639-45. https://doi.org/10.1101/gr.092759.109 PMID: 19541911

138. Jain C, Rodriguez-R LM, Phillippy AM, Konstantinidis KT, Aluru S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat Commun. 2018;9(1):5114.

139. Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods. 2012;9(7):671-5. https://doi.org/10.1038/nmeth.2089 PMID: 22930834

140. Frisvad JC, Filtenborg O. Classification of Terverticillate Penicillia Based on Profiles of Mycotoxins and Other Secondary Metabolites. Appl Environment Microbiol. 1983;46(6):1301-10.

141. Maechler M. diptest: Hartigan's Dip Test Statistic for Unimodality - Corrected. 2016. Available from: https://cran.r-project.org/package=diptest
142. Kieslich PJ, Henninger F, Wulff DU, Haslbeck JMB, Schulte-Mecklenbeck M. Mouse-tracking: A practical guide to implementation and analysis. In: A handbook of process tracing methods [Internet]. PsyArXiv; 2019 [cited 2023 Sep 5]. Available from: https://osf.io/zuvqa

143. Pfister R, Schwarz KA, Janczyk M, Dale R, Freeman JB. Good things peak in pairs: a note on the bimodality coefficient. Front Psychol. 2013;4:700. https://doi.org/10.3389/fpsyg.2013.00700 PMID: 24109465

144. Hall B, Hall M, Statificat L, Brown E, Hermanson R, Charpentier E, et al. Package "LapplacesDemon" [Internet]. 2021. Available from: https://github.com/LapplacesDemonR/LapplacesDemon

145. Adler D, Kelly T, Elliott T, Adamson J. vioplot: violin plot. [Internet]. 2022. Available from: https://github.com/TomKellyGenetics/vioplot

146. WickhamH. ggplot2: Elegant Graphics for Data Analysis. New York; 2016. Available from: https://ggplot2.tidyverse.org

147. Wickham H. Reshaping Data with the {reshape} Package. J Stat Softw. 2007;21(12):1–20.

148. Lê S, Josse J, Husson F. FactoMineR: An R Package for Multivariate Analysis. J Stat Softw [Internet]. 2008 [cited 2023 Aug 23]; 25(1). Available from: http://www.jstatsoft.org/v25/i01/

149. Kassambara A, Mundt F. factoextra: Extract and Visualize the Results of Multivariate Data Analyses [Internet]. 2020. Available from: https://CRAN.R-project.org/package=factoextra

150. Wei T, Simko V. R package "corrplot": Visualization of a Correlation Matrix [Internet]. 2021. Available from: https://github.com/taiyun/corrplot

151. Endelman JB. Ridge Regression and Other Kernels for Genomic Selection with R Package rrBLUP. Plant Genome. 2011;4(3):250–5. https://doi.org/10.3835/plantgenome2011.08.0024

152. Bloom JS, Ehrenreich IM, Loo W, Võ Lite TL, Kruglyak L. Finding the sources of missing heritability in a yeast cross. Nature. 2013;494(7436):234–7. https://doi.org/10.1038/nature11867 PMID: 23376951

153. Cheeseman K, Ropars J, Renault P, Dupont J, Gouzy J, Branca A, et al. Multiple recent horizontal transfers of a large genomic region in cheese making fungi. Nat Commun. 2014;5:2876. https://doi.org/10.1038/ncomms3876 PMID: 24407037

154. Katoh K. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30(14):3059–66. https://doi.org/10.1093/nar/gkf436 PMID: 12136088

155. Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, et al. Primer3--new capabilities and interfaces. Nucleic Acids Res. 2012;40(15):e115. https://doi.org/10.1093/nar/gks596 PMID: 22730293

156. Schuelke M. An economic method for the fluorescent labeling of PCR fragments. A poor man's approach to genotyping for research and high-throughput diagnostics. Nat Biotechnol. 2000;18:233–4. https://doi.org/10.1038/72708 PMID: 10657137

157. Missiaggia A, Grattapaglia D. Plant microsatellite genotyping with 4-color fluorescent detection using multiple-tailed primers. Genet Mol Res. 2006;5(1):72–8.

158. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9. https://doi.org/10.1038/nmeth.1923 PMID: 22388286

159. Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. GigaScience. 2021;10(2):giab008. https://doi.org/10.1093/gigascience/giab008 PMID: 33590861

160. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303. https://doi.org/10.1101/gr.107524.110 PMID: 20644199

161. Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011; 27(21):2987–93.

162. Taylor J, Butler D. R Package ASMap: Efficient Genetic Linkage Map Construction and Diagnosis. J Stat Softw [Internet]. 2017 [cited 2023 Aug 18]; 79:6. http://www.jstatsoft.org/v79/i06/

163. Broman KW, Wu H, Sen S, Churchill GA. R/qtl: QTL mapping in experimental crosses. Bioinformatics. 2003;19:889–90. https://doi.org/10.1093/bioinformatics/btg112 PMID: 12724300

164. Tang H, Zhang X, Miao C, Zhang J, Ming R, Schnable JC, et al. ALLMAPS: robust scaffold ordering based on multiple maps. Genome Biol. 2015;16(1):3. https://doi.org/10.1186/s13059-014-0573-1 PMID: 25583564

165. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2.

166. Li H. Minimap2: pairwise alignment for nucleotide sequences. Birol I, Bioinformatics. 2018;34(18):3094–100.

167. Gel B, Serra E. karyoplateR: an R/Bioconductor package to plot customizable genomes displaying arbitrary data. Hancock J, editor. Bioinformatics. 2017; 33(19):3088–90.