Dynamic Shifts in the Root Microbiota of Cultivated Paphiopedilum armeniacum during Different Stages of Growth

: Paphiopedilum armeniacum S. C. Chen et F. Y. Liu is an endangered lady’s slipper orchid species with high horticultural value. As observed for other orchids, mycorrhizal fungi and endophytic bacteria play important roles in the growth and development of P. armeniacum . In the present study, the community structure dynamics across three growth and development stages of cultivated P. armeniacum were investigated. The potential interactions between Tulasnellaceae fungi and core bacterial genera on one hand and the stability of the presumed mycorrhizal fungi communities on the other were analyzed in three growth stages of P. armeniacum to enhance our understanding of endophytic microbial community structure dynamics in the roots at different development stages. Based on sequencing, 3 and 16 phyla and 59 and 269 genera were identiﬁed in the fungal and bacterial communities, respectively. The predominant fungi and bacteria were Basidiomycota (62.90%) and Proteobacteria (43.98%), which exhibited changes in abundance and diversity depending on the growth stage of P. armeniacum . Assessment of the entire microbial communities from different growth stages showed that the seedling stage had the highest richness and diversity. The microbial communities recruited by P. armeniacum at the seedling stage were different from those recruited at the vegetative and reproductive growth stages, and the microbial communities recruited in the latter two stages overlapped. Tulasnellaceae were the only dominant fungal symbionts during P. armeniacum growth. Brevibacillus , Mycobacterium , and Sphingomonas , the three core genera, showed signiﬁcant interactions with the main OTUs of Tulasnellaceae. Putative mycorrhizal fungi in P. armeniacum were relatively stable across different growth environments, and the core mycorrhizal fungi were uncultured Tulasnellaceae (OTU1). This could facilitate the ex situ conservation and commercial development of the endangered orchid.


Introduction
Paphiopedilum Pfitzer is the largest genus of the subfamily Cypripedioideae, which comprises about 96-100 species around the globe [1,2]. Paphiopedilum armeniacum S. C. Chen et F. Y. Liu is an endangered orchid distributed in Baoshan, Gongshan, Fugong, Yunlong, Lanping, and Lushui counties in western Yunnan, China, as well as in neighboring Myanmar [3]. It is an attractive and popular horticultural species due to its unique labellum and its bright yellow flowers. Despite their high horticultural value, wild P. armeniacum populations are under a constant threat of extinction due to overcollection and habitat destruction. As a result, P. armeniacum is listed as endangered in the IUCN Red List of Threatened Species, as well as in the Convention on International Trade in Endangered Species of Wild Fauna and Flora Appendix I, and its international trade is prohibited [4].

Sample Collection and Surface Sterilization
To assess endophytic microbial community structures across different development stages of cultivated P. armeniacum, 15-20 plants in each stage were collected randomly. We collected 1-2 healthy roots from each plant. All plants were obtained from seed germination by tissue culture in vitro. The seedling-stage plants were tissue culture plants transplanted into the greenhouse for 1 month, with a height of about 7-9 cm, and are referred to as phase S. The vegetative growth-stage plants had been planted in the greenhouse for two years and are referred to as phase M. Finally, the reproductive growth-stage plants had experienced their first bloom and are referred to as phase L ( Figure 1). domly selected. After microscopy, they were surface-sterilized immediately (30 s submergence in 75% ethanol, followed by a 30 s rinse step in sterile distilled water, 5 min submergence in 1% sodium hypochlorite, and five 30 s rinse steps in sterile distilled water). The sterilized roots from each of the growth stages were then cut separately into 2-3 mm small pieces, evenly mixed in equal proportions, and divided into three replicates for molecular analysis. To confirm the sterilization process was successful, the final water rinse was also collected for molecular analysis as a negative control.

