Analysis of Ciliate Community Diversity in Decaying Pinus nigra Logs

: Ciliates are an important component of the detritus and energy ﬂow in forest ecosystems. The present study aims to provide an early insight into the abundance and composition of the ciliate community inhabiting deadwood in relation to the different decay classes. We took advantage of a mesocosm experiment of black pine deadwood already underway to evaluate changes in chemical properties, microbial communities, and potential CO 2 emissions over time. The abundance and the number of ciliate taxa increased as wood decay progressed. Greater diversity was observed in the early stages of decomposition, while similarity in community composition increased along the decomposition gradient with several taxa commonly found in the more decomposed classes 3–5. The identiﬁed species were related to soil-inhabiting ciliates and mainly belonged to Colpodea and Spirotrichea classes. Ciliate abundance correlated positively with bacterial abundance, total nitrogen (N), and CO 2 potential production, while it correlated negatively with the C/N ratio. Through grazing activity, ciliates contribute to regulate the degrading activity of microbial communities inhabiting deadwood and CO 2 emission, enhancing soil fertility. Looking ahead, speciﬁc ciliate taxa may be used as indicators of the stage of decomposition and their biodiversity may provide knowledge into deadwood decay activity.


Introduction
Forest ecosystems harbor an immensurable percentage of terrestrial biodiversity, playing a significant role in global biogeochemical cycles, influencing climate, and providing ecosystem services for the well-being of society. Therefore, the maintenance of forest ecosystems currently constitutes an important instrument for biodiversity conservation and ecosystem functioning [1]. In the context of forest ecosystem maintenance, sustainable forest management (SFM) has become an important guiding principle for sustaining the biodiversity and functionality of forests [2], allowing for the diffusion of management practices in order to consider all components of forest ecosystems including deadwood [3]. In particular, in close-to-nature forest management, deadwood is considered a key component of forests, for both functional and structural point of view, because it supports habitats and species diversity [4,5].
In forest ecosystems, most of the specific diversity, even at intraspecific level, is found among bacteria, fungi, protists, lichens, and arthropods that live mainly in the soil and are essential for maintaining the efficiency of ecological processes [2]. In this sense, the investigations on deadwood abundance, typology, and decomposition mechanisms are of particular interest as deadwood has long been identified as one of the most suitable indicators for assessing the state of naturalness and therefore of forest biodiversity [3].

Experimental Design and Deadwood Sampling
A mesocosm experiment was set up consisting in 20 deadwood samples collected from lying deadwood randomly spread in the Pratomagno study area and distributed among five decay classes (4 samples within each decay class), incubated at controlled moisture and temperature, as previously described [20]. A system based on the main visual characteristics of deadwood, such as the occurrence of small branches and bark, the softness of wood, the rot development, and the extension of fungal mycelium (Table S1), was used to assign decay classes [22][23][24].
In addition, the DNA extracted from deadwood samples of decay classes 4 and 5 from the Monte Morello study area and used for a similar mesocosm experiment [21] were used to compare ciliate communities of two different locations.

DNA Extraction and Molecular Analysis
From each deadwood sample, the genomic DNA was extracted by using the FAST DNA SPIN kit for soil (Biomedicals, Santa Ana, CA, USA) as previously described [21].

