Characteristics and Metabolic Patterns of Soil Methanogenic Archaea Communities in the High Latitude Natural Wetlands of China

Background: Soil methanogenic microorganisms are one of the primary methane-producing microbes in wetlands. However, we still poorly understand the community characteristic and metabolic patterns of these microorganisms according to vegetation type and seasonal changes. Therefore, to better elucidate the effects of the vegetation type and seasonal factors on the methanogenic community structure and metabolic patterns, we detected the characteristics of the soil methanogenic mcrA gene from three types of natural wetlands in different seasons in the Xiaoxing'an Mountain region, China. Result: The results indicated that the distribution of Methanobacteriaceae (hydrogenotrophic methanogens) was higher in winter, while Methanosarcinaceae and Methanosaetaceae accounted for a higher proportion in summer. Hydrogenotrophic methanogenesis was the dominant trophic pattern in each wetland. The results of principal coordinate analysis and cluster analysis showed that the vegetation type considerably inuenced the methanogenic community composition. The methanogenic community structure in the Betula platyphylla – Larix gmelinii wetland was relatively different from the structure of the other two wetland types. Indicator species analysis further demonstrated that the corresponding species of indicator operational taxonomic units from the Alnus sibirica wetland and the Betula ovalifolia wetland were closer. Network analysis showed that cooperative and competitive relationships exist both within and between the same or different trophic methanogens. The core methanogens with higher abundance in each wetland were conducive to adaptation to environmental disturbances. Conclusions: This information is crucial for the assessment of metabolic patterns of soil methanogenic archaea and future uxes in the wetlands of the Xiaoxing'an Mountain region given their vulnerability. divergent

Forested wetlands are an important component of wetland environments. They are located at the interface between forest and water ecosystems, and have more complex service functions than either wetlands or forests [36,37]. The wetlands of Xiaoxing'an Mountain, China, whose total area is 1,069,600 m 2 , are representative of the main mountainous wetlands found in high latitude temperate areas, and play a critical role in the study of Chinese wetlands [38]. Wetlands in this area are diverse and include forested wetland, shrub wetland, herbaceous vegetation wetland, moss bog wetland, and shallow marshland, most of which are primarily forested wetland types. Currently, there is a lack of understanding regarding the indicative groups and potential metabolic patterns of the methanogenic community in the high latitude wetlands of China's forested areas.
In this study, three typical types of natural wetland along the forest-wetland transition zone of Xiaoxing'an Mountain, namely the Betula platyphylla -Larix gmelinii, Alnus sibirica, and Betula ovalifolia wetlands were selected for analysis. The protein encoding gene mcrA was chosen as a molecular marker using high-throughput sequencing to: (i) explore the distribution characteristics of soil core methanogenic communities in the three wetland types; (ii) elucidate the divergence of the methanogenic community among different vegetation types and predict the potentially methanogenic metabolic patterns in wetlands; and (iii) determine the indicator species of soil methanogens according to vegetation type and seasonal variability. Exploring the community distribution characteristics of methanogens in the transitional zone of forested wetland is important for the prediction of methane release during wetland vegetation succession in forest regions.

Richness and diversity indices
The α-diversity indices were analyzed based upon the mcrA gene data for the different vegetation types and seasons ( Table 1). The Good's coverage estimator ranged from 99.95% to 99.99% for all samples, indicating that the sequences obtained represent the dominant phylotypes. In total, 147-496 OTUs were detected in different treatments, which implies that methanogen species are distributed widely in these wetland soils. Overall, the numbers of methanogen OTUs in winter were generally slightly higher than those in summer, compared with corresponding vegetation type. Except for B. platyphylla -L. gmelinii wetland in winter, the number of methanogenic OTUs in the A. sibirica wetland was slightly above that found in the other wetlands, in either winter or summer.
Similar trends were observed in the results of the Shannon, ACE, and Chao indices, which also showed that, in each vegetation type, the diversity and richness of methanogens clearly change with season. In addition, the Simpson's evenness index revealed that the distribution of the methanogen community in the B. platyphylla -L. gmelinii wetland was very uneven compared with that of the other wetlands. A one-way ANOVA of each index demonstrated that the results of the diversity and richness indices (including the Shannon, ACE, and Chao indices) were not signi cantly different (Table 1), but the results of the Simpson index was signi cantly different across vegetation type (p=0.007).

