Dynamic distribution of Massilia spp. in sewage, substrate, plant rhizosphere/phyllosphere and air of constructed wetland ecosystem

Introduction Massilia bacteria are widely distributed and have various ecological functions. Preliminary studies have shown that Massilia is the dominant species in constructed wetland ecosystems, but its species composition and distribution in constructed wetlands are still unclear. Methods In this paper, the in-house-designed primers were used to construct a 16S rDNA clone library of Massilia. The RFLP sequence analysis method was used to analyze the diversity of Massilia clone library and the composition of Massilia in sewage, substrate, plant rhizosphere, plant phyllosphere and air in a constructed wetland sewage treatment system. Redundancy analysis (RDA) and canonical correspondence analysis (CCA) were used to analyze the correlation between environmental factors and the population characteristics of Massilia in the corresponding environment. The dominant species of Massilia were analyzed for differences. Results The results showed that the 16S rDNA clone library in primer 5 worked well. According to the clone library diversity index analysis, the richness of Massilia varied significantly in different environments in different seasons, where the overall summer and autumn richness was higher than that in the spring and winter. The relative abundance of 5 Massilia in the constructed wetland ecosystem was greater than 1% in all samples, which were M. alkalitolerans, M. albidiflava, M. aurea, M. brevitalea, and M. timonae. The seasonal variation of dominant genera was significantly correlated with environmental factors in constructed wetlands. Discussion The above results indicated that the species of Massilia were abundant and widely distributed in the constructed wetland ecosystem, and there were significant seasonal differences. In addition, the Massilia clone library of constructed wetland was constructed for the first time in this study and the valuable data of Massilia community structure were provided, which was conducive to the further study of microbial community in constructed wetland.


Introduction
Massilia spp. are Gram-negative bacteria that are widely distributed. For example, four species of Massilia were isolated from Antarctic streams, lakes and weathering layers (Holochová et al., 2020). Massilia aquatica was isolated from a subtropical stream in China (Lu et al., 2020). Massilia sp. YMA4 was isolated from the open ocean off Lamay island in Taiwan (Ho et al., 2021), and M. buxea was found on rock surfaces (Sun et al., 2017). In addition, Massilia spp. can also be found in water, organisms and human clinical specimens. Among those, soil is the main environment of Massilia spp (Hye et al., 2022).
Numerous studies have shown that Massilia spp. have many functions. Some Massilia are able to synthesize violacein (Agematu et al., 2011;Baek et al., 2014), which has a good therapeutic effect on tumors, malaria and diarrhea and has great application potential in the pharmaceutical field. PHAs (as reserve material displayed as intracellular bright granules) is a kind of natural polymer biomaterial, which has good biocompatibility and biodegradability, and has become the most active research focus in the field of biomaterials in recent years. For the first time, Cerrone et al. find majority of the Massilia strains were capable to produce significant amounts of PHA under batch culture mode using soluble starch as carbon source (Cerrone et al., 2011). Studies have shown that some Massilia spp. were beneficial to environmental modification. For example, M. aromaticvorans ML15P13 T can degrade toxic gaseous organic air pollutants, such as benzene, toluene, ethylbenzene and xylene isomers (m, p and o-x) (BTEX) (Son et al., 2021); M. tieshanensis TS3 and M. putida 6NM-7 are resistant to multiple heavy metals (Du et al., 2012;Feng et al., 2016); M. phosphatilytica 12-OD1 has phosphate solubilization and can improve available phosphorus content in soil and improve the soil environment (Zheng et al., 2017); Massilia sp. WF1 and Massilia sp. WG5 can degrade phenanthrene and repair contaminated soil (Lou et al., 2016;Gu et al., 2021); and M. chloroacetimidivorans TA-C7e degrades the chloroacetamide herbicide (Lee et al., 2017). Some Massilia bacteria can also promote plant growth and increase the survival rate of plants in harsh environments. Krishnamoorthy et al. found that the combined application of Massilia and arbuscular mycorrhizal fungi alleviated the effects of salt stress on plant growth, root colonization, nutrient accumulation and growth in studying the effects of three different salinity on the growth of corn plants in coastal reclamation area (Krishnamoorthy et al., 2016). Turnbull et al. isolated a strain of Massilia from potato roots that can produce heteroauxin and cellulosedegrading enzymes to promote plant growth (Turnbull et al., 2012).
In recent years, research on Massilia spp. made significant progress. However, the community composition and distribution of Massilia spp. in the constructed wetland sewage treatment system has been neglected. Preliminary research on this research group showed that Massilia spp. are the dominant bacteria in the air of constructed wetland sewage treatment systems. The pollutant degradation ability of artificial wetland sewage treatment systems is higher than that in natural wetlands and microbes play an important role (Liu et al., 2019;Hu et al., 2022). In consideration of the functions of Massilia spp., it is of great significance to study the community composition of Massilia bacteria in constructed wetland ecosystems. The culturable microorganism accounts for 1% ~ 5%; therefore, molecular biology has been widely used in microbial community research. Lang et al. and Zhang et al. studied the community structure of constructed wetlands using molecular biology methods Lang et al., 2021). This study was the first to investigate the composition and distribution characteristics of Massilia in sewage, substrate, plant rhizosphere, plant phyllosphere and air in a constructed wetland sewage treatment system by two molecular biology methods: construction of a 16S rDNA clone library and RFLP sequence analysis. At the same time, redundancy analysis (RDA) and canonical correspondence analysis (CCA) were used to analyze the correlation between environmental factors and population characteristics of Masilia. The compositional differences and dynamic changes of Massilia in different environments and seasons were also explored, which is of great significance to reveal the mechanism of degradation of pollutants in the constructed wetlands and is beneficial to the operation of constructed wetlands.