PCR-DGGE
To estimate the composition and diversity of ciliate community in deadwood samples, we used the polymerase chain reaction-denaturing gradient gel electrophoresis (PCR-DGGE) method as described by Pastorelli et al. [20,21]. The ciliate 18S rRNA gene amplifications were carried out in a T100 Thermal Cycler (Bio-rad Laboratories, Hertfordshire, UK) by using a semi-nested PCR protocol. All PCR reactions were carried out in 25 µL containing 1X Xpert Taq reaction buffer (GRiSP Research Solutions, Porto, Portugal), 2 mM of MgCl 2 , 250 µM of deoxynucleotide triphosphates (dNTPs), 400 nM of each primer, 0.4 ng/µl of BSA, and 1U Xpert Taq DNA polymerase (GRiSP Research Solutions, Porto, Portugal). The first PCR step was performed with the primers CilF (5 -TGGTAGTGTATTGGACWACCA-3 ) and Cil959R (5 -TCTRAYCGTCTTTGATCCCYTA-3 ) [25]; 1 µL of the resulting product was diluted and re-amplified in a second PCR step using the same forward primer CilF carrying a 40 bp GC clamp [26] and the primer Cil940R (5 -TGAAAACATCCTTGGCAAATG-3 ) [18]. Cycling conditions consisted of an initial denaturation at 95 • C for 3 min followed by 35 cycles of denaturation at 94 • C for 30 s, 30 s of annealing (at 58 • C for the first step and at 52 • C for the second step), elongation at 72 • C for 45 s, and a final step of 5 min at 72 • C. The presence of amplified fragments of ca. 600 bp was verified by agarose gel electrophoresis (1.2 w/v), and the products of three independent PCR amplifications were pooled to minimize the effect of PCR biases.
The assessment of band migration distance and intensity within each lane of the DGGE was performed using Gel Compare II software v. 4.6 (Applied Maths, Saint-Martens-Latem, Belgium). The Shannon index was calculated considering the relative abundance of each band within a DGGE profile and was used as a proxy of the dominant taxon diversity.

Sequence Analysis
Selected dominant DGGE bands were excised from gels and sequenced, in order to detect representative taxa in the protistan deadwood community. The middle portion of each band was placed in 30 mL of distilled water and stored at −20 • C overnight. DNA fragments were eluted from gel through freezing and thawing, and 1 µL of each elution was used as a template in an amplification reaction using the primers pair CilF and Cil940R. The re-amplified PCR products were purified using PureLink TM Quick PCR Purification kit (Invitrogen-Life Technologies, Carlsbad, CA, USA) and cloned using the pGEM-T Easy Vector System (Promega) according to the manufacturer's guidelines. The Escherichia coli JM109 high-efficiency competent cells (Promega) were used for transformation of recombinant plasmids, and transformants were selected on LB/ampicillin plates. Five positive clones were chosen, screened by DGGE, and one for each type was subjected to DNA sequencing by BMR Genomics (Padova, Italy).
The absence of ambiguous peaks in sequence chromatograms was verified using Chromas Lite software (v2.1.1; Technelysium Pty Ltd.; Tewantin, Old, Australia) [27]. To obtain closely related nucleotide sequences, the converted FASTA files were compared against all sequences stored within the NCBI database using the Web-based BLAST tool [28]. For phylogenetic analysis, nucleotide sequences were aligned with other sequences of equivalent length from environmental clones and reference protists retrieved from the GenBank database, using the ClustalX 2.0.11 multiple sequence alignment software [29]. Distance analysis was performed according to the method of Jukes and Cantor [30], followed by phylogenetic tree construction performed by the neighbor-joining algorithm [31] using TREECON 1.3b [32]. The robustness of association nodes among samples was evaluated by bootstrap analysis with 1000 replicates.
Protistan taxonomic identification was achieved by means of different sequence similarity thresholds as described by Webster et al. [33].