Taxonomic characteristics of methanogenic communities
Using taxonomic analysis at the class level, all OTUs obtained were identi ed as belonging to four taxonomic classes known to be present in all wetlands (Fig. 1A). The dominant classes in each of the three wetland types were Methanobacteria, Methanomicrobia, and Thermoplasmata, all of which belong to Euryarchaeota. Many differences in the community characteristics of these dominant classes were present in each wetland. The proportion of Methanobacteria was higher in winter in all wetlands, whereas Methanomicrobia had a higher distribution in summer compared with the corresponding wetland. Furthermore, the proportion of the most dominant class Methanobacteria was lower in the B. platyphylla -L. gmelinii wetland than in the other wetlands. In this habitat, its proportion dropped to 41.3% in summer, during which the proportion of Methanomicrobia remarkably increased (Fig. S1). Moreover, the highest number of unclassi ed OTUs was detected in the B. platyphylla -L. gmelinii wetland in winter, while the proportion was lowest in the B. ovalifolia wetland in summer.
At the family level, the effective sequences could be mostly classi ed into eight taxonomic families (Fig. 1B). Methanobacteriaceae, Methanosarcinaceae, and Methanosaetaceae were the most abundant families in all soil samples. Methanobacteriaceae accounted for a higher proportion (54.2-78.1%) in winter, while the proportions of Methanosarcinaceae and Methanosaetaceae increased gradually in summer. Noticeably, the total number of effective sequences assigned to Methanosarcinaceae in the B. platyphylla -L. gmelinii wetland in summer was signi cantly higher than that of the other wetlands (Fig. S2). Additionally, several small parts of the sequences could not be identi ed at the family level, which implies that unknown methanogens may exist in the wetlands studied.
Despite the clear differences at class and family levels, heterogeneity need further to be analyzed at lower taxonomic levels. A Venn diagram was drawn to illustrate the common and unique methanogenic OTUs present in all three wetland types studied (Fig. 1C). The methanogenic communities of all soil samples shared only 47 common OTUs, most of which can be assigned to Methanobacteria (59.6%), Thermoplasmata (19.1%), and Methanomicrobia (14.9%). These 47 OTUs were considerably abundant, together taking up 34.5-77.7% of the relative abundance at the class level in the whole community of each wetland, while the relative abundance of these core OTUs in B. platyphylla -L. gmelinii wetland in summer was the highest ( Fig. 1D and Table S1). Furthermore, the Venn diagram showed that the B. platyphylla -L. gmelinii wetland had the highest number of unique OTUs in winter, while the B. ovalifolia wetland had only 86 unique OTUs in the summer. Overall, every wetland observed a higher number of unique methanogenic OTUs in winter than in summer.
Network analysis was perform to assess the relationships among methanogenic communities at OTU level (

Comparison of methanogenic community similarities
To analyze the similarity of the methanogenic communities, cluster analysis was applied based upon the Bray-Curtis dissimilarity distance test (Fig. S3). In the results, soil samples from the B. ovalifolia and A. sibirica wetlands were gathered in a cluster, which implies that the methanogenic communities of these wetland types are similar. Soil samples from the B. platyphylla -L. gmelinii wetland formed an independent branch, indicating that the methanogen community in this wetland is distinctly differentiated from the other cluster. According to PERMANOVA analysis, there was a signi cant difference between the two clusters (r 2 =0.384, p=0.001).
To further verify the results from cluster analysis, principal component analysis was performed using unifrac PCoA based on the OTUs of the mcrA gene (Fig. 3). The rst two PC axes were estimated on account of the results from the soil samples, with 44.79% and 16.32% for the variance, respectively. Two groups were detached from each other along the PC1 axis in the PCoA diagram. The B. platyphylla -L. gmelinii wetland samples from both seasons were clustered together, while the samples of the B. ovalifolia and A. sibirica wetlands were combined into a group and spread along the positive direction of the PC1 axis. The pattern of the PCoA further strengthened the result of the cluster analysis. At the same time, it can be seen that the distribution of the methanogen communities was primarily characterized according to vegetation type.
Considering the aforementioned results, LEfSe analysis was used to distinguish the potential discriminating taxa from the methanogenic communities in each wetland type. The larger differences in the taxa were depicted using a cladogram (Fig. 4). The LEfSe analysis showed that there were 29 signi cantly differential taxa across the three vegetation types (LDA score = 2). Most taxa were attributed to the A. sibirica wetland (14 taxa), followed by the B. ovalifolia wetland (eight taxa), and the B. platyphylla -L. gmelinii wetland (seven taxa). In this cladogram, most of the taxa in the A. sibirica wetland were no rank or unclassi ed except for the members of Methanobacteria. The taxa in the B. platyphylla -L. gmelinii wetland can be primarily assigned to class Thermoplasmata, while those in the B. ovalifolia wetland belonged to the Methanocellaceae and Methanosaetaceae families.

Responses of methanogenic indicator species to vegetation type and season
To determine the in uence of vegetation type and season on methanogens, indicator species analysis was performed for all OTUs. This analysis demonstrated that 28 OTUs were signi cantly associated with vegetation type and seasonal pattern, as listed in Table  2. Five OTUs represented the A. sibirica wetland in winter, and were assigned to the Methanosaetaceae and Methanobacteriaceae families, while two indicator OTUs were representative of the summer. Only two indicator OTUs were found in the winter B. platyphylla -L. gmelinii wetland data, whereas there were nine for the summer. The indicator OTUs that were associated with the B. ovalifolia wetland for the summer all belonged to the family Methanosaetaceae (three OTUs), while the seven associated with the winter were classi ed into Methanosaetaceae and Methanobacteriaceae.
Coevolution of the methanogens is clearly demonstrated by the similar tree topologies of the mcrA gene (Fig. 5). Based on phylogenetic analysis, the methanogenic species of different wetland soils were categorized into four clusters (I, II, III, and IV). Cluster I contains the representative species assigned to Methanobacterium lacus, Methanobacterium sp. MB1, and Methanobacterium sp. SMA-27. Cluster II represents an unclassi ed cluster of methanogens within Methanosarcinaceae. Cluster III includes Methanoregula boonei, Methanosphaerula palustris, and Methanoculleus marisnigri, and the unclassi ed Thermoplasmata cluster. Cluster IV included not only the divergent hydrogenotrophic Methanocellales but also the unclassi ed Methanosaetaceae and Methanosarcinales.

Discussion
In this study, we demonstrated that the structure and diversity of methanogenic communities in different wetland soils varies with season. Shifts in vegetation type and season correspond to variation in the indicator taxa and metabolic patterns of methanogenic communities. Here, we focused on the roles that vegetation type and seasonal mixing played in structuring methanogenic communities.
Forty-seven common OTUs, mostly assigned to the known classes of Methanobacteria (59.6%), Thermoplasmata (19.1%), and Methanomicrobia (14.9%), were detected in all soil samples in the Venn diagram. These common OTUs can be seen as the core microbiota. They seem to share functional similarity in different environments and to provide a robust microbial ecosystem under environmental uctuation [39,40]. We paid close attention to the core species because they were deemed to re ect the stability of soil communities in wetlands [41]. This result indicates that the methanogenic community has high potential to adjust in response to changes in vegetation type and season. Although there were only 47 core OTUs, the corresponding species accounted for a large proportion of methanogenic community composition. The similarity in the soil methanogenic populations found to exist in different wetland types and seasons signi es that these populations could be resilient and were able to quickly return to their original abundance distributions as soon as disturbances occurred in the environment. This also indicates that these methanogens in dynamic wetland environments were likely to generate the ability to adapt to such disturbances [42]. Although similar sequences were detected among different wetlands, each wetland has its own characteristic and distinctive OTU types, which suggests that methanogen communities were disturbed either by changes in vegetation type or by season.
Seasonal variations in climate may have a considerable effect on methanogenesis in wetland ecosystems and is also an underlying driving factor affecting the structure of methanogenic communities [32,43,44]. Methanogens are sensitive to seasonal variation due to changes in oxygen availability and temperature [45], and thus the abundance and composition of methanogens is very variable in wetland ecosystems. In the current study, the results showed that the distribution of Methanobacteriaceae (hydrogenotrophic methanogens) was higher in winter in each wetland [33,46], while Methanosarcinaceae and Methanosaetaceae accounted for a higher proportion in summer. The variation in methanogenic community structure also gave rise to a shift in potential methanogenic metabolic pattern. Previous studies have demonstrated that hydrogenotrophic methanogen can exist under psychrophilic conditions [47], which agrees with the ndings of this research. Season altered methanogenic metabolic pattern. It can be clearly seen that the hydrogenotrophic pathway played a fundamental role in methanogenesis in each wetland. This also indicated that Cluster I, phylogenetically assigned to Methanobacteriaceae and its closest relatives, was the major population in all three wetland sites. Meanwhile, although most of the results of the α-diversity indices were higher in winter than in summer, we found no statistically signi cant difference in the richness and diversity of the methanogenic communities among the different samples. Given the toxic effect of oxygen on methanogens, it is possible that snow and ice cover formed a diffusive barrier with the air in winter [42,48].
Wetland soils were therefore in anaerobic conditions, which was bene cial to raise the diversity of the soil methanogenic community.
In this study, more unique OTUs were detected in winter, which further supports this conclusion.
Wetland vegetation species could also be an important factor driving the variation in methanogen community composition and methanogenesis [49][50][51]. We found differences in the abundance distribution and composition of the methanogenic communities according to wetland type. The soil methanogen populations in each vegetation type were unique, while differences in the vegetation type resulted in the observed discrepancies in methanogenic composition. The higher relative abundance of Methanobacteriaceae and Methanosaetaceae in A. sibirica and B. ovalifolia wetlands presented an obvious difference in vegetation pattern compared with the B. platyphylla -L. gmelinii wetland, in which higher relative abundance of Methanosarcinaceae and Methanomassiliicoccaceae was found. Species of Methanosarcinaceae, many of which can possess both hydrogenotrophic and acetoclastic methanogenesis [52], had a markedly higher abundance in the B. platyphylla -L. gmelinii wetland during the summer than other wetlands. Zhang et al also con rmed that the number and percentage of methanogens in a constructed wetland was affected by the plant species present [49]. Different characteristics of plants, such as the root exudates, litters, and root porosity can induce competition in microbial communities [49,53]. Furthermore, the available electron donors and acceptors in the soil are in uenced by plants species, which causes changes in the structural variation of the methanogenic community and that of other microbes [53].
The results of β-diversity analysis, such as cluster analysis and PcoA, demonstrated that vegetation type was the primary factor affecting the variation observed in the composition of the methanogenic community. LEfSe analysis revealed differences in the methanogenic taxa from phylum to genus level according to wetland vegetation type. The result demonstrated that the number of Methanosaetaceae members in the B. ovalifolia wetland were signi cantly higher than those in the other wetlands, indicating that the aceticlastic methanogenesis pathway (which is regulated by family Methanosaetaceae) occupied a signi cant place in this wetland. Meanwhile, the order Methanobacteriales, which regulates the hydrogenotrophic methanogenesis pathway, had a signi cantly higher distribution in the A. sibirica wetland. In addition, there was a considerable number of unknown taxa in this wetland, which need to be further explored. The relative abundance of the Thermoplasmata group (seven taxa) in the B. platyphylla -L. gmelinii wetland was signi cantly higher than in the other wetlands. Methanomassiliicoccales, as a subgroup of Thermoplasmata, is phylogenetically distant from the other methanogens [54]. Also known as Methylotrophic methanogens, Methanomassiliicoccales can metabolize methanol [55,56], and mostly originated from animal intestinal and rumen tracts [57,58]. As supported by the ndings of this study, they are generally less distributed in environmental ecosystems due to limited methylic precursors [59,60]. Nevertheless, Methanomassiliicoccales or unclassi ed Thermoplasmata-like species can utilize the noncompetitive methyl compounds as a preferred substrate to participate in the synthesis of methane [61], when acetoclastic and hydrogenotrophic methanogenesis are hindered [62]. The contribution of Thermoplasmata with less distribution in B. platyphylla-L. gmelinii wetland soil as a characteristic metabolic pathway can be inestimable, compared with other wetlands. These results indicate that wetland vegetation type signi cantly in uenced the composition of soil methanogens, and further created the conditions for different biochemical metabolic patterns for methanogenesis.
To investigate correlations between key methanogenic species, a network analysis was built and was used to determine whether their relationships are based on sharing or competition within complex co-occurrence patterns. Network analysis indicated strong links between methanogens, with signi cant positive or negative associations (276 edges) among 51 OTUs (nodes) from eight families (p < 0.05). These OTUs were correlated with hydrogenotrophic (Methanobacteriaceae and Methanocellale), acetoclastic (Methanosaetaceae), facultative (Methanosarcinaceae), and methylotrophic (Methanomassiliicoccales and unclassi ed Thermoplasmata) metabolizing methanogens. In the network, hydrogenotrophic methanogenic lineages, such as members of the Methanobacteriaceae with 31 nodes (OTUs), were the most abundant, and strong positive and negative connections both within and between families were also present. OTU 1823, as the most abundant OTU in Methanosarcinaceae (facultative methanogenesis) was correlated negatively with some OTUs in Methanosaetaceae (acetoclastic methanogenesis) and a portion of OTUs in Methanobacteriaceae (hydrogenotrophic methanogenesis). Meanwhile, OTU 2174, OTU 1595, and OTU 72 in Methanosaetaceae also showed a negative correlation with some OTUs in Methanobacteriaceae. The previous study discovered that acetate could be degraded to produce CH 4 by either aceticlastic methanogenesis or by syntrophic acetate oxidation [25]. The aceticlastic pathway was only regulated by aceticlastic methanogens, whereas the syntrophic acetate oxidation pathway was mediated by the synergistic effect of two kinds of microbes [63,64]. First, acetate was oxidized to H 2 and CO 2 by syntrophic acetateoxidizing bacteria, and was then further metabolized to methane by hydrogenotrophic methanogens. Therefore, competition not only occurred among different trophic microorganisms for population niches, but could also exist between species on the common substrate. Nevertheless, coexistence among various microorganisms with similar trophic niches had an important in uence in maintaining the functional stability and resilience of microbial ecosystems [65,66]. In this study, the stable coexistence of methanogenic core members in different seasons and in different vegetation types also showed that both positive and negative interspecies interactions among members supports the stability of these functional communities.
The characteristics of methanogenic community structure and metabolism according to different vegetation types and seasons may be further represented by methanogenic indicator species. In the current study, methanogenic indicator OTUs can be largely divided into four types of methanogenic metabolic pattern: hydrogenotrophic, acetoclastic, facultative, and methylotrophic. Nine indicator OTUs represented the B. platyphylla -L. gmelinii wetland in summer, and belonged to three methanogenic metabolic patterns. Unclassi ed Thermoplasmata, as an indicator OTU, was detected only in the B. platyphylla -L. gmelinii wetland in summer, which induced more diverse methanogenic metabolic patterns. A. sibirica is a vegetation type in the transition zone from forested wetland to shrub wetland. In the current study, the soil methanogenic community composition of A. sibirica was similar to that of B. ovalifolia, which represents shrub wetland. All indicator OTUs of A. sibirica and B. ovalifolia for the summer were classi ed to Methanosaetaceae. The number of indicator OTUs in these two wetlands was higher in winter than in summer and mainly belonged to Methanobacteriaceae. Some researchers also consider that methanogenesis is implemented by catabolic interactions among microorganisms at different trophic levels [25,67]. Thus, each vegetation type has manifold indicator OTUs which represent the characteristic methanogenic metabolic patterns in different seasons. This paper reveals a new view of the transformation of soil methanogenic metabolic patterns across the transition zone of high latitude forested wetland.

Conclusion
In summary, we detected that soil methanogenic communities in high latitude natural wetlands in China have high sensitivity to season and vegetation type. Compared with season, vegetation type displayed a stronger in uence on the establishment of a methanogenic community. Changes in methanogenic communities further affected the shifts in methanogenic metabolic patterns.
Our results also re ected that the co-existence of multiple metabolic patterns and a steady core community enhanced the resilience of methanogenesis processes to environmental changes. The construction of a network was also useful for assessing symbiotic or competitive interspecies interactions in wetland environments. Thus, we suggest that vegetation type should be taken into consideration when selecting appropriate sites for the restoration of forested wetland in order to optimize environmental bene ts and reduce methane production.

Study site and sampling
The study area is located in Wuyiling Wetland Nature Reserve of Xiaoxing'an Mountain in Northeastern China (E129°00′-129°30′, N48°33′-48°50′). The terrain of this area is characterized by modestly undulating, low hills with an altitude of 350-550 m. The zonal vegetation type is conifer-broadleaf forest dominated by Pinus koraiensis. Its climate belongs to the temperate continental monsoon climate, with an average annual precipitation of 584 mm. The annual average temperature is -1.1 °C, and annual accumulated temperature (≥10 °C) is 1700-2000 °C. The frost-free period is approximately 97 days, and the freezing period is 180 days. The snow accumulation period is approximately 160 days and the average snow depth is 27 cm. The zonal soil type is dark brown soil, while azonal soils include meadow soil, swamp soil, and peat soil [68].
One study site was set up in each of the following wetland types: B. platyphylla -L. gmelinii (a forested wetland, abbreviated as BLW), A. sibirica (a forested wetland; ASW), and B. ovalifolia (a shrub wetland; BOW

Illumina sequencing and data analysis
The raw sequencing data of the methanogenic mcrA gene amplicons were optimized, ltered, and assessed using Trimmomatic [70]. Forward and reverse reads were merged using Flash (version 1.2.11). Sequences with an average quality score less than 20 were ltered, while sequences that contained ambiguous (N) bases were screened out [71]. Exact barcode matching was employed, with a two-nucleotide mismatch permitted for primer matching. The trimmed sequences with chimeras were removed using the Uchime algorithm [72]. To equalize the read sizes of the samples, each sample was downsized to 48 929 reads by random selection for further standardization analysis. The sequences were clustered into operational taxonomic units (OTUs) by setting a sequence similarity standard of 97 %. A representative sequence for each OTU was assigned to analyze taxonomic information based on the functional gene database (http://fungene.cme.msu.edu/). α-diversity indices including Shannon index, Simpson, Chao estimator, ACE, and the Good's coverage estimator was calculated using Mothur at 97% identity.

Statistical analysis
One-way analysis of variance (ANOVA) was used to test the signi cance among the treatments on these alpha-diversity indices in SPSS 16.0 (SPSS Inc., Chicago, IL, United States). A Tukey's multiple comparison test was then applied to determine the signi cance of the differences at a p-value of 0.05. A signi cance test of relative abundance at class and family level was conduceted using STAMP Software [73]. A Venn diagram was used to show the common and unique OTUs using R software (Venn diagram package).
A principal coordinate analysis (PCoA) of OTUs of all samples was generated based on the Bray-Curtis dissimilarity distance test. The soil methanogenic feature of the three vegetation types was determined in order to identify the signi cant biomarkers using linear discriminate analysis effect size (LEfSe) analysis (LDA scores = 2) [74].
A co-occurrence network analysis was built using the relative abundance of top 60 OTUs with a Spearman's rank correlation (coe cient r ≥0.7 and p < 0.05). The network was visualized using Cytoscape (version 3.7.1) [75]. The hierarchical cluster of all samples was constructed alongside the Bray-Curtis dissimilarity distance matrix at the OTU level using the R software's Vegan package. The permutational multivariate analysis of variance (PERMANOVA; 999 permutations) further assessed the signi cant differences within the hierarchical cluster of methanogenic communities.
Indicator species analysis was conducted using the "labdsv" and "vegan" packages in R software. For this analysis, all OTUs were selected. Species with a high indicator value (IndVal; established using a Monte Carlo permutation with 1000 replicates) were considered the most pivotal indicators within the community [76]. Finally, a phylogenetic tree was constructed using the relative abundance of the top 500 species. This was calculated based on the maximum likehood method and computed using FastTree software (version 2.1.3).

Figure 2
Network analysis revealing the connection pattern of the soil methanogenic community according to wetland type and season. The top 60 OTUs were selected according to mean abundance from all samples. The size of each node is proportional to the relative abundance of methanogenic OTUs. Nodes from the same family in the network are marked with identical colors. The thickness of each edge with r ≥ 0.7 and p < 0.05 is proportional to the coe cient value (r).

Figure 3
Principal coordinate analysis of the methanogenic communities found in three different wetland types. Circles represent winter, while squares represent summer. Red represents B. platyphylla -L. gmelinii wetland, blue represents A. sibirica wetland, green represents B. ovalifolia wetland.

Figure 4
Taxonomic cladogram of LefSe analysis from phylum to genus level of the three different vegetation types.

Figure 5
Phylogenetic a liations of methanogenic species from different wetland soils.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. Supplementary.pdf