Sampling sites
The free water surface constructed wetland, located in the interior of the Huangdao district (Qingdao city, Shandong province, China), at a latitude of 35°35′ to 36°08' North and a longitude of 119°30′ to 120°11′ East, is a part of an integrated sewage purification system ( Figure 1). This region has a warm temperate continental monsoon climate with a mean annual temperature of 12°C and a mean annual precipitation of 794 mm. The constructed wetland wastewater treatment system had a total area of 76.7 hm 2 and a treatment capability of 3.0 × 10 4 m 3 ·d −1 and was surrounded by the Yellow Sea to the east and south. It consisted of 99 treatment beds and received secondary unchlorinated wastewater from the Jiaonan Municipal Wastewater Treatment Facility with an A 2 O application as the secondary treatment. All beds were planted with common reed (Phragmites australis) and a number of naturally germinated wetland plants (Typha orientalis, Scirpus validus, and Lemna minor). To facilitate the harvest progress of aboveground biomass, sewage did not enter the constructed wetland bed from December to March of the following year.
A total of five sampling points were selected in this experiment, as shown in Figure 2. The soil, wastewater, leaves and air were sampled from April 2018 to January 2019, and the sampling collection frequency was once every month.

Sampling method
Fifty grams of soil sample and 10 L of water sample were collected from each sample site using sterile sealed bags and sterile bottles, respectively. After removing the fine roots in the soil samples, the water and soil samples were transferred to the laboratory immediately. After samples were dewatered by centrifugation, a fraction of the soil samples was stored at −20°C for molecular analysis.
The rhizosphere samples of 3 reed plants were randomly selected at each sampling point. Large soil was removed from the roots, the roots were soaked in sterile water, and the difficult soil on the roots was rinsed with sterile water. The collected water samples were mixed evenly, concentrated by centrifugation and placed in a sterile centrifuge tube to extract DNA.  Schematic diagram of the geographic location of the study area.

FIGURE 2
Schematic diagram of the sampling location of the constructed wetland. For the specific location of 5 sampling points, air samples, leaf interval, sewage, substrate and rhizome samples.
Frontiers in Microbiology 04 frontiersin.org Leaves of the reed plants were aseptically collected from different locations within the experimental unit and placed in Ziploc bags or Whirl-Pak bags by using new gloves for each replicate and by applying an ethanol disinfection of the pruning shears between samples. Samples were then transported back to the laboratory at 4°C. One hundred milliliters of sterile water was added to the bags, and the samples were agitated for 1 min by hand and then sonicated for 2 min. The microfloral wash was then transferred to polypropylene tubes and centrifuged at 30,000 ×g overnight at 4°C. The pellet was then transferred to a microcentrifuge tube and stored at −80°C until DNA extraction was performed.
The bacterial activity samples in TSP were collected via two KC-1000 high-flow air samplers with fixed Teflon membranes. The samples were placed 1.5 m from the ground, and the TSP samples were collected on the membranes at a sampling flow rate of 1250 L/min. After 24 h of sample collection, the membranes of the samples were quickly sent back to the laboratory and stored at −80°C until subsequent analysis.

DNA extraction
Sterile distilled water was used to rinse the surface of the Teflon membrane to obtain the particulate matter. After centrifugation, the particulate matter was transferred to a sterile 2 mL centrifuge tube. The sample quality was maintained at 0.2 ~ 0.3 g. Using the E.Z.N.A. VR Soil DNA Kit (Omega Bio-Tek, USA), DNA samples were obtained according to the manufacturer's instructions. An ultramicro spectrophotometer (IMPLEN, Germany) was used to detect the concentration.