Real-Time PCR
To determine the abundance of ciliate 18S rRNA gene in the deadwood samples, quantitative real-time PCR was performed as previously described [20]. Amplification reactions were carried out using the ABI StepOne TM real-time PCR system (Applied Biosystems, Foster City, CA, USA) in a 10 µl mixture containing 1X Blas Taq qPCR MasterMix (Applied Biological Materials Inc., Richmond, BC, Canada), 400 nM of each primer (CilF and Cil940R), distilled water (RNase/DNase free, Promega), and 2 µL of ten-or thousand-fold diluted DNA extracts or standard DNA. Cycling conditions consisted of incubation at 95 • C for 10 min, followed by 40 cycles of denaturation at 95 • C for 15 s, annealing at 55 • C for 60 s, and elongation at 72 • C for 45 s. All samples were run in triplicate in optical 96-well plates together with negative control and standard curve. Fluorescent light outputs were collected during each elongation step and analyzed with ABI Step One Plus Real-Time System SDS software (Version 2.3). Optimization assay conditions were performed for each template DNA concentration and a linear regression of the threshold cycle (C T ) for different DNA dilution vs. the log dilution of pooled DNA was used to estimated PCR efficiency (E = 10 −1/slope ) [34]. Melting curves were generated after amplification to verify the absence of primer dimers or artifacts. In addition, the amplification products were also checked on 1% agarose gel (w/v).
For absolute quantification, a standard curve was created using plasmid containing a target gene fragment obtained by the direct amplification of a sample of deadwood DNA with primer pairs CilF-Cil940R as described above. Four clones were randomly chosen and sequenced (BMR Genomics, Padova, Italy). The web-based BLAST tool available at the NCBI website [28] was used to check specificity.
Plasmid DNA was extracted using the QIAprep ® Spin Miniprep kit (Qiagen, Hilden, Germany), and plasmid concentrations were determined by spectrophotometry using a Nan-oDrop™ 1000 spectrophotometer (Thermo Scientific, Wilmington, NC, USA). Standard curves were freshly prepared with ten-fold dilutions ranging from 5.52 × 10 6 to 5.52 × 10 3 gene copies µL −1 . The abundance of ciliates in deadwood samples was expressed as 18S rRNA gene copy number in extracted DNA per mg of deadwood dry weight.

Physical and Chemical Properties
The physical and chemical properties of deadwood samples collected from Pratomagno study areas were previously determined by Pastorelli et al. [20] and are reported in Table 1. Data of the cumulative values of CO 2 potential production measured in deadwood samples are also reported in Table 1 [20]. Table 1. Moisture, total C (TC) and total N (TN) contents, C/N ratio, and pH of the black pine deadwood cores collected from the different decay classes in Pratomagno study area calculated as the percentage of the total dry mass of deadwood. Maximum values of CO 2 production by deadwood cores under closed laboratory systems are also reported [20]. Standard errors are in parentheses. Different letters in a column indicate significant differences at p < 0.05 (Tuckey test) among means. Briefly, moisture was determined by measuring dry weight of deadwood samples after incubation at 50 • C for 48 h. The pH was determined using a glass electrode meter (Hanna-pH211, Hanna Instruments, Padova, Italy) with a deadwood sawdust/water ratio of 1:10. The percentages of total carbon (TC) and the total nitrogen (TN) were measured by dry combustion using a Thermo Flash 2000 NC soil analyzer (Fisher Scientific, Waltham, MA, USA). The analysis was performed on 10 to 20 mg fine coarse sawdust deadwood weighted into Ag-foil capsules [20].

Statistical Analyses
To assess the significant differences in relation to decay class, data of deadwood physico-chemical properties, CO 2 production from deadwood, Shannon diversity index, and ciliate 18S rRNA gene copy number were analyzed by one-way analysis of variance (ANOVA) followed by Tukey post hoc test (p < 0.05) using Statistica 12 software (StatSoft, Palo Alto, CA, USA). The normality and the variance homogeneity of the data were tested prior to the ANOVA.
Pearson correlation analysis was performed among deadwood physico-chemical properties and maximum values of CO 2 emission from deadwood (Table 1), as well as abundance of total bacteria (data reported in Pastorelli et al. [20]) and ciliate (this study) by using the PAST4.03 software [35].
The banding patterns of each DGGE extracted as band intensity-matching tables were normalized by calculating the relative intensity of each band (ratio of the intensity of each band divided by the sum of the intensities of all bands in the same lane) and imported into PAST4.03 for multivariate statistical analysis. Non-metric multidimensional scaling (nMDS) was used to construct maps highlighting the similarity/dissimilarity of ciliate community structure among deadwood samples of each decay class. The accuracy of the nMDS plots was determined by calculating a stress value. The three-dimensional nMDS was chosen for Pratomagno-DGGE because the stress value (i.e., the distortion between actual similarity rankings and the corresponding distance ranking in the map), was less than 2.00 in contrast to the two-dimensional map, which showed a stress value close to 2.5. One-and two-way analysis of similarity (ANOSIM) followed by pairwise comparisons were conducted to determine the extent of differences in microbial communities among decay classes and/or location. nMDS and ANOSIM were performed using the Bray-Curtis distance measure and 9,999 permutational tests.
Correlation between ciliate communities and deadwood properties was determined by canonical correspondence analysis (CCA) performed with PAST4.03. DGGE bands were used as "species" data (filled symbols), while deadwood chemical properties (moisture, TN, TC, C/N, and pH) and CO 2 potential production were used as "environmental" variables (vectors); the statistical significance was assessed using a 999 permutational tests. The distances between the symbols reflects their dissimilarity; a symbol's position in relation to a vector head is a function of the correlation between the environmental parameter and the ciliate species. The length of vector reflects the relative importance of the environmental parameter in discriminating the microbial community of each decay class [36].