DNA Extraction and PCR Amplification
DNA was extracted from samples and negative control using the QIAamp DNA Stool Mini Kit (QIAGEN, Hilden, Germany), according to the manufacturer's instructions. DNA quality was checked by 1% agarose gel electrophoresis, and DNA concentrations and purity were determined using a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, NC, USA). Primers 799F (5′-AACMGGATTAGATACCCKG -3′) and 1193R (5′-ACGTCATCCCCACCTTCC-3′) targeting the V5-V7 regions of 16S rRNA genes and Primers ITS3F (5′-GCATCGATGAAGAACGCAGC-3′) and ITS4-OF (5′-GTTACTAGGGGAATCCTTGTT-3′) targeting the fungal ITS-2 regions were used for PCR [30,31]. The library was constructed by two-step PCR amplification. The first step involved amplifying the target fragment with specific primers. The initial PCR reactions were performed in 50 μL reaction volumes with 1-2 μL of DNA template, 1 μL of dNTPs at 10 mM, 5× reaction buffer, and 1U of Phusion DNA Polymerase (New England Biolabs, Ipswich, MA, USA). PCR conditions consisted of initial denaturation at 94 °C for 2 min, All the samples were collected on 3 April 2019 from a glass greenhouse with a water curtain and blowers for cooling and ventilation in Xingyi, Guizhou, southwest China. The plants were potted in a bark substrate for orchids, under no more than 800 µmol·m −2 ·s −1 of natural light maintained by a sunshade net. The average temperature and relative humidity were 10-32 • C and 70-90%, respectively. In the process of cultivation, stability and proportion of the substrate were maintained as far as possible. Meanwhile, part of the original substrate was retained during the repotting process.
After careful excavation from the culture substrate, the roots were transported to the laboratory at 4 • C. All the vegetative-and reproductive-stage roots were checked microscopically for mycorrhizal colonization (Figure 1), and the root segments with mycorrhizal colonization were selected for molecular analyses. The distinct intracellular fungal pelotons were not obvious in the seedling stage, hence fifteen seedling-stage roots were randomly selected. After microscopy, they were surface-sterilized immediately (30 s submergence in 75% ethanol, followed by a 30 s rinse step in sterile distilled water, 5 min submergence in 1% sodium hypochlorite, and five 30 s rinse steps in sterile distilled water). The sterilized roots from each of the growth stages were then cut separately into 2-3 mm small pieces, evenly mixed in equal proportions, and divided into three replicates for molecular analysis. To confirm the sterilization process was successful, the final water rinse was also collected for molecular analysis as a negative control.

DNA Extraction and PCR Amplification
DNA was extracted from samples and negative control using the QIAamp DNA Stool Mini Kit (QIAGEN, Hilden, Germany), according to the manufacturer's instructions. DNA quality was checked by 1% agarose gel electrophoresis, and DNA concentrations and purity were determined using a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, NC, USA). Primers 799F (5 -AACMGGATTAGATACC-CKG -3 ) and 1193R (5 -ACGTCATCCCCACCTTCC-3 ) targeting the V5-V7 regions of 16S rRNA genes and Primers ITS3F (5 -GCATCGATGAAGAACGCAGC-3 ) and ITS4-OF (5 -GTTACTAGGGGAATCCTTGTT-3 ) targeting the fungal ITS-2 regions were used for PCR [30,31]. The library was constructed by two-step PCR amplification. The first step involved amplifying the target fragment with specific primers. The initial PCR reactions were performed in 50 µL reaction volumes with 1-2 µL of DNA template, 1 µL of dNTPs at 10 mM, 5× reaction buffer, and 1U of Phusion DNA Polymerase (New England Biolabs, Ipswich, MA, USA). PCR conditions consisted of initial denaturation at 94 • C for 2 min, followed by 29 (16S rRNA) or 35 (ITS) cycles of denaturation at 94 • C for 30 s, annealing at 55 • C for 30 s, and extension at 72 • C for 30 s, with a final extension at 72 • C for 5 min. The PCR products were recovered using an AxyPrep DNA gel recovery kit (AXYGEN Inc., Union city, CA, USA) and quantified using an FTC-3000TM Real-Time PCR instrument; the samples were mixed according to a mole ratio similar to that for secondary PCR amplification. The role of the secondary PCR amplification was to add the connector and the sequence primer and barcode required by the Illumina platform to both ends of the target fragment. The PCR reactions were carried out in 40 µL reaction volumes with 5 µL of DNA template, 1 µL of dNTPs at 10 mM, 8 µL of 5× reaction buffer, and 1U of Phusion DNA polymerase. The cycling conditions consisted of one cycle at 94 • C for 2 min, followed by eight cycles at 94 • C for 30 s, 56 • C for 30 s, and 72 • C for 30 s, followed by a final extension cycle at 72 • C for 5 min. The subsequent sequencing analysis was not performed on the negative control, as the target band was not detected in the first PCR amplification.