PCR amplification of the gene fragment
Five in-house-designed primers were amplified for 16S rDNA, and the community diversity of Massilia was analyzed by RFLP (Table 1). The PCR amplification procedure was as follows: 94°C for 5 min; 94°C for 30 s, annealing temperature Tm°C for 30 s, and 72°C for 1 min (36 cycles), and 72°C extension for 10 min.
2.3.3. Construction of a 16S rDNA clone library and RFLP sequence analysis and sequencing PCR products were purified by a PCR product Purification Kit (OMEGA). Purified PCR products were inserted into the T-vectors by a T4 ligase cloning kit (MBI Fermen-tas). The recombinant was transformed into Escherichia coli (Top10) cells. The cells were cultured on LB plate medium with X-Gal/IPTG resistance screening, and the white bacterial colonies were selected and a cloned library was constructed. The vector primers M13 were used for the sequencing reactions according to the manufacturer's instructions. The new PCR products were digested by the restriction endonucleases HhaI and RsaI, and the digested products were detected by agarose gel electrophoresis. The DNA fingerprints were analyzed, and the types were designated.
After analyzing the enzyme digestion map, an OTU clone was selected to prepare the puncture tube, which was cultured for 16 h and sent to Shenggong Bioengineering (Shanghai) Co., LTD for sequencing. The sequencing results were inputed into the BLAST program in NCBI to compare with the sequences in the database, and the coverage (C), Shannon diversity index (H′), Simpson index (D) and richness (R) of the clone library were calculated. Physicochemical indexes of sewage: Ammonia nitrogen (NH 3 -N) was determined by Nessler 's reagent spectrophotometry. Nitrite nitrogen (NO 2 − -N) was determined by N-(1-naphthyl) ethylenediamine spectrophotometry. Nitrate nitrogen (NO 3 − -N) was determined by ultraviolet spectrophotometry. The pH, DO, water temperature and redox potential were measured by a multi-parameter water quality analyzer (HQ30d, HACH, USA).
The soil environmental factors of constructed wetland: pH, Organic matter, Water content, TP, TN, TK. The physical and chemical indexes of the substrate: The pH was determined by HI2221 pH meter. Organic matter was determined by potassium dichromate volumetric method. The water content was measured by drying at 105°C to constant weight. Total phosphorus (TP) was determined by molybdenum antimony colorimetric method. Total nitrogen (TN) was determined by micro-Kjeldahl method. Total potassium (TK) was determined by flame photometric method.
The ambient temperature and humidity of the sampling points were measured by TASI-620 digital thermometer and hygrometer. The wind speed was measured by TASI-641 anemometer. Gaseous pollutants O 3 , SO 2 , NO and NO 2 were measured in real time online by 49C ozone analyzer, 43C sulfur dioxide analyzer and 42C nitrogen oxide analyzer of American thermoelectric company, and CO was monitored in real time by 300 EU CO analyzer.
The above environmental factors were measured using different instruments and methods, as shown in Tables 2-4. 2.3.5. The structural differential analysis method of the Massilia species group OriginPro 9.1 was used to draw the analysis map differences of Massilia species number and dominance, and the R language tool was used to draw and analyze the Venn diagram of Massilia's population structure correlation.
Redundancy analysis (RDA) and canonical correspondence analysis (CCA) were used to analyze the correlation between environmental factors and the population characteristics of Massilia in the corresponding environment. Canoco 5 software was used to analyze the discriminant components of the data. If the length of the first axis was greater than 4.0, CCA should be selected. If between 3.0 and 4.0, RDA and CCA can be selected; if less than 3.0, RDA results are better than CCA. Finally, the RDA or CCA analysis chart is drawn. The contribution rate of environmental factors to the change of population structure was obtained by VPA analysis of R language. The experimental data were analyzed by SPSS 20.0 software.

Effects of different pairs of primers on the 16S rDNA amplification of Massilia spp.
The 16S rDNA amplification products of Massilia 16S rDNA amplified by different primers were detected by 1% agarose gel electrophoresis ( Figure 3 in the Additional files). The result shows that Primers 2, 3 and 4 all had nonspecific amplication production. Primers 1 and 5 had better amplification effects and can be used for 16S rDNA clone library construction.   Electrophoresis analysis of PCR products with the five different primers.