Nucleotide Sequence Accession Number
The nucleotide sequences determined in this study were deposited in the GenBank database under Accession numbers OK663130-OK663163. For sequences of different samples that resulted at 100% of nucleotide identity, only one was deposited in Gen-Bank database.

PCR-DGGE Analysis
In order to explore the effects of decay class and location on deadwood ciliate communities, two different DGGE analyses were performed: one including the black pine deadwood samples of the five decay classes collected from Pratomagno study area, and the other including deadwood samples of decay classes 4 and 5 collected from black pine deadwood in Pratomagno and Monte Morello. Differences in band position and intensity, depending on the decay class and/or location, were revealed within the DGGE profiles ( Figure S1).
A first visual inspection of the Pratomagno-DGGE profiles showed that several bands were common to deadwood samples from the most decomposed classes 3-5, with the number of common bands increasing as the decay progressed ( Figure S1). Significant differences in ciliate community diversity were found between decay class 1 (Shannon index = 2.16) and the more decomposed decay classes 2-5 (Shannon index ranging from 2.55 and 2.69). The nMDS showed an overlap of the symbols representing the ciliate communities of the five decay classes ( Figure 1); albeit, it is evident that the arrangement is not random and that ciliate communities tended to form clusters following in the progress of deadwood decomposition. In fact, symbols of early decay classes 1 and 2 were noticeably spaced apart and scattered in the nMDS maps, while those of decay classes 3-5 showed decreasing intersymbol distance and dispersion as well (Figure 1). To quantify community scatters, the intra-community distance was calculated as the Euclidean distance between symbols in the bi-dimensional cartesian plane of the nMDS graphic. Intra-community distance was higher in early decay classes (ranging from 0.35 ± 0.04 to 0.25 ± 0.02 of decay classes 1 and 2, respectively) and strongly decreased in the more decomposed deadwood classes (0.08 ± 0.01, 0.10 ± 0.04, and 0.07 ± 0.01 of decay classes 3, 4, and 5, respectively). differences in ciliate community diversity were found between decay class 1 (Shannon index = 2.16) and the more decomposed decay classes 2-5 (Shannon index ranging from 2.55 and 2.69). The nMDS showed an overlap of the symbols representing the ciliate communities of the five decay classes ( Figure 1); albeit, it is evident that the arrangement is not random and that ciliate communities tended to form clusters following in the progress of deadwood decomposition. In fact, symbols of early decay classes 1 and 2 were noticeably spaced apart and scattered in the nMDS maps, while those of decay classes 3-5 showed decreasing inter-symbol distance and dispersion as well (Figure 1). To quantify community scatters, the intra-community distance was calculated as the Euclidean distance between symbols in the bi-dimensional cartesian plane of the nMDS graphic. Intra-community distance was higher in early decay classes (ranging from 0.35 ± 0.04 to 0.25 ± 0.02 of decay classes 1 and 2, respectively) and strongly decreased in the more decomposed deadwood classes (0.08 ± 0.01, 0.10 ± 0.04, and 0.07 ± 0.01 of decay classes 3, 4, and 5, respectively). The distance between the same symbols is inversely proportionated to the degree of similarity between the ciliate community within a decay class. Three-dimensional (3D) nMDS was preferred since it led to lower distortion than twodimensional nMDS. The distance between the same symbols is inversely proportionated to the degree of similarity between the ciliate community within a decay class. Threedimensional (3D) nMDS was preferred since it led to lower distortion than two-dimensional nMDS.
The outcomes of the ANOSIM global test (Table 2) showed a significant effect of the decay class on ciliate community composition (p < 0.001), although the R value was close to 0.5; therefore, it was interpreted as separated but overlapping [37]. The pairwise comparison showed that the assemblages of ciliate communities of the early decay classes 1 and 2 significantly differed from those of the late and intermediate decay classes 3, 4, and 5 (Table 3). Regarding the nMDS performed on DGGE profiles ( Figure S2) obtained from 18S rRNA gene fragments amplified from deadwood samples of decay classes 4 and 5 collected from each study area, broad groupings by location and/or decay class can be highlighted ( Figure 2). The average intra-location distances between the symbols showed Pratomagno area symbols were more dispersed in the transformed plane than those of Monte Morello (0.36 ± 0.03 and 0.21 ± 0.02, respectively). DGGE profiles of Pratomagno were distinctly separate on the basis of decay class, with symbols of decay class 5 more closely grouped than those from decay class 4 (0.18 ± 0.04 and 0.25 ± 0.05, respectively). The opposite occurred for Monte Morello, which showed overlap between decay class 4 and class 5 with average distances of 0.19 ± 0.06 and 0.21 ± 0.10, respectively. . The distance between the same symbols is inversely proportionated to the degree of similarity between ciliate community within a decay class.