Covariation Network Analysis of Core Bacterial Genera and Mycorrhizal Fungi
To explore the interactions between root core bacterial communities and mycorrhizal fungi of P. armeniacum, a covariation network diagram was constructed based on the correlation coefficients between core bacterial genera and mycorrhizal fungi.
The OTUs that appeared in all samples were defined as "core OTUs" [36,37]. Only core bacterial genera accounting for more than 1% of the total reads were used in the calculations. The Spearman rank correlation coefficient between two different populations was calculated using R v3.6.0, and the co-variation networks of different populations were constructed using igraph v1.2.4.2 [38]. The correlation P value less than 0.05 was plotted.

Comparative Analysis of Tulasnella in Paphiopedilum armeniacum
To analyze the stability of the Tulasnellaceae fungal communities in P. armeniacum samples, we generated a presumed Tulasnellaceae dataset of P. armeniacum using data gathered from a previous study [14], in addition to data obtained from the present study. Phylogenetic trees were generated using ITS sequences of each of the sequences obtained from BLAST searches (Fasta S1). All sequences were aligned in MAFFT v7.311 [39] and manually adjusted in MEGA X [40]. MEGA X was used to construct a maximum likelihood  16 November 2021)). Branches that received bootstrap support for ML greater than 50% or equal 0.70 were considered significantly supported.

Results
A total of 423,337 and 379,593 high-quality reads were obtained from the ITS2 and 16S rRNA gene sequencing data, respectively. After discarding non-bacterial, non-fungi, mitochondrial, chloroplast, and low-abundance OTUs, we obtained 199 and 828 fungal and bacterial OTUs, respectively, at 97% similarity. After subsampling, 813 bacterial OTUs and 195 fungi OTUs were observed (Table S1). The observed OTUs and the numbers of unique or shared OTUs identified in different development stages are illustrated in Figure 2. The seedling stage had the highest number of observed OTUs ( Phylogenetic trees were generated using ITS sequences of each of the sequences obtained from BLAST searches (Fasta S1). All sequences were aligned in MAFFT v7.311 [39] and manually adjusted in MEGA X [40]. MEGA X was used to construct a maximum likelihood (ML) tree under the K2+G+I model, with 1000 bootstrap replicates. Phylogenetic trees were visualized using FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/ (accessed on 16 November 2021)). Branches that received bootstrap support for ML greater than 50% or equal 0.70 were considered significantly supported.