Analysis of the enzyme digestion products for the amplification of two pairs of primers
The 16S rDNA of Massilia spp. was amplified by using primers No. 5 and No. 1. Then, the amplified products were digested with the endonucleases RsaI and HhaI, and the digested products were detected by 3% agarose gel electrophoresis. The results shows that the digested products of 16S rDNA in the clone library amplified by primer No. 1 had similar fragments ( Figure 4A in the Additional files), while the digested products of 16S rDNA amplified by primer No. 5 had discrepancies ( Figure 4B). Figure 5 shows the OTUs in the constructed 16S rDNA clone libraries amplified by primers No. 1 and No. 5. There were 6 OTUs in the clone library amplified by primer No. 1, while there were 24 OTUs in the clone library amplified by primer No. 5. This indicated that the 16S rDNA amplified products of primer No. 5 were better in distinguishing the diversity of Massilia spp. Therefore, primer No. 5 was chosen for 16S rDNA clone library construction of Massilia spp.

Diversity analysis of the Massilia spp. in constructed wetlands
The specific alpha diversity indices are shown in Table 5. Through the coverage index, the coverage of the 16S rDNA clone libraries of Massilia spp. in the constructed wetland was greater than 86%, which indicated that most of the Massilia spp. were included in the constructed clone libraries. Large differences in the richness of Massilia spp. could be observed every month. The sampling in summer had the highest abundance, while the lowest abundance was in the samples collected in the winter. The richness of Massilia spp. was obviously different in the constructed wetland in different seasons and was higher in the summer and autumn than in the spring and winter. The Shannon index and Simpson index were selected to estimate the diversity of Massilia spp. According to the Shannon index and Simpson index, the samples in the summer had the highest community diversity.
3.4. The community structure of Massilia spp. in a constructed wetland sewage treatment system 3.4.1. The community structure of Massilia spp. in sewage of constructed wetland A total of 23 Massilia species were found in constructed wetland sewage ( Figure 6). There were significant differences in the number of dominant species of Massilia spp. in different months. There was a significant temporal variation in the distribution of dominant Massilia spp. The relative abundance of M. albidiflava was the highest in November. From June to August, the relative abundance was less than 1%, which was different from that in other seasons. The relative abundance of M. alkalitolerans was the highest in September, which was no more than 0.5% relative abundance from June to August. The relative abundance of M. aurea was more than 7% over the year, and there was no significant difference in different months. The relative abundance of M. brevitalea was highest in the summer and lowest in the autumn, which was similar to the variation in M. brevitalea in the substrate and rhizosphere. The relative abundance of M. timonae was highest in the summer, with approximately 20% relative abundance in other seasons. The relative abundance of M. umbonata was lowest in May and highest in December, and the seasonal variation was insignificant.
The data were analyzed with Canoco 5 software. The findings indicated that the RDA linear model should be chosen for analysis because the length of the first axis was less than 3.0. The findings reveal that the first axis' interpretation rate is 40.35%, the second axis' interpretation rate is 23.63%, and the combined interpretation rate of the two axes is 63.9%. It demonstrates that the structure of Massilia species is significantly influenced by environmental factors in constructed wetland sewage. According to Figure 7 and Table 6, the dominant species M. timonae was highly significantly positively correlated with T (p < 0.01), and highly significantly negatively correlated with NH 3 -N and NO 3 − -N (p < 0.01). M. brevitalea was highly significantly positively correlated with NH 3 -N, NO 3 − -N, DO and φ (p < 0.01). M. albiflava was highly significantly negatively correlated with T (p < 0.01). M. aurea was significantly positively correlated with φ (p < 0.05). M. plicata and M. aerilata were significantly negatively correlated with NH 3 -N and φ (p < 0.05). The results of the VPA analysis showed that the four factors with the highest contribution rates were NH 3 -N, T, DO and NO 3 − -N, which were 26.5, 22.2, 21.8 and 21.8%. The findings revealed that the aforementioned four factors had the greatest influence on the composition and distribution of the Massilia community in constructed wetland sewage.
3.4.2. The community structure of Massilia spp. in constructed wetland substrate A total of 24 Massilia species were found in the constructed wetland substrate (Figure 8). The dominant species in the constructed wetland substrate were M. albidiflava, M. alkalitolerans, M. aurea, M. brevitalea, M. timonae and M. umbonata. The relative abundance of M. albidiflava was the lowest in the summer, which was significantly different from other seasons. There was no significant difference in the relative abundance of M. albidiflava in the other seasons. The relative abundance of M. alkalitolerans in the substrate decreased in order from autumn to winter to spring to summer. The relative abundance of M. aurea was higher in the spring and summer than in the autumn. The relative abundance of M. brevitalea was highest in the spring and winter and lowest in the autumn, and there were obvious seasonal differences. The relative abundance of M. timonae in the summer was significantly higher than that in the other three seasons. There was no significant difference in the relative abundance of M. timonae between months in the same season. The relative abundance of M. umbonata was lower all year.   The RDA linear model was used to analyze the data. The results show that the interpretation rate of the first axis is 31.57%, the second axis is 21.18%, and the common interpretation rate of the two axes is 52.75%. It shows that environmental factors in the constructed wetland substrate have a significant effect on the structure of Massilia spp. According to Figure 9 and Table 7, the dominant species M. timonae was very significantly positively correlated with T and TP content (p < 0.01). M. alkalitolerans was significantly negatively correlated with TP (p < 0.05). M. brevitalea was highly significantly positively correlated with TN (p < 0.01). M. albiflava was very significantly positively correlated with TN (p < 0.01), and highly significantly negatively correlated with TP and T (p < 0.01). M. aurea was significantly positively correlated with TP (p < 0.05). M. plicata and M. aerilata were positively correlated with water content and T. The results of the VPA analysis showed that the three factors with the highest contribution rates were T, TP and TN, which were 19.1, 14.1 and 9.7%. The findings revealed that the aforementioned three factors had the greatest influence on the composition and distribution of the Massilia community in constructed wetland substrate. The relative abundance of M. alkalitolerans was higher in the autumn and winter than in the spring and summer, and the relative abundance was lowest in the summer (less than 1%). The variation in relative abundance for M. aurea was relatively slight. The relative abundance of M. aurea was higher in spring and summer than in autumn and winter, and the relative abundance was lowest in January and highest in May. The relative abundance of M. brevitalea was higher in the spring and summer than in the autumn and winter, and the relative abundance was highest in the summer, at nearly 50%. In the autumn, the relative abundance of M. brevitalea was lower in September and October, at approximately 8%. The relative abundance of M. timonae was the highest from June to August, up to approximately 47%. The relative abundance in the Analysis of the 16S rDNA clone library based on microbes in sewage.   Analysis of the 16S rDNA clone library based on soil microbes.
other months was approximately 20%. The variation in the relative abundance of M. timonae in the rhizosphere sample of plants was similar to that in the substrate sample of constructed wetland. The relative abundance of M. umbonata was lower than 2.5% all year, and was higher in the spring and winter than in the summer and autumn.