Sequence Analysis
By BLAST search, most of the bands excised from Pratomagno-DGGE were assigned to the Ciliophora phylum. Additional hits were obtained for some bands that were assigned to the Fungi and Nematoda phylum (Table S2). One band excised from a DGGE profile of decay class 1 resulted mostly closely related to P. nigra small subunit rRNA gene (Table S2), while two bands were associated to an uncultured soil eukaryote (accession number MK946118).
The phylogenetic tree reported in Figure 3 was obtained by align our ciliate sequences with other retrieved from GenBank. Overall, our sequences clustered within five ciliate classes. The predominance of the sequenced clones was related to sequences representative of the most common ciliate taxa found in soil, such as Colpodea and . The distance between the same symbols is inversely proportionated to the degree of similarity between ciliate community within a decay class.
The two-way ANOSIM global test (Table 2) showed that both location and decay class had significant effects on ciliate community composition (p < 0.01).

Sequence Analysis
By BLAST search, most of the bands excised from Pratomagno-DGGE were assigned to the Ciliophora phylum. Additional hits were obtained for some bands that were assigned to the Fungi and Nematoda phylum (Table S2). One band excised from a DGGE profile of decay class 1 resulted mostly closely related to P. nigra small subunit rRNA gene (Table S2), while two bands were associated to an uncultured soil eukaryote (accession number MK946118).
The phylogenetic tree reported in Figure 3 was obtained by align our ciliate sequences with other retrieved from GenBank. Overall, our sequences clustered within five ciliate classes. The predominance of the sequenced clones was related to sequences representative of the most common ciliate taxa found in soil, such as Colpodea and Spirotrichea. The remaining sequences were related to ciliate classes that have been described as being less present in soil, such as Oligohymenophorea, Nassophorea, and Litostomatea.    The nearest matches obtained by BLAST search and the overall taxonomical affiliation of the sequenced clones are shown in Table 4. The sequences were mainly taxonomically classified at species or genus level, showing a maximum identity with the nearest ciliate sequences of GeneBank ranging from 95.74% to 100%. Other clones were affiliated to the family level within the Oligohymenophorea and Spirotrichea classes. Table 4. BLAST analysis of 18S rDNA fragments selected from Pratomagno DGGE and belonging to the Ciliophora Phylum. The sequence coverages of the nearest matches were always 100%. The taxonomical identification was achieved by using different sequence similarity thresholds: a similarity of ≥97%, ≥95%, ≥90%, ≥85%, ≥80%, and ≥75% for assignment at the species-, genus-, family-, order-, class-, and phylum-level identification, respectively [33].