Results
A total of 423,337 and 379,593 high-quality reads were obtained from the ITS2 and 16S rRNA gene sequencing data, respectively. After discarding non-bacterial, non-fungi, mitochondrial, chloroplast, and low-abundance OTUs, we obtained 199 and 828 fungal and bacterial OTUs, respectively, at 97% similarity. After subsampling, 813 bacterial OTUs and 195 fungi OTUs were observed (Table S1). The observed OTUs and the numbers of unique or shared OTUs identified in different development stages are illustrated in Figure  2. The seedling stage had the highest number of observed OTUs (

Diversity Analysis in Different Samples
Microbiota diversity and richness in different samples were investigated using Chao1 and Shannon indexes based on ITS and V3-V4 sequence data. According to the Chao1 index results, the seedling stage samples had higher species richness than the

Diversity Analysis in Different Samples
Microbiota diversity and richness in different samples were investigated using Chao1 and Shannon indexes based on ITS and V3-V4 sequence data. According to the Chao1 index results, the seedling stage samples had higher species richness than the samples from the other stages in both the fungal and the bacterial dataset ( Figure 3A,C). In addition, significant differences were observed in fungal richness among the samples ( Figure 3C, p < 0.05). Similar results were observed based on the Shannon index values, with higher diversity in fungal and bacterial communities in seedling-stage samples than in samples from the other two stages, with significant differences in bacterial diversity among the samples ( Figure 3B, p < 0.05).
In addition, significant differences were observed in fungal richness among the samples ( Figure 3C, p < 0.05). Similar results were observed based on the Shannon index values, with higher diversity in fungal and bacterial communities in seedling-stage samples than in samples from the other two stages, with significant differences in bacterial diversity among the samples ( Figure 3B, p < 0.05).
In the bacterial microbiomes, measures of within-sample diversity (α-diversity indices) showed decreases in microbial richness and diversity with the development of P. armeniacum based on the Chao1 index and Shannon index, respectively, showing that in both bacterial and fungal communities, higher richness and diversity were observed in the seedling stage than in the other development stages. The relationships among the microbial communities in the different samples were investigated using Principal Coordinate Analysis (PCoA) based on Bray-Curtis distance matrix and unweighted UniFrac distances. The PCoA results indicated that the microbiota in the seedling stage samples was distinct from that in the other stages.
The two axes explain 76.2% of the total variation in V3-V4 sequences ( Figure 4A), as well as 81.9% of the total variation in ITS sequences ( Figure 4C). In the heatmap, different distances were observed between the microbial communities of individual samples (Figure 4). In the case of the fungal community, the seedling-stage samples tended to cluster closer, while samples in the other stages formed another cluster. It was not obvious in hierarchical clustering analysis of the bacterial community, but the PCoA also highlighted differences between the seedling stage and the other two stages. In the bacterial microbiomes, measures of within-sample diversity (α-diversity indices) showed decreases in microbial richness and diversity with the development of P. armeniacum based on the Chao1 index and Shannon index, respectively, showing that in both bacterial and fungal communities, higher richness and diversity were observed in the seedling stage than in the other development stages.
The relationships among the microbial communities in the different samples were investigated using Principal Coordinate Analysis (PCoA) based on Bray-Curtis distance matrix and unweighted UniFrac distances. The PCoA results indicated that the microbiota in the seedling stage samples was distinct from that in the other stages.
The two axes explain 76.2% of the total variation in V3-V4 sequences ( Figure 4A), as well as 81.9% of the total variation in ITS sequences ( Figure 4C). In the heatmap, different distances were observed between the microbial communities of individual samples (Figure 4). In the case of the fungal community, the seedling-stage samples tended to cluster closer, while samples in the other stages formed another cluster. It was not obvious in hierarchical clustering analysis of the bacterial community, but the PCoA also highlighted differences between the seedling stage and the other two stages.

Taxonomic Assignment of the Microbial Community Composition
In total, 3 and 16 phyla and 9 and 27 classes were identified in the fungal and bacterial communities, respectively. In the bacterial communities ( Figure 5A, Table S1), the dominant phyla (>1%) were Proteobacteria, Firmicutes, and Actinobacteria, with average relative abundances of 43.98%, 38.51%, and 15.81%, respectively. In the fungal communities ( Figure 5B, Table S1), the dominant phyla across all samples were Ascomycota and Basidiomycota, with average relative abundances of 62.94% and 18.60%, respectively. The average relative abundance of Proteobacteria in the seedling stage was higher than that in the other stages, while the average relative abundance of Actinobacteria displayed an opposite trend ( Figure S1). Fungal phyla's relative abundances changed considerably in the course of P. armeniacum development. The average relative abundance of Basidiomycota increased during P. armeniacum development, whereas the average relative abundance of Ascomycota decreased significantly in the course of development ( Figure S1).

Taxonomic Assignment of the Microbial Community Composition
In total, 3 and 16 phyla and 9 and 27 classes were identified in the fungal and bacterial communities, respectively. In the bacterial communities ( Figure 5A, Table S1), the dominant phyla (>1%) were Proteobacteria, Firmicutes, and Actinobacteria, with average relative abundances of 43.98%, 38.51%, and 15.81%, respectively. In the fungal communities ( Figure 5B, Table S1), the dominant phyla across all samples were Ascomycota and Basidiomycota, with average relative abundances of 62.94% and 18.60%, respectively. The average relative abundance of Proteobacteria in the seedling stage was higher than that in the other stages, while the average relative abundance of Actinobacteria displayed an opposite trend ( Figure S1). Fungal phyla's relative abundances changed considerably in the course of P. armeniacum development. The average relative abundance of Basidiomycota increased during P. armeniacum development, whereas the average relative abundance of Ascomycota decreased significantly in the course of development ( Figure S1).
A detailed analysis of the relative distributions of the sequences from the phylum to genus levels revealed significant variations in the course of Paphiopedilum armeniacum development (Table S2). Among the bacterial genera, the sequences assigned to Caldibacillus, Paenibacillus, Uruburuella, and Vibrio were significantly enriched in the seeding stage. Conversely, among the fungal genera, unclassified Tulasnellaceae showed a significantly increase during P. armeniacum development (Table S2).

Covariation Network Analysis of Core Bacterial Genera and Mycorrhizal Fungi
Co-occurrence networks can display interrelationships among various microorganisms and reveal effects of microorganisms overlooked in difference comparisons based on abundance [41]. The core microbiota potentially performs critical functions within its habitats. In the present study, there were obvious interactions between bacteria and fungi based on the covariation networks of core bacterial genera and mycorrhizal OTUs ( Figure  6). Mycobacterium was significantly correlated with seven genera, whereas Sphingomonas  (Table S2). Among the bacterial genera, the sequences assigned to Caldibacillus, Paenibacillus, Uruburuella, and Vibrio were significantly enriched in the seeding stage. Conversely, among the fungal genera, unclassified Tulasnellaceae showed a significantly increase during P. armeniacum development (Table S2).

Covariation Network Analysis of Core Bacterial Genera and Mycorrhizal Fungi
Co-occurrence networks can display interrelationships among various microorganisms and reveal effects of microorganisms overlooked in difference comparisons based on abundance [41]. The core microbiota potentially performs critical functions within its habitats. In the present study, there were obvious interactions between bacteria and fungi based on the covariation networks of core bacterial genera and mycorrhizal OTUs ( Figure 6). Mycobacterium was significantly correlated with seven genera, whereas Sphingomonas and Brevibacillus were significantly correlated with six genera. Excluding Cupriavidus, there were significant correlations among other core bacterial genera. OTU1, OTU3, OTU11, and OTU33 were significantly correlated with bacterial genera, with OTU33 being significantly correlated with six bacterial genera. Mycobacterium, Sphingomonas, Brevibacillus, and OTU33, the four core genera, interacted with each other and influenced seven bacterial genera and three OTUs of Tulasnella species. The covariation network analysis results showed significant interactions between core bacterial genera and OTUs of Tulasnella species.

Covariation Network Analysis of Core Bacterial Genera and Mycorrhizal Fungi
Co-occurrence networks can display interrelationships among various microorganisms and reveal effects of microorganisms overlooked in difference comparisons based on abundance [41]. The core microbiota potentially performs critical functions within its habitats. In the present study, there were obvious interactions between bacteria and fungi based on the covariation networks of core bacterial genera and mycorrhizal OTUs ( Figure  6). Mycobacterium was significantly correlated with seven genera, whereas Sphingomonas and Brevibacillus were significantly correlated with six genera. Excluding Cupriavidus, there were significant correlations among other core bacterial genera. OTU1, OTU3, OTU11, and OTU33 were significantly correlated with bacterial genera, with OTU33 being significantly correlated with six bacterial genera. Mycobacterium, Sphingomonas, Brevibacillus, and OTU33, the four core genera, interacted with each other and influenced seven bacterial genera and three OTUs of Tulasnella species. The covariation network analysis results showed significant interactions between core bacterial genera and OTUs of Tulasnella species.

Comparison of the Mycorrhizal Fungi of P. armeniacum in Xingyi Greenhouse and in Other Habitats
Seven Tulasnellaceae OTUs, designated OTU1, OTU3, OTU11, OTU33, OTU67, OTU74, and OTU178, which accounted for 39.54% of the total reads, were observed in the present study. The sequences were BLAST searched against the GenBank database (National Center for Biotechnology Information) for identification and the sequences deposited in the GenBank database (Accession No: OL891761-67) ( Table 1). Phylogenetic analysis showed that all the sequences were roughly divided into four groups (Figure 7). OTU1 (type1) and OTU3 (type7) appeared closely related to mycorrhizal fungi from different habitats of P. armeniacum observed in previous studies [14]. OTU11, OUT33, and FJ786646(type20) formed one branch; OTU178 resulted similar to OTU1, whereas OTU67 formed an independent branch. The seven OTUs were distributed distinctly among the different growth phases of P. armeniacum (Table 1). OTU1 was observed in all developmental stages of P. armeniacum and was the only mycorrhizal fungus in the seedling stage. OTU11 and OTU178 were only observed in the reproductive growth stage, whereas OTU74 was only observed in the vegetative growth stage.  Scale bar is shown to infer evolutionary distances. PAB represent Paphiopedillum armeniacum from Baoshan (wild); PAK represents P. armeniacum from Kunming (greenhouse). The sequences with NCBI serial numbers are mycorrhizal fungi sequences obtained from P. armeniacum by Huang [14] and typing of the sequences.

Discussion
In the present study, the authors investigated microbial community structure dynamics across three growth and development stages of cultivated P. armeniacum. In addition, Tulasnellaceae fungi of P. armeniacum at different growth stages were analyzed to Scale bar is shown to infer evolutionary distances. PAB represent Paphiopedillum armeniacum from Baoshan (wild); PAK represents P. armeniacum from Kunming (greenhouse). The sequences with NCBI serial numbers are mycorrhizal fungi sequences obtained from P. armeniacum by Huang [14] and typing of the sequences.

Discussion
In the present study, the authors investigated microbial community structure dynamics across three growth and development stages of cultivated P. armeniacum. In addition, Tulasnellaceae fungi of P. armeniacum at different growth stages were analyzed to gain insights into their mycorrhizal features. Most studies on the root biomes of orchids have mainly focused on specific stages of growth, with little attention directed at trends during the development of cultivated orchids [5,6,[42][43][44][45].
The composition of endophyte communities is governed by host genotype and developmental stage, as well as by the environment. However, the significate of the effects may differ between studies [42]. When compared to environmental factors, greater influence of the host plant was observed in root endophytes in some research. Xiong et al. demonstrated that microbiome assembly along the soil-plant continuum is shaped predominantly by compartment niche and host species rather than by site or fertilization practice [43]. Strong cultivar-dependent variations in the fungal and bacterial microbiome were found for Cannabis sativa [44,45]. In order to minimize the impact of the substrate on endophytes, the stability and proportion of the substrate were maintained as far as possible. Meanwhile, part of the original substrate was retained during the potting process in this study. Host effects have been reported in many studies. The common OTUs of different growth stages and the genus of the core bacteria indicated the influence of host selection by P. armeniacum to a certain extent.
Higher richness and diversity were observed in the seedling stage than in the other development stages for both bacterial and fungal datasets, indicating that P. armeniacum may rely on interactions with more fungi or bacteria during its early growth than in the later growth stages, which is similar to what observed for Gymnadenia conopsea, in which microbial diversity in root samples was lower in the reproductive growth stage than in the vegetative growth stage [46]. The PCoA analysis results revealed that microbial community taxa were clustered at the seedling stage; no obvious clustering was observed at the other two stages. The results suggest that the microbial communities recruited by P. armeniacum at the initial growth stages were distinct from those in the following two stages, and there was a certain degree of overlap in microbial communities during the vegetative growth stage and the reproductive growth stage; some of these microorganisms could participate in essential processes.
Many endophytic bacteria isolated from orchid tissues play important functions for their hosts, including plant growth promotion [25,[47][48][49][50]. However, little is known about orchid-associated microbial community structure in the course of growth and development of P. armeniacum. Across the three development stages of P. armeniacum defined in the present study, Proteobacteria, Firmicutes, and Actinobacteria were the dominant groups, which is consistent with findings reported for other orchids [51,52]. The bacterial genera that exhibited the greatest differences among the three development stages were Caldibacillus, Paenibacillus, Uruburuella, and Vibrio, which also represent the genera most commonly found as bacterial endophytes [53,54]. The genus Caldibacillus includes thermophilic and aerobic or facultative anaerobic bacilli, and some species in the genus exhibit cellulase activity [55]. Bacteria in the genus Paenibacillus, which can produce indole-3-acetic acid, are commonly reported to promote plant growth and increase root development [47,56]. In addition, some Vibrio strains have been reported to participate in nitrogen fixation or salicylic acid synthesis [57,58]. Changes in relative abundance at the genus level indicated that P. armeniacum may selectively modulate microbial communities to meet its requirement during growth.
Interactions between fungi and bacteria are beneficial for the growth of orchids [27,49,59]. The obvious interactions between OTUs from Tulasnellaceae and core genera indicated that fungi are also key factors influencing bacterial community composition, especially in the orchid microbiota. Whereas orchids fully or partially rely on mycorrhizal fungi, some bacteria, such as MHB, commonly occur in ectomycorrhiza and in arbuscular mycorrhizal associations, which facilitate mycorrhiza formation [60,61]. In the co-occurrence network analysis in the present study, Brevibacillus, Mycobacterium, and Sphingomonas were core bacteria. In previous studies, Brevibacillus spp. were considered an MHB, whereas Mycobacterium spp. and Sphingomonas spp. were reported to have plant growth promotion effects [62][63][64].
As show in the results, a lot of non-mycorrhizal microorganisms were detected in the roots of P. armeniacum (Tables S1 and S2). The composition of non-mycorrhizal fungi, which were highly abundant in the seedling stage, began to converge throughout the vegetative stage and stabilized during the reproductive phase. This indicated that P. armeniacum screened for non-mycorrhizal fungi depending on the development of growth. Some non-mycorrhizal microorganisms are helpful for seed germination or seedling growth of orchids [9,21]. Fusarium, which was found in the root of P. armeniacum, is reported to form hyphae inside the cells and presented germinative potential, favoring seed germination or seedling growth of orchids [65,66]. Whether these strains are strictly mycorrhizal fungi requires further verification by methods such as stable isotope analyses.
Tulasnellaceae is considered a dominant group among the mycorrhizal fungi of orchids [67,68]. Researchers identified Tulasnellaceae through uncultured techniques as potentially mycorrhizal in orchids. In the present study, Tulasnellaceae were the only dominant fungal symbionts in P. armeniacum and were associated with the three growth and development stages. The compositions and levels of specificity of fungal partners can vary [69][70][71]. In the present study, different growth stages of P. armeniacum harbored different types and amounts of Tulasnellaceae OTUs. The seedling stage was only associated with one Tulasnellaceae OTU (OTU1), while five and six Tulasnellaceae OTUs were observed in the vegetative growth and reproductive growth stages, respectively. Among the seven Tulasnellaceae OTUs detected in the present, only OTU1 (Type7) was observed in all three development stages. Considering the findings of previous studies on mycorrhizal fungi of P. armeniacum [12][13][14], OTU1 may be the core group in P. armeniacum mycorrhizal fungi in different growth environments, geographical locations, and developmental stages and is also the first group to be recruited from the environment. OTU3 (Type25) is the main mycorrhizal group that was present in other habitats. The closest matches in GenBank for OTU11 and OTU33 were isolated from Paphiopedilum micranthum in Xingyi greenhouse, which indicated that OTU11 and OTU33 in P. armeniacum might be derived from horizontal transfer in the environment. OTU11, OTU74, and OTU178 were only associated with one stage of development. OTU67 and OTU178 may be novel taxa recruited from the environment, or they may not have been detected before due to the limited techniques used in previous research. Mycorrhizal fungi were similar in wild and cultivated P. armeniacum, indicating that the cultivated P. armeniacum may survive when transplant in the wild. Meanwhile, the similarity in mycorrhizal fungi makes ex situ conservation or even propagation by means of mycorrhization of axenically grown seedlings possible. Considering the important roles of bacteria and fungi in orchid growth and development [25,49], a more comprehensive study should be carried out in the wild P. armeniacum.
Due to overharvest, habitat loss, and degradation, wild populations of P. armeniacum have declined drastically. Conservation and restoration of the populations of this endangered orchid could facilitate the maintenance of its biodiversity as well as ecosystem stability. It is essential to comprehensively consider the structures of the resident microbial communities in this endangered orchid [72][73][74]. Meanwhile, although Illumina-based sequencing has some certain benefits for understanding the non-culturable endophytes of P. armeniacum, culture-depended methods are also necessary. In the future, more culturomics research is needed to provide strain resources in ecology and conservation research and practice.

Conclusions
In the present study, the diversity of microbial communities at different cultivated P. armeniacum development stages was examined based on Illumina sequencing approaches. Microbial community diversity and composition varied across the development stages. The seedling stage presented the highest microbial community richness and diversity, which indicated that P. armeniacum may need the interaction with more fungi or bacteria during the seedling stage. In addition, the core bacterial genera and the main OTUs of Tulasnellaceae showed significant interactions, according to the covariation network analysis results. Tulasnellaceae were the only dominant fungal symbionts during P. armeniacum growth. The results of the present study could facilitate the ex situ conservation and commercial development of endangered P. armeniacum populations.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/d14050321/s1, Fasta S1: Sequences of Tulasnellaceae displayed in phylogenetic trees; Table S1: All OTU sequences, profiling, and annotations. Table S2: Significant variations in genera in the course of Paphiopedilum armeniacum development; Figure S1: Shifts in microbial (bacterial and fungal) taxonomic composition at the main phylum level in the course of P. armeniacum development.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used in the study is available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.