The community structure of Massilia spp. in the phyllosphere sample of constructed wetland plants
A total of 18 Massilia species were found in the phyllosphere sample of constructed wetland plants (Figure 11). The relative abundance of M. albidiflava in the phyllosphere sample decreased in an order from autumn to spring to summer. The relative abundance of M. albidiflava was highest in the autumn and lowest in the summer. The variation regularity of the relative abundances of M. alkalitolerans was similar to that of M. albidiflava in the phyllosphere sample. There was no obvious seasonal variation in the relative abundance of M. aurea. The relative abundance of M. brevitalea exhibited significant differences between seasons, with the highest relative abundance observed in the spring and the lowest relative abundance observed in Redundancy analysis of Massilia and environmental factors in soil. Microbiology  11 frontiersin.org autumn. There were obvious seasonal differences in the relative abundance of M. timonae; the highest relative abundance was approximately 50% in the summer, and the lowest relative abundance was approximately 19% in the spring and autumn. The relative abundance of M. umbonata was similar to that in the other environments of constructed wetlands, and its relative abundance was lower.

Frontiers in
3.4.5. The community structure of Massilia spp. in the air sample of constructed wetland plants A total of 16 Massilia species were found in the air sample of the constructed wetland (Figure 12). The relative abundances of M. albidiflava and M. alkalitolerans exhibited significant seasonal variation; the relative abundance was higher in the autumn and winter, while the species were not detected in the summer. The variation in the relative abundance of M. aurea between seasons was not obvious; the relative abundance was highest in April, May, July and December (approximately 16%), while it was lowest in January. The relative abundance of M. brevitalea was similar to that in other environments of constructed wetlands, with the highest relative abundance observed in January and the lowest relative abundance observed in September and October. The relative abundance of M. timonae exhibited obvious seasonal variation, with the highest relative abundance observed in the summer, while there were no significant differences among the spring, autumn and winter.
The data in Table 4 and Figure 12 were examined using the RDA model. The results show that the first axis interpretation rate is 48.03%, the second axis interpretation rate is 29.67%, and the combined interpretation rate of the two axes is 77.69%. It was discovered that the air environmental factors of the constructed wetland had a substantial impact on the community structure of Massilia. As shown in Figure 13 and Table 8, the dominant species of the genus M. timonae in the air of constructed wetlands was highly significantly positively correlated with temperature, humidity and O 3 (p < 0.01), highly significantly negatively correlated with SO 2 , NO 2 , CO, PM 2.5 and PM 10 (p < 0.01), and significantly negatively correlated with wind speed (p < 0.05); M. alkalitolerans showed a highly significant positive correlation with wind speed (p < 0.01), a significant positive correlation with NO 2 (p < 0.05), and a highly significant negative correlation with humidity (p < 0.01); M. albidiflava showed highly significant positive correlations (p < 0.01) with SO 2 , NO 2  Analysis of the 16S rDNA clone library based on rhizosphere microbes.
Frontiers in Microbiology 12 frontiersin.org Analysis of the 16S rDNA clone library based on phyllosphere microbes. and PM 10 , significant positive correlations (p < 0.05) with PM 2.5 and CO, and highly significant strong negative correlations (p < 0.01) with temperature and humidity; M. aurea showed highly significant negative correlations (p < 0.01) with wind speed and PM 2.5 and significant negative correlations (p < 0.05) with CO and PM 10 . The results obtained from the VPA analysis showed that humidity, NO 2 , temperature and SO 2 were the four factors that explained the greatest contribution with values of 40.1, 11.4, 15.2 and 5.1%. The findings revealed that the community composition and distribution of Massilia in the air of the constructed wetland is most influenced by the above four factors.
3.5. The influence of environmental factors on the structure of the Massilia spp. community in constructed wetlands The structure of Massilia spp. in different seasons and different environments was analyzed by a Venn diagram (Figure 14). The species numbers of Massilia spp. in the constructed wetland sewage, substrate, plant rhizosphere, plant phyllosphere and air samples were 15, 16, 16, 10 and 8 species, respectively in the spring, 17, 20, 21, 14 and 7 species, respectively in the summer, and 13, 19, 19, 13 and 7 species, respectively Analysis of the 16S rDNA clone library based on airborne microbes. Microbiology  13 frontiersin.org in the autumn. The species number of Massilia spp. in the substrate and rhizosphere in the winter was 18, while it was 7 in the air. The above results show that the species of Massilia in the air samples were all found in the other environments of the constructed wetland. Therefore, it can be concluded that Massilia spp. in the air were probably from the underlying surface of the constructed wetland. PCA was employed to analyze the composition of Massilia species in each environmental sample of the constructed wetland, and the distribution characteristics among the samples were described by two-dimensional images, as shown in Figure 15. The results showed that there were significant differences in the composition of Massilia spp. among different environmental samples, while the differences in different seasons were small. Most interestingly, the study found that the similarity of Massilia spp. in different environmental samples within the same season was higher than that in the same environmental sample within different seasons.