Clone
Nearest

Correlation Analysis
The CCA diagram ( Figure 5) showed the ciliate communities from the intermediate and late decay classes (classes 3, 4, and 5) clearly grouped separately from those of the early decay classes (classes 1 and 2). As decomposition progressed, ciliate communities increased their relationship with the environmental variables: moisture, TC, TN, and CO 2 flux from deadwood. A larger number of ciliate taxa were involved in separating communities of the early decay classes from the more decomposed ones. Six bands mostly contributed to separate the ciliate communities of decay classes 4 and 5 from those of the other classes, which were related to species belonged to classes Spirotrichea (namely Anteholosticha sigmoidea; Paraholosticha pannonica; Holosticha multistylata) and Colpodea (namely Platyophrya vorax; Colpoda aspera; Mykophagophrys terricola).

Figure 5.
Canonical correspondence analysis (CCA) ordination diagram of ciliate community and environmental variables defined by the first and second axes. DGGE band scores were also plotted as light grey dots; bands that significantly contribute to separate last decay classes were plotted as stars: Holosticha multistylata (Hm); Platyophrya vorax (Pv); Mykophagophrys terricola (Mt); Anteholosticha sigmoidea (As); Colpoda aspera (Ca); Paraholosticha pannonica (Pp). Vectors represent deadwood moisture, pH, deadwood total C content (TC), deadwood total N content (TN), deadwood C/N ratio, and CO2 flux from deadwood cores. The percentage of the variation in the data and significance is reported at each axis.

Discussion
We proved that the primer set developed to investigate ciliate diversity in soil is also suitable for studying this group of microorganisms in deadwood. Previous studies indicated a high level of specificity of this primer set and that other eukaryotic organisms, Figure 5. Canonical correspondence analysis (CCA) ordination diagram of ciliate community and environmental variables defined by the first and second axes. DGGE band scores were also plotted as light grey dots; bands that significantly contribute to separate last decay classes were plotted as stars: Holosticha multistylata (Hm); Platyophrya vorax (Pv); Mykophagophrys terricola (Mt); Anteholosticha sigmoidea (As); Colpoda aspera (Ca); Paraholosticha pannonica (Pp). Vectors represent deadwood moisture, pH, deadwood total C content (TC), deadwood total N content (TN), deadwood C/N ratio, and CO 2 flux from deadwood cores. The percentage of the variation in the data and significance is reported at each axis.

