species-to-species association networks in wood-inhabiting fungal communities

1833 –––––––––––––––––––––––––––––––––––––––– © 2020 The Authors. Oikos published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Wolf Eiserhardt Editor-in-Chief: Dries Bonte Accepted 10 August 2020 129: 1833–1843, 2020 doi: 10.1111/oik.07502 129


Introduction
Interspecific interactions are often hard to observe directly in species-rich communities. Consequently, inferring interspecific interactions indirectly from community data is of central interest in community ecology (Kissling et al. 2012, Morales-Castilla et al. 2015, Dormann et al. 2018) and fundamental for predicting community dynamics under the global environmental change (Araújo and Luoto 2007, Gilman et al. 2010, Wisz et al. 2013, Urban et al. 2016. However, data on species communities can be surveyed using different methods, each of which may differ in the amount and type of species that are detected, and thus produce diverging information on interaction networks (Wirta et al. 2014, Zhu et al. 2019. Data based on direct observations of species are prone to introducing false negatives due to imperfect detection (Guillera-Arroita 2017), whereas false positives due to species misidentifications are typical of data derived from DNAbased methods (Ficetola et al. 2015). Such observation error associated with each type of data is one of the reasons for the ongoing debate on which data types are valid for drawing ecological inferences (Yoccoz et al. 2001, Tyre et al. 2003, Lahoz-Monfort et al. 2014, Guillera-Arroita 2017. Given the exponential rise of DNA-based community datasets in ecology (Bush et al. 2017), assessing how the survey method affects the inferred interaction networks is a premise for advancing research on community ecology.
In species-rich communities where observation of biotic interactions is challenging, interactions are often estimated from snapshot data using network and co-occurrence analyses (Steele et al. 2011, Barberán et al. 2012, Weiss et al. 2016. Because species can co-occur more than expected by random not only due to interactions but also due to shared habitat requirements, the development of methods allowing to disentangle between the environmentally constrained and non-environmentally constrained co-occurrences is gaining increasing interest in the ecological literature (Peres-Neto et al. 2001, Ovaskainen et al. 2010a, Barberán et al. 2012, Legendre and Legendre 2012, Blois et al. 2013, Pollock et al. 2014, Morueta-Holme et al. 2016, Zurell et al. 2018. Recently developed joint species distribution models with latent variable structures enable inferring the residual associations (i.e. the co-occurrence patterns remaining after accounting for the effects of environmental predictors) from datasets consisting of large amounts of species (Warton et al. 2015, Ovaskainen et al. 2017). Since the residual species-tospecies associations can reflect both biotic interactions and shared responses to missing environmental variables, it is important to note that they should not be directly interpreted as interspecific interactions (Ovaskainen et al. 2010a, Dormann et al. 2018. However, compared to raw associations, residual association provide a closer approximation of species interactions.
Fungal species comprise megadiverse communities (Blackwell 2011) in which species interact in a number of ways, ranging from strong antagonisms to mutualism (Boddy et al. 2008). In particular, wood-inhabiting fungal communities are highly structured by interspecific interactions. Laboratory experiments have shown that the community development of wood-inhabiting fungi is largely affected by biotic interactions (Holmer and Stenlid 1997, Heilmann-Clausen and Boddy 2005, Fukami et al. 2010, Hiscox et al. 2015a. Likewise, field studies based on snapshot fruit-body observations suggest that during the decay process of dead wood, the identity of a successor fungal species is affected by the identity and abundance of the preceding fungal species (Renvall 1995, Pouska et al. 2013, Ottosson et al. 2014, Norberg et al. 2019). In addition, environmental conditions may influence the outcome of interspecific interactions  and conversely, fungal interactions may influence ecosystem functions such as wood decomposition (Boddy 2001, Fukami et al. 2010, Hiscox et al. 2015b, Maynard et al. 2018. To date, most studies exploring wood-inhabiting fungal communities have applied fruit-body surveys ). The method is relatively fast and low-cost, and thus enables large scale comparative studies with conspicuous species , Runnel et al. 2015. Nevertheless, several fungal groups produce ephemeral fruit bodies during limited fruiting seasons, making the number and timing of surveys critical for species' detectability (Halme and Kotiaho 2012, Abrego et al. 2016, Purhonen et al. 2017. Moreover, some wood-inhabiting fungi are microscopic (Lumley et al. 2000) and some individuals may not produce fruit bodies at all (Moore et al. 2008). Therefore, only a fraction of species detected by DNA-based metabarcoding methods are detected by fruit body-based (Allmér et al. 2006, Ovaskainen et al. 2010b, Kubartová et al. 2012, Ottosson et al. 2015, generating the debate on whether fruit body-based surveys are a valid method for studying fungal community dynamics compared to DNA-based surveys (Allmér et al. 2006, Porter et al. 2008, van der Linde et al. 2012, Runnel et al. 2015, Frøslev et al. 2019). Since fruit body-based surveys are able to capture only the part of the fungal community that is in the reproductive phase, these methods are clearly less optimal than DNA-based surveys for assessing the fungal diversity in terms of the number of species present. While Ovaskainen et al. (2013) showed that the ecological signal on the effects of dead wood properties on wood-inhabiting fungal communities is highly consistent between data based on fruit-body observations and DNA metabarcoding, it remains unclear whether the data collected by these two survey methods carry the same (or complementary) information on species-to-species associations for those species that can be detected with both methods.
In this paper, our aim is to assess whether the species-tospecies associations inferred from fruit body-based data reflect the ones inferred from DNA-based data. We hypothesize that DNA-based data allows estimating more associations compared to fruit body-based data, yet those associations detected from both data types are consistent. For comparing data collected by these survey methods, we use snapshot fruit-body and sequence data published in Ovaskainen et al. (2013) and Mäkipää et al. (2017). The data consist of 100 Norwegian spruce Picea abies logs of similar size that were simultaneously sampled by fruit-body and DNA-based surveys, the latter based on ITS metabarcoding of wood samples. We estimate and compare the species-to-species associations using a joint species distribution model. The model partitions the variation in species occurrences and co-occurrences to that explained by the abiotic predictors and variation remaining after accounting for these abiotic predictors, namely residual associations (Ovaskainen et al. 2016(Ovaskainen et al. , 2017.

Data collection
The data were collected from a protected, natural-like spruce dominated forest located in the municipality of Sipoo in Helsinki-Uusimaa Region, Finland (Ovaskainen et al. 2013). In this site, 100 large-sized logs (DBH 20-43 cm) were randomly selected. The selected logs varied in their decay stage (decay classes 1-4; see Hottola and Siitonen (2008) for class descriptions, and Supplementary Table 1 in Ovaskainen et al. (2013) for a table of decay stages), fall type (uprooted or broken) and ground contact but as minimally as possible in their size.
In November 2008, saw dust samples were taken from the 100 spruce logs for subsequent fungal DNA sequencing. Saw dust samples were taken separately from the basal and the middle part of the log with an electric drill. With strongly decayed logs, a sampling cylinder was used. The fruit-body surveys were carried out three times during the peak fruiting season. In September and November 2008, all logs were sampled twice for fruit bodies of all polypores and certain corticioid and hydnoid basidiomycetes that were easily identified in the field. During a third survey in October 2009, all polyporoid, corticioid and hydnoid basidiomycetes and some easily identified ascomycetes were sampled. If the species could not be reliably identified in the field, specimen samples were taken for later microscopic identification. During the latter two surveys, the abundance of species on a log was measured as the area covered by fruit bodies (cm 2 ) for all polypores and those corticioids and hydnoid species that were successfully identified in the field. For more detailed information on the sampling scheme, see Ovaskainen et al. (2013).

DNA extraction and sequencing
The protocols for DNA extraction and sequencing are described in detail in Ovaskainen et al. (2013). Briefly, DNA extractions for the 200 saw dust samples (two for each log) were done using Power Soil DNA isolation kit (MoBio Laboratories, Inc., Carlsbad, CA, USA). Taq polymerase enzyme (Thermo Fisher Scientific, Waltham, MA, USA) with the primers ITS1 and ITS4 (White et al. 1990) was used to amplify the target DNA region. Real-time quantitative PCR (qPCR) assays were conducted on the 18S ribosomal RNA gene. Samples with a DNA concentration ≥ 0.05 ng μm −1 were included in the qPCR analyses (163 out of 200 samples). CGX384 Real-Time PCR detection system was used to run the reactions. As a results, around 400 PCR fragments were obtained for the primers ITS1 and ITS2 to be used in the 454-sequencing.
For 297 samples with a DNA concentration > 1.43 ng μm −1 , PCR reactions for ITS1 region were done with a primer pair ITS1F (Gardes and Bruns 1993) and ITS2 (White et al. 1990). The primer pair ITS4 and ITS3 (White et al. 1990) was used for ITS2 region. For 97 samples with a DNA concentration < 1.43 ng μm −1 , a shorter form of ITS1F and ITS4 primers without the A and tag sequences was used for both ITS1 and ITS2 regions. Then, the PCR reactions were diluted and amplified again using ITS1F and ITS4 primers with tag and A sequences. Finally, sequencing was done using a Genome Sequencer FLX. In the case of six samples, no PCR fragments could be obtained.
We note that the DNA extraction and sequencing methods applied for Ovaskainen et al. (2013) represent the early stages of ITS metabarcoding techniques. Therefore, in these data the fungal community might be underrepresented, especially for those less abundant species requiring a deeper sequencing depth for being detected. However, Ovaskainen et al. (2013) showed a good correspondence between fruit-body and DNA abundance, suggesting a high signal-to-noise ratio in these data despite the limitations.

Molecular species identification
For this paper, we renewed the original molecular species identifications done for Ovaskainen et al. (2013) using PROTAX-fungi (Abarenkov et al. 2018). PROTAX-fungi is a software performing probabilistic taxonomic placement of fungal ITS sequences to known and unknown fungal lineages. The tool is based on a taxonomic classification system of Index Fungorum + Species Fungorum (Royal Botanic Gardens Kew et al. 2019) that covers a total of 131 484 species. The reference sequences are obtained from UNITE database (Kõljalg et al. 2013). PROTAX-fungi uses a database that consists of 420 319 reference sequences out of which 217 663 are annotated at the species level (Abarenkov et al. 2018). In our sequence data used for the identifications, samples for the ITS1 and ITS2 regions and for the basal and the middle part of the log were treated separately (ca 400 samples altogether, four per log). We included here those taxonomic placements that could be assigned to the species level with at least 50% probability.

Data preparation
We treated individual logs as sampling units, and thus we first pooled the data at the log level. For the data based on fruit-body observations, the area covered by fruit bodies of each species in the basal and the middle part of the logs were summed for each log, and then averaged over the surveys to yield a measure of fruit-body abundance for those species that were measured in the field (Ovaskainen et al. 2013). In the DNA-based data, we pooled the ITS1 and ITS2 sequences as well as the basal and the middle part of the logs by summing the species-specific sequence counts for all species detected. The sequence counts were used as a proxy of DNA abundance. Thus, the resulting data consisted of occurrence (i.e. presence-absence) and abundance (i.e. fruit-body coverage and number of sequences mapped to the focal species) of each species as fruit bodies and DNA on each log. The species nomenclature in fruit body-based data was updated according to Index Fungorum (Royal Botanic Gardens Kew et al. 2019) to match the nomenclature in the DNA-based data. Sequencing depth (i.e. the total number of sequences in a sample) was used as an estimate for sampling effort in the DNA-based data (Tedersoo et al. 2014).

Statistical methods
Since our aim was to test whether association networks estimated using each data type are consistent with each other, for the statistical analyses we selected only those fungal species that were detected both in fruit-body and DNA-based data. Namely, we did not estimate the associations among those species detected by only one of the survey methods, as those associations can obviously not be compared between datasets.
To obtain sufficient statistical power, only those species that occurred at least on 10 logs in at least one of the datasets (based on fruit-body observations or DNA) were included in the analyses. Such prevalence threshold was applied because if a given species occurs in a very small fraction of the sampling units, there is not sufficient statistical power to reliably estimate the species-to-species associations (Ovaskainen and Abrego 2020).
We conducted our analyses with the R package Hmsc (<www.r-project.org>, Tikhonov et al. 2020) which belongs to the class of joint species distribution models (Warton et al. 2015, Ovaskainen et al. 2017. Joint species distribution models allow inferring the species-to-species associations after controlling for the effects of the variables included in the model. The response matrix Y (see Ovaskainen et al. 2017 for notation) consisted of the occurrences and log-transformed abundances of each species as fruit bodies and DNA on the 100 sampling units (i.e. logs). Each species was thus included in the Y matrix in four columns representing the different types of data. We applied a hurdle approach that modelled the occurrences (presence-absences) with probit regression and log-transformed abundances conditional on presence with linear regression. Fruit-body and DNA abundances were standardized to zero mean and unit variance and considered conditional on presence, i.e. considered as missing data for cases in which the species were absent. Decay stage (categorical variable with the four levels), ground contact (present or absent), fall type (uprooted or broken) and log-transformed sequencing depth were included as the sampling unit level explanatory variables in the matrix X. To estimate the speciesto-species association networks, we included communitylevel random effects at the scale of sampling unit. The data type was included as a categorical variable with four levels (fruit-body occurrence, fruit-body abundance, DNA occurrence or DNA abundance) in the species trait matrix T. To disentangle between the raw species associations (i.e. those arising mainly from the shared log characteristics) and the residual associations (i.e. the ones remaining after controlling for log characteristics), we fitted two model variants. The first model variant with raw associations included sequencing depth as the only covariate (model 1), and the second one with residual associations included sequencing depth, decay stage, ground contact and fall type as covariates (model 2).
To sample the posterior distribution, we applied the default priors in the R-package Hmsc (Tikhonov et al. 2020) and performed 1 500 000 Markov chain Monte Carlo (MCMC) iterations with four chains. The first 500 000 iterations were removed as burn-in and the remaining ones were thinned by 1000, resulting in 1000 posterior samples per chain, and thus 4000 posterior samples in total. Before fitting the final model, we explored the rate of MCMC converge by fitting otherwise identical models but with 1500 (burn-in 500, thin 1), 15 000 iterations (burn-in 5000, thin 10) and 150 000 iterations (burn-in 50 000, thin 100). We followed Tikhonov et al. (2020) by evaluating MCMC convergence by computing the effective sample size and the potential scale reduction factor for the model parameters related to species-level (β) and community-level (γ) responses to environmental covariates, residual variation in species responses to environmental covariates (V) and the species association matrices (Ω).
We assessed the explanatory power by computing for each species and each data type the AUC for the occurrence data and R 2 for the abundance data. Then, we averaged for each data type the AUC values over all species and the R 2 values over those species for which abundance data was available for at least five logs. To examine the relative importance of environmental covariates and random effects (i.e. the associations) for the species communities, we performed variance partitioning among the explanatory variables.
To compare the species-to-species association networks identified from the fruit-body and DNA-based data, we classified the estimated associations as positive, negative or neutral, depending whether the posterior probability for the association was greater than 0.75, smaller than 0.25 or between 0.25 and 0.75, respectively. We note that the neutral case includes both those species pairs for which there is no association in reality, as well as those for which there may be an association in reality but the data do not have the sufficient statistical power for identifying it. We computed the numbers of species pairs for which the fruit-body and DNAbased data showed consistent results (both positive, both negative or both neutral), contradictory results (one positive and the other one negative), and for which the results were not consistent nor contradictory (one positive or negative and the other one neutral).

Results
Altogether, 142 species were observed in the fruit body-based surveys. The DNA sequencing resulted in a total of 1 349 965 sequences, out of which 18.5% (249 153) were assigned to a species level with at least 50% probability, yielding 413 species. Fruit-body and DNA-based data shared 80 species, out of which 37 species occurred at least 10 times in either fruitbody or DNA-based data and were consequently included in the HMSC-analyses. Out of these 37 target species, fruitbody abundance was available for 16 species. On average, each log hosted 6.2 target species as fruit bodies (ranging from 0 to 15) and 6.7 species as DNA (ranging from 1 to 15). The mean number of occurrences for the target species was 16.7 for fruiting species (ranging from 1 to 73) and 18.2 for species detected as DNA (ranging from 1 to 61). For a list of the 37 species included in the analyses and information on their prevalence and abundance, see Supplementary material Appendix 1 Table A1.
The MCMC convergence was satisfactory, as the effective sample sizes were close to the theoretical optimum of 4000, and the potential scale reduction factors were close to the theoretical optimum of one (Supplementary material Appendix 1 Fig. A1). Measured by the AUC statistic, the fitted model with residual associations explained on average 82% (respectively, 84%) of the variation in the occurrence of fruiting (respectively, present as DNA) species. Measured by R 2 , the fitted model with residual associations explained on average 27% (respectively, 29%) of the variation in species abundance (conditional on presence) of fruiting (respectively, present as DNA) species. Decay stage was the most important abiotic predictor, explaining 47% of the variance (on average over all species and data types). Ground contact, fall type and sequencing depth explained on average 9, 11 and 12% of the variance, respectively. The 22% of the variance belonged to the random effects at the log level. As expected, including the environmental covariates in the model improved the model fit (Supplementary material Appendix 1 Table A2) and decreased the proportion of variance explained by the random effect (Supplementary material Appendix 1 Table A3).
The total number of detected species-to-species associations was higher in the DNA than fruit body-based data (Fig. 1, 2). However, the directions of estimated speciesto-species associations were generally consistent between the data types, especially in the case of residual associations (Fig. 1b) that showed no contrary associations. In other words, species pairs showing positive (respectively, negative) associations in the fruit body-based data showed positive (respectively, negative) associations also in the DNA-based data (Fig. 1). In the occurrence data, a total of 666 residual species-to-species associations were detected. Out of these, 61.5% were consistent (both positive 11.4%, both negative 8.1% or both neutral 42.0%) between fruit-body and DNAbased data (Supplementary material Appendix 1 Table A4). No contrary associations were found (i.e. where one data type would support a positive association and the other one negative) (Supplementary material Appendix 1 Table A4). The remaining 38.4% represent cases in which the association was detected only in one data type (i.e. one data type giving support for a positive or negative association, and the other one for a neutral association) (Supplementary material Appendix 1 Table A4, Fig. 2). Such unique, positively and negatively associated species pairs were detected with both methods but the number of such pairs was larger for DNA-based data ( Fig. 2a-b). Out of all positive and negative residual associations in the occurrence data, 33.7% were detected by both methods (34.9% for raw associations), 18.7% only by fruitbody surveys (23.1% for raw associations) and 47.6% only by molecular surveys (41.9% for raw associations) ( Fig. 2a-b). Considering the type of detected association, both fruit-body and DNA-based data recorded more positive than negative associations in the occurrence data (Fig. 1, 2).
Compared to residual associations, raw associations in the occurrence data recorded a somewhat lower proportion (52.4%) of consistent associations (both positive 17.7%, both negative 7.2% or both neutral 27.5%) (Supplementary material Appendix 1 Table A4). They showed also some contradictory associations (1.2%), and more associations that were not consistent nor contradictory (46.5%) (Supplementary material Appendix 1 Table A4). In the case of abundance data (conditional on presence), only few associations were recorded in general (Fig. 1, 2c-d), making the comparison between fruit-body and DNA-based data less informative (Supplementary material Appendix 1 Table A4).

Discussion
In this study, we demonstrate that even if data derived from fruit body-based surveys are particularly prone to introducing false negatives compared to data derived from DNA-based surveys, these data types are valid for inferring species-tospecies association networks in wood-inhabiting fungal communities. Our results show that the directions of estimated residual species-to-species associations were consistent between the data derived from direct fruit-body observations and DNA-based surveys. However, there were differences in the amount of association links detected by each of the survey methods, and consistency of the links varied depending on whether the associations were estimated from models that accounted or did not account for the effects of environmental conditions. Below, we will explain each finding in turn. Figure 1. The estimated species-to-species associations measured as (a) raw and (b) residual correlations, separately for each data type (species occurrences (occur.) and abundances (abun.) in data based on fruit-body observations and DNA metabarcoding). The matrices show species pairs with positive (red) and negative (blue) associations with at least 75% posterior probability. The color shade indicates the strength of the correlation (the darker the shade, the stronger the correlation). The remaining cases, i.e. neutral associations, are shown in white. Each data type contains the same set of 37 species (except for fruit-body abundance data including 16 out of these species) in a corresponding order (for species names, see Supplementary material Appendix 1 Table A1). The order is set to highlight positively and negatively associated groups.
Compared to more classical species surveys, DNA metabarcoding methods yield more species-rich community datasets and represent more inconspicuous species (Creer et al. 2016). This is also the case for wood-inhabiting fungi (Allmér et al. 2006, Kubartová et al. 2012, Ovaskainen et al. 2013, Ottosson et al. 2015. Supporting these previous findings and our hypotheses, our results showed that the DNAbased survey detected more species than the survey based on fruit-body observations. However, with DNA-based methods not only more species are detected per log, but also each species is likely to be detected in a larger number of logs than with fruit-body surveys, providing stronger statistical power for inferring species-to-species associations. Consequently, we estimated a larger absolute number of associations from data based on DNA than from data based on fruit-body observations. The result is linked to the observational error in fruit body-based surveys. Only a fraction of the species found as DNA are detected as fruit bodies since only the most abundant species as DNA produce fruit bodies and they appear with a time delay from the colonization (Allmér et al. 2006, Kubartová et al. 2012, Ovaskainen et al. 2013, Ottosson et al. 2015. Furthermore, most of the species produce ephemeral fruit bodies emerging at different times during the fruiting season (Halme and Kotiaho 2012, Abrego et al. 2016, Purhonen et al. 2017, fruit bodies of some fungal species are too small to be detected with naked eye Figure 2. Euler diagrams showing the number of positively (red), negatively (blue) and neutrally (grey) associated species pairs estimated from the data based on fruit-body observations (FB) and DNA metabarcoding (DNA) (numbers on the left-and right-hand sides of the diagrams, respectively), and the number of species pairs recorded by both survey methods (numbers within the overlapping areas). The first and second rows of panels show (a) raw and (b) residual associations in the occurrence data, and the third and fourth rows of panels show (c) raw and (d) residual associations in the abundance data. The sizes of the circles correspond to the total number of associations detected by each survey method. If a data type detected no associations, this is shown with a cross. (Cannon andSutton 2004, Abrego andSalcedo 2015) and some fungal individuals may never produce fruit bodies (Moore et al. 2008). Moreover, supporting our results, interspecific interactions in fungal communities take place mainly at the mycelial stage (Heilmann-Clausen and Boddy 2005, Fricker et al. 2008, Hiscox et al. 2018, which can be surveyed using the DNA-based methods only. Before the potential fruiting, species have to occupy space and resources within dead wood as mycelium and while doing so, they necessarily interact with other species (Boddy 2000, Boddy andHiscox 2016).
There were also associations that were estimated only from the data based on fruit-body observations. Most likely these species pairs represent cases in which sequences have been misidentified in the DNA-based data. A large body of literature has demonstrated that when applying molecular species identification methods to data originating from high-throughput sequencing of environmental samples, species can easily be misidentified (Meyer and Paulay 2005, Lou and Golding 2012, Ficetola et al. 2015, Lahoz-Monfort et al. 2016). There are myriad reasons for this, such as sequencing errors, mislabeling errors in the reference databases, and the fact that some species are simply missing from the reference databases (Lou andGolding 2012, Somervuo et al. 2017). As a consequence, some species that actually occurred in sampling units might be missing from the data, while some absent species might be falsely recorded to be present. To take these issues into account as much as possible, we applied a probabilistic taxonomic placement method which accounts for these sources of uncertainty while estimating the reliability of taxonomic identifications (Somervuo et al. 2016, Abarenkov et al. 2018. Also, by including sequencing depth as a covariate in the models, we accounted for false absences due to limited observation effort (Ovaskainen et al. 2017). We also note that some species might remain undetected by the DNA-based survey also because saw dust samples are limited in volume (Allmér et al. 2006, Runnel et al. 2015, affecting the number of species detected (Ovaskainen et al. 2010b, Kubartová et al. 2012. Furthermore, the data applied in this study originates from the early phases of ITS metabarcoding, and thus it is not fully state-of-the-art in terms of the DNAextraction and sequencing methods. More modern methods enable, for example, a higher sequencing depth which decreases the likelihood for missing species that are present in low DNA abundances. However, as here we focused only on the dominating part of the species community, we expect our results to be robust even if they can be somewhat conservative due to the noise generated by relatively modest sequencing depth. The level of consistency between the species-to-species associations inferred from different data types depended on how the associations were estimated. By adding the environmental covariates in the model, the proportion of consistent associations increased compared to the model without the covariates. This supports the use of residual species-tospecies associations as a hypothesis of species interactions, as has been done in previous studies (Ovaskainen et al. 2010b, Ottosson et al. 2014, Abrego et al. 2017. The result also demonstrates that some species co-occur merely due to shared responses to certain abiotic conditions. In specific, the production of fungal fruit bodies depends strongly on the environmental conditions (Moore et al. 2008). Therefore, it is highly relevant to account for the effects of these conditions when inferring species-to-species associations from data based on fruit-body observations.
We emphasize that while some of the residual associations detected in our study might represent biotic interactions, these can only be confirmed through experimental tests (Dormann et al. 2018). The detected positive residual species-to-species associations might reflect facilitative biotic interactions (Tiunov andScheu 2005, Fukami et al. 2010) or parasitism (Niemelä et al. 1995), whereas the negative associations might reflect competitive interactions among wood-inhabiting fungi (Boddy 2000, Heilmann-Clausen andBoddy 2005). Species-to-species associations might also reflect community succession and priority effects (Ottosson et al. 2014, Norberg et al. 2019. For instance, the predecessor species Fomitopsis rosea showed a negative association with the successor Phellopilus nigrolimitatus both in the fruit-body and DNA occurrence data, indicating that the later species replaces the former during species succession (Niemelä et al. 1995). Given the high predictive power of our model including environmental covariates, it is plausible that the detected residual species-to-species associations represent interspecific interactions. However, some of the associations most likely reflect a shared response to some missing environmental covariate that is relevant for the wood-inhabiting fungal species included in our study. For instance, we included decay stage of the logs as a proxy of the physiochemical properties of the wood. Yet, a direct measurement of physiochemical properties could have explained the occurrence of some species better. We also note that while competition is known to be the most common interaction type among wood-inhabiting fungi (Boddy 2000), we found more positive than negative species-to-species associations. This is partially because when inferring species-to-species associations using statistical methods, more data are needed for inferring negative than positive associations (Ovaskainen et al. 2016). In addition, also competing species might co-occur on a sampling unit level if they are spatially separated within a log, leading to a positive association at the log level (Ovaskainen et al. 2010a, Ottosson et al. 2014. Finally, while a large number of residual species-to-species associations were found from occurrence data, only few associations were found from the abundance data. Namely, the fact that a given species was abundant did not affect the abundances of other species. This suggests that wood-inhabiting fungal interactions take place when the potentially interacting species are present, the interactions not depending on the abundance of interacting species. On the other hand, there might not have been enough statistical power to estimate associations from the abundance data as most of the abundance data from fruit-body surveys were missing. Also, in the case of sequence data, read counts are affected by a range of factors and cannot be directly translated to abundance data (Deagle et al. 2019).
We conclude that even though fruit-body surveys focus only on the reproductive part of the wood-inhabiting fungal community, data derived from fruit-body and DNA-based surveys yield consistent association networks of wood-inhabiting fungi. This supports the robustness of our results and therefore, both methods should lead to comparable ecological inference on the association network structure of woodinhabiting fungi. DNA-based surveys were more informative on the total number of interacting species but since unique interactions were also estimated from the fruit-body data, a combination of the methods would be ideal to obtain the most comprehensive insight into the network structure of wood-inhabiting fungi. Finally, our study demonstrates that accounting for the effects of relevant environmental covariates is of particular importance when inferring species-tospecies associations in wood-inhabiting fungal communities.

Data availability statement
The original data used in this study are published in Ovaskainen et al. (2013) and Mäkipää et al. (2017). The data formatted for this study as well as the R scripts for reproducing the results are available from the Dryad Digital Repository: < http://dx.doi.org/10.5061/dr yad.tmpg4f4wm> (Saine et al. 2020).