Discussion
Since Massilia was first discovered by La Scola in 1998, scholars from all over the world have successfully discovered other Massilia from different regions and environments and studied them. In 2003, Wery et al. studied the microbial community in Antarctic soil, isolated a strain capable of producing protease, used universal primers to amplify 16S rDNA by PCR, and determined the phylogenetic relationship of the strain by gene sequencing, and classified it as belonging to the Massilia genus (Wery et al., 2003). Weon et al. and Feng et al. also used the same molecular biology approach to study Massilia found in other environments (Weon et al., 2009;Feng et al., 2016). It can be seen that there are many studies on Massilia, but its community structure composition and dynamic changes in different environments are rarely reported. This study was the first to systematically study the composition and distribution of Massilia spp. in constructed wetlands using the method of cloning library construction. Five pairs of primers designed in-house were innovatively used to amplify the 16S rDNA of Massilia spp. and further analyzed using agarose gel electrophoresis. The detection showed that the amplification effect of primer No. 5 was the best, and the electrophoresis bands of the restriction map of RFLP were abundant. Therefore, primer No. 5 was finally selected to construct a Massilia clone library and conduct RFLP sequence analysis. The coverage rate of the clone library constructed by primer No. 5 was greater than 86%, indicating that the clone library constructed by the primer contained most of the Massilia spp. in the wetland, which could more fully reflect the diversity of the Massilia community. The construction of the clone library can obtain a more systematic and comprehensive understanding of the community structure of Massilia in the constructed wetland, which provides valuable experience for studying the community composition of Massilia in other environments and provides data support for the construction and operation of the constructed wetland.
Previous studies have shown that the higher the H′ is, the lower the D is, and the higher the microbial diversity. The Shannon index and Simpson index of Massilia in the constructed wetland were analyzed, and the results showed that the variation range of Shannon index and Simpson index was 2.97-5.97 and 0.86-1.00, which were the highest in summer and significantly higher than that in winter, and the difference between spring and autumn was not significant. According to the clone library diversity index analysis, the abundance of Massilia in different environments of constructed wetlands in different seasons was significantly different, and the overall abundance in the summer was higher than that winter, indicating that there is a significant seasonal fluctuation in the structure and diversity of microbial communities in wetland ecosystems. When Chazarenc studied the vertical flow constructed wetland in France, he found that the degree of bacterial activity varied with seasonal temperature, the highest in June, and the   lowest in late winter and autumn (Chazarenc et al., 2010), which was similar to the results of this study. This may be due to the low temperature in winter and low microbial activity, so the microbial abundance is low. According to the correlation analysis of environmental factors, T was the most significant factor affecting the community structure of Massilia in constructed wetlands. After the beginning of spring, with the increase of temperature, the microbial activity increases, and the reproduction speed also accelerates, so the number of microorganisms gradually increases, reaching the peak in summer. The temperature in autumn is also high, so the microorganisms continue to grow, and then the number of microorganisms starts to decrease after the temperature decreases in late autumn, until it drops to the lowest point in winter. Buckeridge et al. research found that the number of bacteria in the soil microbial community of arctic tundra decreased significantly from the peak in the late thawing period to the low point in spring, and then increased to a level similar to that in autumn in midsummer, which is different from this study (Buckeridge et al., 2013). The reason may be that the climate types and ecosystems of the two places are different. The Arctic has a long winter and freezes all the year round; Qingdao is located in the mid latitude region, with four distinct seasons and short winter. This indicates that the seasonal variation of the composition of Massilia in the constructed wetland may be related to the environment and climate of the area where the constructed wetland is located, and the microbial community structure of the constructed wetland in different areas is different. This study found that Massilia were distributed and abundant in sewage, substrate, plant rhizosphere, phyllosphere and air in constructed wetlands, but their proportions in each environment were different. A total of 24 species of Massilia were detected in the constructed wetland system, and the number of species was in the following order: substrate and rhizosphere samples > sewage samples > phyllosphere samples > air samples. Moreover, the species of Massilia in air samples were found in other environments of the constructed wetland, indicating that the substrate (soil) and plant rhizosphere of the constructed wetland were important sources of Massilia. It may diffuses into the air from the underlying surface of the constructed wetland. In addition, previous studies have reported that the microbial community in the rhizosphere of constructed wetland plants is the most abundant and that plant root exudates can have a positive effect on microorganisms indicating that the constructed wetland substrate (soil) was an important source of Massilia. It diffuses into the air from the underlying surface of the constructed wetland (Chaparro et al., 2013;Wu et al., 2017).
Through Venn diagram comparison, we found that Massilia and dominant bacteria were similar in sewage samples, substrate samples, plant rhizosphere and phyllosphere samples and showed that the number of Massilia species in plant rhizosphere samples was slightly greater than that in substrates. This may be due to artificial wetlands using substrates, plants, microorganisms of physical, chemical and biological sewage treatment triple synergy (Stottmeister et al., 2003), and aquatic plants as an indispensable part of the purification process. This increase in the type of bacteria in constructed wetland system substrates and plant rhizospheres enhances the growth of microbial communities to create a favorable environment (Cao et al., 2017); for example, a variety of useful compounds are released into the rhizosphere Frontiers in Microbiology 15 frontiersin.org (Olanrewaju et al., 2019). Massilia can secrete auxin, accelerate the growth of plants and increase the concentration of cations such as iron (Kuffner et al., 2010), calcium and magnesium in plant roots (Hrynkiewicz et al., 2010), improve the nutrient supply of plant roots, and improve the survival rate of plants under adversity. Massilia and aquatic plants complement each other, and the two act synergistically, so there are more species of Massilia in the rhizosphere than in the substrate. This study found that there are 6 dominant species of Massilia in different environmental samples of constructed wetlands, with obvious seasonal differences, namely, M. albidiflava, M. alkalitolerans, M. aurea, M. brevitalea, M. timonae and M. umbonata. This shows that the dominant species of Massilia have multiple ecological functions. Among them, M. timonae was the most important of the six dominant bacteria, with the highest relative abundance of 58.75%. Because of its significant positive correlation with T, the relative abundance of M. timonae was higher in summer and autumn, and lower in winter and spring. In 2009, Faramarzi et al. found that the most dominant species, M. timonae, could produce chitinase in a medium containing colloidal chitin as the only carbon and nitrogen source, which greatly reduced the production cost of chitinase (Faramarzi et al., 2009). The relative abundance of the second largest species M. brevitalea was higher in winter and spring, reaching 54.76%. The correlation analysis of environmental factors showed that it was significantly negatively correlated with T, so it was a low temperature resistant Massilia. In 2008, Zul et al. isolated a strain of M. brevitalea from the substrate and found that it can hydrolyze starch and casein and has salt tolerance, and in this study, it was found that its relative abundance was high in the winter (Zul et al., 2008). In addition, M. brevitalea was once isolated from organic fertilizers as the dominant strain of ammonification, indicating that it is more suitable for survival in environments with higher organic matter and nitrogen content, and its abundance was positively correlated with total nitrogen and organic matter. There is little difference between M. alkalitolerans and M. aurea, but M. alkalitolerans is slightly dominant, and the relative abundance can reach 44.12%, while the highest relative abundance of M. aurea is 17.24%. M. alkalitolerans has alkali resistance (Xu et al., 2005), and M. aurea can produce auxin (indoleacetic acid) (Gallego et al., 2006). The relative abundances of M. albidiflava and M. umbonata were low. M. umbonata can produce poly-β-hydroxybutyrate (Rodriguez-Diaz et al., 2014). M. albidiflava was positively correlated with total nitrogen, which may be related to its nitrogen fixation activity (Zhang et al., 2006).
In addition, the environmental factors in the constructed wetland have a significant impact on the composition, distribution and seasonal variation of the community structure of the Massilia dominant bacteria. In this study, we found that M. albidiflava and M. brevitalea, the dominant species of Massilia, were positively correlated with NH 3 -N, NO 2 − -N and NO 3 − -N, consistent with their having nitrate reductase activity; M. timonae and M. aurea were positively correlated with TP content, which proved that Massilia had phosphorus solubilization; M. aurea absorbed the fungal secretion; The photosynthetic products secreted by the fungus were absorbed by M. aurea, which promoted the mineralization and transformation of organic phosphorus in the substrate. The substrate with high phosphorus content was more suitable for its survival . When studying the correlation between microbial community structure and environmental factors in Poyang Lake wetland profile, it was found that water content and pH were the main factors affecting bacterial community structure (Ma et al., 2018), which was consistent with the results of this study. Since airborne bacteria are usually attached to the surface of atmospheric particles and can easily spread with the wind. Many air environment factors affect atmospheric particulate matter, which in turn affects the composition and function of the microbial community attached to its surface (Hwang and Park, 2014). Therefore, the air environmental factors investigated in this study, such as temperature, humidity, wind speed, PM 2.5 and PM 10 , all had an impact on the seasonal variation and community structure of the Massilia dominant strains (Maron et al., 2006;Cao et al., 2014). In summary, the study of Massilia dominant species in the constructed wetland will help us to analyze the main functions of the constructed wetland and adjust its operation status, which can make it treat sewage more efficiently. Although the genus Massilia has multiple environmental functions, clinical medical studies have found that infection with one of the dominant species, M. timonae, has the potential to lead to multiple human health-threatening diseases, such as septicemia, osteomyelitis, and encephalitis (Scola et al., 1998;Van Craenenbroeck et al., 2011), but the infection mechanism is unknown. Song recently isolated a Massilia oculi sp. nov. of type strain CCUG 43427 T from the eyes of an endophthalmitis patient, and analyzed its genome. It was found that it contains pathogenic genes, indicating that the strain has posed a threat to human health (Song et al., 2019). The above results indicate that some Massilia bacteria have evolved a wide range of pathogenicity, and the air mosaic bacteria in constructed wetlands can spread with the air to the surrounding environment, which is one of the important ways of Massilia transmission. Therefore, studying the community composition of Massilia bacteria is important for the safety evaluation of artificial wetlands, as well as for public health and ecological safety.

Conclusion
In conclusion, self-designed primers were used to construct a 16S rDNA clone library, and RFLP sequence analysis was used to determine the composition of the Massilia community in sewage, substrate, plant rhizosphere, plant leaf and air in an artificial wetland sewage treatment system. The overall community of Massilia bacteria in the artificial wetland sewage treatment system was higher in the summer and autumn than in the spring and winter, while the relative abundance of the dominant species had obvious seasonal differences. There were significant seasonal differences in the relative abundance of dominant species in sewage, substrate, plant rhizosphere and plant leaf sphere. Environmental factors in the sewage, substrate and air of the constructed wetland had a great influence on the composition and seasonal variation of Massilia's community structure. The diversity of Massilia bacteria in the artificial wetland air was less different than in the environmental samples, and the overlap rate of the total OTUs with the remaining samples was 100%, indicating that Massilia spp. in the air were probably from the underlying surface of the constructed wetland. In different seasons, the composition of Massilia spp. was significantly different in the same environment. In the same season, the similarity of Massilia spp. in different environments was high, and there was homology.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.