Discussion
We proved that the primer set developed to investigate ciliate diversity in soil is also suitable for studying this group of microorganisms in deadwood. Previous studies indicated a high level of specificity of this primer set and that other eukaryotic organisms, such as fungi, yeast, and algae, were not amplified in the second nested reaction [18,26]. Nevertheless, we also obtained DGGE bands referable to the fungal domain that were excised from the bottom of early decay classes 1 and 2 profiles, indicating a higher GC content of these bands than those referred to the ciliate phylum. Therefore, it is necessary to be careful in choosing the denaturation gradient of the gel, favoring gradients with a lower concentration of denaturing agents so as not to evidence DNA fragments richer in GC. One band excised from a DGGE profile of decay class 1 was attributed to P. nigra small subunit ribosomal RNA gene. In addition, some sequenced clones were referable to the nematode domain, suggesting that these retrieved bands might be faint and hidden in the background smear of DGGE profile and were amplified by cloning. Anyway, DGGE techniques demonstrated to be a valuable and inexpensive tool for rapidly monitoring the succession of the ciliate community in deadwood ecosystems and determining the taxonomic affiliation of the predominant species.
Observing the DGGE of the Pratomagno study area, the progress of wood decay reflected an increase in the number of ciliate taxa. In addition, an increasing abundance of total ciliates was also evidenced by real-time PCR, as decomposition progressed. These findings suggested that late decaying classes may contain a large number of niches and food resources to support and accommodate a large abundance of ciliates and species diversity. Moreover, we also observed an increase in the number of bands in common among profiles, indicating an increase in the similarity of the ciliate community structure along the decay gradient. This statement is supported by the nMDS analysis, which showed the symbols referring to the various DGGE profiles moving nearer in the three-dimensional space with the increasing of decay class. The same trend was observed in the bacterial and fungal communities of the same deadwood samples [20].
We hypothesized that the early stages of decomposition are distinguished by large differences in wood chemical-physical characteristics, as a probable consequence of stochastic events [38], which have led to colonization by bacterial and fungal taxa with different degrading abilities [21] and ciliate taxa with different prey preferences. Successively, as decomposition progresses, the deadwood environment becomes increasingly harsh and the microbial and ciliate communities assemble mainly according to deterministic mechanisms [39] controlled by the metabolic processes required for the degradation of complex wood residues and palatability of degrading microorganisms [6], respectively. According to our findings, previous research showed that beta-diversity can be increased in environments with high rates of stochastic extinction and immigration events, while it can be decreased where deterministic processes dominate [40,41].
In soil, ciliates are fairly affected by prevailing environmental factors such as overall habitat quality (climate, elevation), pH value, and food supply [2,10,25]. The ciliate community composition typically reflects fluctuations in moisture and temperature, as protists need water to be active, and temperature regulates moisture through drought and wetting [2,42]. Since we conducted a mesocosm experiment under controlled temperature conditions and constant internal humidity [20], our deadwood samples did not experience fluctuations in these two parameters during testing. Deadwood pH showed a not significant decreasing trend from the early decay classes 1 and 2 to the more decomposed ones (decay classes [3][4][5]. Given these claims, we might suppose that ciliate species in our deadwood samples were predominantly affected by nutrient availability changing along the decay gradient. Moreover, the abundance of ciliates was found to be positively correlated with CO 2 production and bacterial abundance. Several trophic functional groups can be highlighted within the deadwood saprotrophic community. Bacteria and fungi are the primary saprophytes. They breakdown the wood structural biopolymer, contributing to the release of nutrients [9]. Their number and activity increase with decay progression, and in the late decay classes, much of the deadwood biomass is immobilized into bacteria and fungi, especially as cell wall molecules (chitin and murein) [9]. The primary saprophytes are preyed upon by secondary saprophytes including specialized ciliates [10]. Therefore, the increase of bacteria and fungi is followed by the increase of ciliates' predators. Different ciliate species have specific feeding preferences, thus directly affecting the community dynamics of bacteria and fungi and, consequently, the expressions of associated decomposing enzymes [6].
By grazing, ciliates contribute to transform organic matter and provide nutrients to consumers inhabiting deadwood [6]. The organic molecules not digested by any species become a not recycled end-product and accumulate into soil, contributing to humus formation [9]. On the other hand, the excess N is excreted, contributing to enhance the deadwood TN content [20,21,43] and subsequently soil fertility [12].
In addition, we found ciliates to be negatively correlated to C/N consistently with the findings that bacterial propagation favors a low C/N ratio [44,45].
Anyway, the pattern of species diversity did not change much in the late decay classes. Several dominant bands were also in common between the deadwood samples of classes 4 and 5 collected in the Monte Morello areas, suggesting that ciliates are adapted to physical structure, chemical properties, and varieties in far-gone decomposed deadwood.
The ciliate species identified by sequence analysis were related to soil-inhabiting ciliates, and most of them belong to the Colpodea and Spirotrichea classes. A small number of excise bands were attributed to Litostomatea and Oligohymenophorea classes, while only one was attributed to Nassophorea. Colpodids are abundant in soil, especially in polluted soil [18,25], and are typically bacterivores [46]. Interestingly, Jia et al. [6] showed vesicular Colpodea strongly correlated to catalase and polyphenol oxidase, both enzymes involved in the degradation of refractory C sources, such as lignin, indicating a potential significant role of this group in deadwood degradation. Spirotrichs are common in soil, freshwater, and marine environments. Litostomatea mainly includes carnivores with some species preying primarily ciliates. Consistent with our results, Bartošová and Tirjaková [10] identified colpodids, spirotrichs, and litostomes as dominant systematic groups in bark and decaying wood. They found that the trophic guilds were very similar in deadwood and soils with the prevalence of bacterivores and predaceous in the more decayed samples. Coherently, we found a high abundance of both bacteria [20] and ciliates in the last decay class 5, suggesting a large presence of bacterivores. The increase of CO 2 production with the progress of wood decay indicates an increasing biomass or activity of organisms inhabiting deadwood and of the interactions among them.
Among sequenced bands identified as effective contributors to separating the ciliate communities of late decay classes, P. pannonica, H. multistylata, and M. terricola were common in decay classes 4 and 5 of both the Pratomagno and Monte Morello study areas (Figures S1 and S2) with the potential to be used as indicators of the state of deadwood decay. Indeed, a central challenge for upcoming research will be to establish a well-defined classification system of boundaries between the decay classes founded on objective and measurable variables, such as wood resistance and biological markers, that overcomes the five-class system based on the subjective visual assessment used so far.
We probably overestimated the ciliate function in deadwood, since we quantified the ciliate community based on extracted total DNA, thus not distinguishing between vegetative cells and cysts. On the other hand, it has been proposed that ciliate diversity in the natural environment is greater than previously proposed on the basis of the described free-living morphospecies [17]. Bartošová and Tirjaková [10] identified deadwood as a new frontier for discovering new species and genera of ciliates. Further high-resolution analyses based on next-generation sequencing (NGS) of both DNA and RNA should promote a detailed study at the species level of total ciliate community structure in deadwood by providing insights into community structures and helping to identify their roles in deadwood decomposition.
Supplementary Materials: The following are available online at: https://www.mdpi.com/article/ 10.3390/f13050642/s1. Figure S1: DGGE image of ciliate 18S rDNA gene fragments' profiles obtained from deadwood samples in the Pratomagno area. Samples were grouped on the basis of decay class, and the letter M indicates the marker used for the normalization of bands. The excised sequenced band related to Pinus nigra and those that significantly contributed to the separate last decay classes were pointed out: Holosticha multistylata, Hm; Platyophrya vorax, Pv; Mykophagophrys terricola, Mt; Anteholosticha sigmoidea, As; Colpoda aspera, Ca; Paraholosticha pannonica, Pp. Figure S2: DGGE image of ciliate 18S rDNA gene fragments' profiles obtained from deadwood decay classes 4 and 5 collected in the Monte Morello and Pratomagno areas. Samples were grouped on the basis of decay class, and the letter M indicates the marker used for the normalization of the band. Table S1: Description of visual characteristics used to assign the decay class. Table S2: BLAST analysis of 18S rDNA fragments selected from PCR-DGGE belonging to the other Eukariota domain. The sequence coverage of the nearest match was also reported. The taxonomical identification was achieved by using different sequence similarity thresholds: a similarity of ≥97%, ≥95%, ≥90%, ≥85%, ≥80%, and ≥75% for assignment at the species-, genus-, family-, order-, class-, and phylum-level identification, respectively.
Author Contributions: Conceptualization, R.P. and I.D.M.; methodology, R.P. and A.P.; formal analysis, R.P., M.A.C. and A.L.; data curation, R.P., M.A.C. and A.L.; writing-original draft preparation, R.P. All authors have read and agreed to the published version of the manuscript.
Funding: This work was supported by the EU's LIFE+ "Nature and biodiversity" program through the SelPiBio LIFE project (Innovative silvicultural treatments to enhance soil biodiversity in artificial black pine stands, i.e., LIFE13 BIO/IT/000282) for demonstration of innovative silvicultural treatments in artificial black pine stands.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.