extfieldsofstudy: Medicine Biology

text:
Baseline microbiota composition modulates antibiotic-mediated effects on the gut microbiota and host

Background Normal mammalian development and homeostasis are dependent upon the gut microbiota. Antibiotics, essential for the treatment and prophylaxis of bacterial infections, can have collateral effects on the gut microbiota composition, which can in turn have far-reaching and potentially deleterious consequences for the host. However, the magnitude and duration of such collateral effects appear to vary between individuals. Furthermore, the degree to which such perturbations affect the host response is currently unclear. We aimed to test the hypothesis that different human microbiomes have different responses to a commonly prescribed antibiotic and that these differences may impact the host response. Methods Germ-free mice (n = 30) humanized with the microbiota of two unrelated donors (A and B) were subjected to a 7-day antibiotic challenge with amoxicillin-clavulanate (“co-amoxiclav”). Microbiome and colonic transcriptome analysis was performed, pre (day 0) and post antibiotics (day 8) and subsequently into recovery (days 11 and 18). Results Unique community profiles were evident depending upon the donor, with donor A recipient mice being dominated by Prevotella and Faecalibacterium and donor B recipient mice dominated by Bacteroides and Parabacteroides. Donor A mice underwent a marked destabilization of their microbiota following antibiotic treatment, while donor B mice maintained a more stable profile. Dramatic and overlapping alterations in the host transcriptome were apparent following antibiotic challenge in both groups. Despite this overlap, donor A mice experienced a more significant alteration in gene expression and uniquely showed correlations between host pathways and key microbial genera. Conclusions Germ-free mice humanized by different donor microbiotas maintain distinct microbiome profiles, which respond in distinct ways to antibiotic challenge and evince host responses that parallel microbiome disequilibrium. These results suggest that inter-individual variation in the gut microbiota may contribute to personalized host responses following microbiota perturbation. Electronic supplementary material The online version of this article (10.1186/s40168-019-0725-3) contains supplementary material, which is available to authorized users.

Background
The mammalian microbiome incorporates a vast and diverse range of microorganisms that have co-evolved to live with us, contributing essential functions through their collective genetic and metabolic repertoire [1]. Ecological diversity and inter-individual variation are hallmarks of the healthy gut microbiome [2] and appropriate assembly of microbial communities in infancy and their maintenance in adulthood is of critical importance for metabolic [3] and immune [4] maturation, with the microbiome playing a role in both immune education and homeostasis. Many common environmental exposures can have demonstrable effects on the microbiome, including mode of birth delivery [5], antibiotic use [6,7], reduced dietary diversity [8], and prescription medications [9], and some associations with these exposures and subsequent development of immune, inflammatory, or metabolic conditions have been made [7,8,[10][11][12][13].
Antibiotic treatment, one of the fundamental medical achievements of the past century, represents a notable example of microbiome perturbation, with rapid and sometimes enduring changes to the community structure [14]. Antibiotic administration during the neonatal "window period"-when the gut microbiota is in part responsible for immune education [15]-may have long-term consequences, and early life antibiotics have been associated with increased risk of asthma, obesity, and Crohn's disease [6,7,11,16]. Furthermore, antibiotic use is directly associated with the development of vancomycin-resistant Enterococcus (VRE) [17], recurrent Clostridium difficile infection (CDI) [18], and the emergence of antibiotic resistance as a global health threat [19]. Notably, response and recovery of human gut microbiota communities following antibiotic treatment can be highly variable between individuals, as described in the context of responses to the antibiotic ciprofloxacin [20] and the routine prescription of antibiotics clarithromycin and metronidazole for the treatment of Helicobacter pylori infection [21]. Whether such changes are restricted to the microbiota or result in a downstream alteration in intestinal host response as well is not currently clear.
Due to the importance of the microbiota for host homeostasis and the marked collateral effects of antibiotics on the microbiota, we posed the questions: (a) can antibiotic administration result in correlated fluctuations of the microbiota and host colonic transcriptome? and (b) if so, is there evidence of a personalized response? To answer these questions, we characterized the microbiota and the colonic transcriptome of germ-free mice (GFM) humanized with fecal samples from two unrelated healthy human donors, before and after treatment with the broadspectrum antibiotic amoxicillin-clavulanate (i.e., co-amoxiclav), commonly used to treat various infections in clinical practice (Fig. 1a).

Results
In total, samples from 30 mice passed quality controls, 15 in each group. One animal from the donor B group at day 0 (D0) was excluded due to low microarray signal intensity values. Four hundred fifty-four OTUs (operational taxonomic units-see the "Methods" section for the definition as we used the dada2 pipeline [22]), were identified following filtering and were assigned to 51 genera (mean sequences per sample following filtering 6089, range 4219-8835, rarefaction curves presented in Additional file 1: Figure S1A).
Baseline microbiome composition modulates susceptibility and recovery following antibiotics exposure Donor A recipient mice (donor A mice) had a markedly different community compared to mice from donor B. Donor A mice had a microbiota composition dominated by Prevotella (mean 26.6% (SD 13.6%)) and Faecalibacterium (mean 24.9% (SD 12.5%)), with Bacteroides (mean 17.3% (SD 7.8%)) being the third most dominant genus (Additional file 2). In donor B mice, the dominant genera were Bacteroides (mean 37.9% (SD 12.2%)) and Parabacteroides (mean 19.6% (SD 7.6%)), with less than 1% mean abundance of both Faecalibacterium and Prevotella (Fig.  1b, c). Samples for each donor were aggregated at each time point and submitted to enterotype clustering with the MetaHIT data set [23,24] (following manual filtering to assure all genera matched) (Additional file 1: Figure  S1B). This confirmed clustering of donor A mice with the Prevotella enterotype and donor B mice with the Bacteroides enterotype. Online reference-based enterotyping [25] similarly classified these samples to the Prevotella and Bacteroides enterotypes.
In contrast to donor A mice, only one genus, Ruminococcus, increased significantly in donor B mice at D18 (Fig. 2b). Overall, there was a much more marked disruption of the community structure in mice from donor A mice when compared to donor B mice, evident in PCoA plots of weighted unifrac distance (Fig. 2c). To determine if there was a significant difference in terms of the effect of donor and day post antibiotics, we performed a two-way PERMANOVA, which indicated a significant interaction between donor and day ( Table 1). A post hoc analysis demonstrated differences only for donor A, at   Table 1. d PCoA of weighted unifrac distances for donor A alone, with genera that correlate significantly with PCoA axes (Spearman's correlation, P value after correction 0.1) and PERMANOVA R 2 and P value for univariate comparison. e PCoA for donor B as in d with weighted unifrac distances. No significant correlations were present. PERMANOVA R 2 and P value for univariate comparison  (Fig. 2d). No correlations were evident for donor B mice (Fig. 2e). While we did not have a noantibiotic control group, we designed the study such that antibiotic-induced changes to the microbiota must occur independently in two cages. To ensure this was the case, we used a distance-based method to determine if withincage differences over time (Additional file 4: Figure S2A) were significantly different from time point-specific distances, including between-cage effects at D8 (post antibiotics) (Additional file 4: Figure S2B and S2C). We observed a significant shrinking in the distance post antibiotics in donor A mice but did not observe a similar reduction in donor B mice, suggesting that the effects of antibiotics in this group were not above the baseline within-cage variability over the course of the study (Additional file 4: Figure S2C). Taken together, these results demonstrated a marked effect of co-amoxiclav administration on the microbiota of donor A mice, with notable shifts in the dominant genus Prevotella and a number of less abundant genera. In contrast, the microbiota profile of donor B mice remained relatively unaffected.
Antibiotics induce dramatic changes on host transcriptome, which overlap independently of baseline microbiota composition Gene expression was assessed initially with PCA for all mice combined (Fig. 3a). This demonstrated a pattern associated with antibiotic treatment, with the main Fig. 3 a PCA plots of gene expression data following antibiotic treatments for all mice combined, clustered by day with respect to antibiotic treatment. b Heatmap of differentially expressed genes at the D18 vs D0 contrast. c GO enrichment analysis, followed by submission to REVIGO, stratified by contrast separation between early (D0 and D8) and recovery (D11 and D18) time points. Figure 3b represents differentially expressed genes between D18 and D0 (1845 in total), which again demonstrate clear clustering between early (D0 and D8) and recovery (D11 and D18) time points. Gene Ontology (GO) analysis with the goana function in limma identified a large number of significant GO terms between different contrasts, so the top 100 were submitted to REVIGO and of the resulting output, those with a − log10(P value) > 5 were plotted. This analysis demonstrated significant changes in many terms, prominently biological processes concerned with extracellular matrix organization, cellular adhesion, signaling, cell migration, developmental processes, and circadian rhythm (Fig. 3c).
A combined analysis of both groups in our study also demonstrated significant differential gene expression alterations at D8 relating to genes associated with circadian rhythm, consistent with previously reported findings following antibiotic treatment, including altered expression of Per genes, Cipc, Prkaa2, NR1D2/Rev-ErbA, Ciart, NPAS2/CLOCK, NFIL3/E4BP4, and ARNTL/ BMAL1 [26] (Additional file 5: Figure S3). Recent work has precisely explored the relationship between circadian rhythm and the microbiome and identified diurnal biogeographic changes in microbial proximity to the epithelium which orchestrate epithelial and host circadian processes [27][28][29].
When the same analysis was repeated for each donor group individually, more significantly differentially expressed genes were detected in the donor A than donor B group (Fig. 4a-c). While only a limited number of differentially expressed genes were evident for donor B mice, there was notable overlap in those genes that were with donor A mice (Fig. 4c). Additionally, heatmaps of contrasts D11 vs D8 (Fig. 4d) and D18 vs D8 (Fig. 4e) including these genes, while dominated by donor A, still demonstrated clear clustering based on early (D0, D8) and late (D11, D18) time points, again suggesting that even though these genes did not meet the threshold for differential expression for donor B, there was a common pattern in response to co-amoxiclav. GO terms with at least four genes were submitted to REVIGO (donor A, Fig. 4f; donor B, Fig. 4g) demonstrated common activation of circadian rhythm terms at D8 compared to recovery time points (D11 and D18), while terms relating to immune system processes, cell adhesion and response to stimulus were evident in donor A mice at D18 when compared to the early time points (D0 and/or D8).
Taken together, these results demonstrate a marked effect of co-amoxiclav administration on the colonic transcriptome of humanized mice. Gene expression analysis identified changes in genes involved in extracellular matrix organization, cellular adhesion and motility, developmental processes, circadian rhythm, and immune system processes. While the pattern of gene expression was similar regardless of the baseline (D0) microbiota, the degree of change appeared to be more marked in mice humanized by donor A than by donor B. As these were the mice that had a more marked alteration in the microbiota composition in response to co-amoxiclav, we next investigated whether correlation between the transcriptome and microbiota was evident.
Host transcriptome-microbiota correlation reveals marked covariance between the microbiota and GO pathways in donor A mice but not in donor B mice Due to the large number of differentially expressed genes across the different contrasts, we used Gene Set Variation Analysis (GSVA) to identify gene pathways associated with antibiotic treatment. To perform correlation with microbial taxa, we used Hierarchical All-against-All significance testing [30], a recently developed technique for determining relationships between multi-omics datasets (Fig. 5, Additional files 6 and 7). This method identifies significant correlations between individual features of both datasets, as well as identifying correlations between clusters of features. Performed using Spearman's correlation values on the genus-level abundance from each donor and the corresponding pathway abundance from GSVA, a large number of correlations were evident for donor A mice following multiple-hypothesis correction, including a number of dominant and/or differentially abundant bacterial taxa from that donor (Prevotella, Dorea, Faecalibacterium, Blautia, and Roseburia). Figure 5 represents a network plot of significant correlations between individual pathways and bacterial genera for donor A mice. The pathways are colored by high-level GO assignments. Furthermore, pathway-bacteria correlations, including additional clusters determined by HAllA, are represented in Additional file 6 (data) and Additional file 7: Figure S4. In contrast, no significant correlations between pathway expression and the microbiota were evident for donor B mice.
Taken together, these results demonstrate correlation between host pathways and microbiome composition was unique to donor A mice. However, as the antibiotic-induced changes in donor B mice were considerably weaker, such correlations are harder to detect and higher numbers of mice would probably be required in this case.

Discussion
We report that GFM humanized with different human donor microbiota maintain a unique community profile and experience individualized responses to oral antibiotic administration. Mice humanized with donor A microbiota were dominated by Prevotella and Faecalibacterium, with a lesser contribution from Bacteroidetes and underwent marked fluctuations in terms of their dominant genera, most prominently Prevotella and Faecalibacterium, as well as in other lower abundance genera within the community. It is notable that by D18, the mean abundance of Prevotella was approximately half that of baseline (D0), while the mean proportional abundance of Faecalibacterium was almost three times the baseline (D0) (Fig. 1b and Additional file 2), suggesting that the donor A group had a lower resilience, without return to baseline by the end of the study.
In contrast to donor A, donor B mice were much more resilient to antibiotic administration, with the notable exception of a late bloom in Ruminococcus sequences at D18. Bacteroides was the dominant genus in this group with Parabacteroides, the second most dominant genus. As mentioned, there was a trend of decreasing proportional abundance in Bacteroides throughout the study in the donor B group, although this did not reach significance. Notably, most dramatic changes in microbiota composition occurred in the recovery period.
We did not observe a significant change in α-diversity over the duration of the study. We note that there were only three baseline (D0) mice in each group and this may have accounted for our inability to detect a significant change. Furthermore, as we did not detect changes in general in the donor B group, it is consistent that we did not detect a change in α-diversity either. In terms of donor A mice, we did note a drop-off in the Simpson diversity at D18 (Fig. 1d) but changes in this donor group appeared to be dominated by compositional changes rather than a reduction in richness (unique OTUs).  Fig. 3a, isolating donor A gene expression. b PCA for donor B expression alone. c Contrasts with differential genes and the relative contribution from donor A and donor B. Red represents donor A, blue donor B, and white overlap. d D11 vs D8 heatmap of differential gene expression. e D18 vs D8 heatmap of differential gene expression. f GO enrichment analysis following submission to REVIGO, stratified by contrast for donor A. g As in f but for donor B Marked effects of antibiotics on the colonic transcriptome were also demonstrated, again most notably in the recovery period. In a pooled analysis, there was a significant increase in many processes related to biological adhesion and sensing (Fig. 3), perhaps due to stimulation from a rebound recovery in bacteria following the withdrawal of antibiotics. When analyzed individually, the gene expression changes were more marked for donor A mice, suggesting that the more dramatic perturbations in composition associated with this donor group were reflected in the transcriptome. Stronger evidence for this conclusion can be drawn from the multiple correlations detected between enriched GO pathways and the microbiota composition of donor A mice (Fig. 5), in contrast to those of donor B. While these correlations only represent an association, it must also be noted that the total bacterial load was not quantified in this study and as such, the changes in composition may reflect significant changes in mucosal-associated bacteria, which were perhaps more dramatically altered in donor A mice.
Recently, a study has elegantly elucidated the effects of an ablative cocktail of antibiotics on the colonic transcriptome, identifying that depletion of the microbiota, direct antibiotic effects on the host, and the effects of antibiotic-resistant microbiota account for the alterations [31]. This study also demonstrated downregulation of immunity with antibiotic treatment, while identifying significant direct toxicity of antibiotics on mitochondrial function. Interestingly, a recent study in healthy adults using a cocktail of antibiotics demonstrated dramatic, rapid shifts in microbiome composition, followed by near-complete recovery within 6 months [32]. In contrast, our study utilized a more gentle perturbation with a single antibiotic combination commonly used in human medicine, in order to determine correlations between disturbances in the microbiota and the transcriptome, in addition to their relationships to baseline microbiota composition. Our data suggest that while the response to antibiotics is very different between groups of GFM humanized by different donors, and therefore enterotypes, gene expression patterns in response to antibiotics are similar. However, the magnitude of the change detected was more evident in the donor A group, suggesting that the magnitude of microbiota perturbation influences the subsequent changes in gene expression. Such a conclusion fits with numerous observations in murine models highlighting the importance of the baseline microbiota, its impact in a range of settings, and the requirement for methods to quantify and standardize microbiome effects in animal models [33]. These data also suggest that in multiomics studies incorporating the microbiome, temporal analysis within individuals may yield more insights than aggregating data between individuals with very different baseline microbiota compositions. This was clearly evident in our results, where aggregation of data led to a complete loss of the covariance structure due to the different levels and responses of shared bacterial genera between donor groups. These findings would also fit more broadly with emerging data indicating that our microbiome represents an important element of our personalized response to environmental factors [34,35].
A number of limitations should be noted in our study. Only two donors groups were available for analysis and therefore findings relating to the behavior of individual genera, both in their response to antibiotics and their association with host transcriptome response, cannot necessarily be extrapolated to include a range of donor microbiome compositions. Furthermore, varying doses of antibiotics, administration protocols, and antibiotic cocktails have been used in previous studies, limiting direct comparisons. It is important to note that we did not sequence the original donor communities and therefore we cannot characterize precisely the true extent of engraftment of strains or the resulting compositional similarity to the original fecal sample. We have based our conclusions on this model protocol as it is established in our facility in multiple studies that replicate phenotypic traits of the donor [36][37][38][39], as well as from other groups, where seminal findings have established that up to 85% of human taxa can be established following oral gavage [40]. However this study and other similar studies have used higher initial concentrations, and other studies have used higher quantities again [41], as well as repeated dosing schedules [42]. We do however note that the microbiota communities in our study clustered closely with human donors from the MetaHIT consortia, suggesting we had recapitulated an adequate proportion of the donor microbiota in each group. In addition to these specific points in relation to our study, humanized GFM also suffer from a number of limitations, including the effects of a non-native microbiota and impaired immune development in the humanized mouse intestine and these must also be acknowledged [43].
In addition, we did not use a no-antibiotic control group but instead housed mice in one of two separate cages for each donor group. This design implicitly required that changes be consistent across cages in each donor group. We acknowledge however that this represents a limitation in our study. We have used the cage effects in the study to demonstrate significant convergence of microbiota communities in donor A post antibiotics (D8) and not in donor B. However, conclusions from both the microbiota and host transcriptome must take into account the lack of a no-antibiotic control group.
Finally, recent work has demonstrated the importance of bacterial load in characterizing the microbiome, and changes in total bacterial counts in both donor groups may have contributed to some of the effects seen [44].

Conclusions
In summary, we have demonstrated that a commonly prescribed antibiotic can exert significant effects on the gut microbiota, in concert with colonic gene expression. Notably, while the host response appeared to be conserved, the degree of activation was related to the degree of microbiota perturbation, resulting in correlations between host transcriptome and the microbiota unique to one donor community. These findings offer a window into the role that our microbiome plays in mediating our personalized response to medications such as antibiotics.

Study design
All procedures in relation to the care and use of laboratory animals were carried out in accordance with European guidelines and the local ethics committee approved this study. Germ-free male 6-week-old mice (n = 31) were transferred into two sterile separate isolators (Table 2; ANAXEM, INRA, Jouy-en-Josas, France) for the 10-week duration of the study and had ad libitum access to irradiated feeds and sterile water. Mice were observed once a day to ensure their welfare. For colonization, GFM were inoculated orally with 250 μl of a 10 −2 dilution of whole fecal homogenate from two healthy adults (termed donors A and B) with no history of intestinal problems or intake of antibiotics 3 months prior to the beginning of the study, resulting in a single donor fecal microbiota transplantation per isolator (donor A (n = 15) and donor B (n = 16)). Mice in each isolator were housed in two separate cages (Table 2), with the aim to have four mice per time point (two per cage) following antibiotics in each group. The human donors who provided a fecal sample gave their verbal consent to do so. Fifty-six days after microbiota inoculation, samples were acquired at baseline (day 0 (D0); donor A (n = 3), donor B (n = 4)), followed by the administration of the penicillin antibiotic amoxicillin in combination with the β-lactamase inhibitor clavulanic acid once a day (80 mg and 10 mg/kg) to all remaining mice. We chose this antibiotic as it is a commonly prescribed, broad-spectrum human antibiotic, and treatment was continued for a total of 7 days, similar to the duration of a typical treatment for a community infection [45], followed by sampling at days 8 (D8), 11 (D11), and 18 (D18) (Fig.  1a). We choose this dose as it is intermediate between doses used in previous rodent studies [46,47] and is the current recommended dosing for children under 40 kg body weight in France [48].

Sampling
At each time-point, a subset of animals was sacrificed (see Table 1), the large intestine (colon) was removed and luminal contents acquired and processed for the extraction of community microbial DNA. RNA extraction proceeded directly from the large intestine. We did not have prior data to estimate the duration of the recovery period to monitor and so early recovery (D11, 3 days post-cessation of antibiotics) and late recovery (D18, 10 days post-cessation of antibiotics) time points were chosen.

Gene expression
Total RNA was isolated from the large intestine with the RNeasy Mini Kit (Qiagen, Courtaboeuf, France) according to the manufacturer's instructions. RNA 6000 Nano chips on a Bioanalyzer 2100 (Agilent Technologies, Les Ulis, France) were used to assess RNA integrity. Gene expression profiling was performed with the SurePrint G3 Mouse GE 8x60 K Microarray (Design ID: 028005, Agilent Technologies). Further details are provided in supplementary methods in Additional file 8.

DNA extraction and sequencing
Fecal samples were weighed and resuspended in 250 μl of 4 M guanidine thiocyanate in 0.1 M Tris at pH 7.5 (Sigma) and 40 μl of 10% N-lauroyl sarcosine (Sigma) for 10 min at room temperature. Five hundred microliters of 5% N-lauroyl sarcosine in 0.1 M phosphate buffer (pH 8.0) was then added and incubated at 70°C for 1 h. A mixture of sterilized silica beads (0.1-and 0.6-mm diameter) was added and the tube shaken three times at 6.5 m/s for 30 s in a FastPrep apparatus (MP Biomedicals). Fifteen milligrams of polyvinylpolypyrrolidone was then added to the tube, and this was vortexed and centrifuged at 20,000×g for 5 min. Following the recovery of the supernatant, the pellets were washed with 500 μl of TENP (50 mM Tris (pH 8), 20 mM EDTA (pH 8), 100 mM NaCl, 1% polyvinylpolypyrrolidone) and centrifuged for 5 min at 20,000×g and the new supernatant was added to the first supernatant. The washing step was repeated twice. The pooled supernatant was centrifuged briefly to remove particles and then split into two 2-ml tubes. Precipitation of nucleic acids was performed by adding one volume of isopropanol for 10 min at room temperature with centrifugation at 20,000×g for 10 min. Pellets were then resuspended and pooled in 450 μl of 100 mM phosphate buffer at pH 8 and 50 ml of 5-M potassium acetate. The tube was placed on ice overnight and then centrifuged for 30 min at 20,000×g. The supernatant was subsequently transferred to a new tube containing 20 μl of RNase (1 mg/ml) and incubated at 37°C for 30 min. Nucleic acid precipitation was performed by adding 50 μl of 3-M sodium acetate and 1 ml of absolute ethanol. The tube was incubated for 10 min at room temperature and the nucleic acids were recovered by centrifugation at 20,000×g for 15 min. The resulting DNA pellet was finally washed with 70% ethanol and then dried and resuspended in 100 μl of Tris-EDTA (TE) buffer.

Bioinformatic and statistical analysis Transcriptome analysis
Following the pre-processing of raw transcriptome data (described in detail in Additional file 8), statistical analysis was performed in the R statistical environment (R version 3.6.0) [49].

Differential gene expression analysis
Principal component analysis was performed using the ade4 (v1.7-13) [50] and factoextra (v1.0.5) packages in R. Differential gene expression was performed using the limma package (v3. 36.5) and the eBayes test [51]. Multiple testing corrections were made using the Benjamini-Hochberg method (adjusted P values < 0.05). Differentially enriched pathways in the Gene Ontology (GO) databases were identified using the goana function in the limma package. To remove redundancy, GO terms that were significantly enriched (adjusted P value < 0.01) were submitted to REVIGO [52] and the output presented by way of -log10 (P value).

Microbiome analysis
Demultiplexing of pyrosequencing fasta and quality files was performed with QIIME (v1.9.1) [53] and FASTQ files were created which were loaded into R and analyzed using the DADA2 package (v1.10.0), optimized for pyrosequencing output [54] (see Additional file 8) to give a table of amplicon sequence variants (referred to as operational taxonomic units (OTUs) throughout). Taxonomic assignments were performed using the RDP v16 (Ribosomal Database Project) training set specially formatted for the DADA2 package [55]. Phylogenetic tree construction utilized the DECIPHER [56] and Phangorn [57] packages. All subsequent microbiome analyses were performed through the phyloseq (v1.24.2) and vegan (v2.5-5) packages in R [22,58,59]. Alpha diversity was presented using the Simpson diversity index. Differential abundance between donor groups and between different time points within donor groups was assessed using DESeq2 (v1.24.0) (FDR P values < 0.01) and significance was confirmed by Mann-Whitney testing [60]. Beta diversity was assessed with principle coordinate analysis (PCoA) of weighted unifrac distances [61] on proportional abundance data (proportion normalized). A two-way PERMANOVA was performed, using donor and day as main effects and the interaction between donor and day as the interaction term. Post hoc pairwise testing was performed using the "pairwise.perm.manova" function from the RVAideMemoire package (v 0.9-73) [62]. Spearman's correlation of genus-level microbial taxa was also performed against the first two principal coordinate axes for each donor group individually to illustrate bacterial genera associated with major variation. P values were adjusted with the Benjamini-Hochberg method and correlations with adjusted P values < 0.1 were retained. Graphs were made using ggplot2 (v3.1.1) [63] (Additional file 9).

Correlation of microbiota and transcriptome
To determine if there was evidence of the correlation between antibiotic-induced changes in the microbiota and the transcriptome, we firstly used Gene Set Variation Analysis (GSVA) to identify enriched GO pathways across the entire transcriptome data set, using the R package GSVA (v1.32.0) [64]. GO pathways were derived from a murine version of the Molecular Signatures Database (MSigDB) [65]. GSVA output was then submitted to limma to determine differentially expressed pathways. This allowed for dimensionality reduction, reducing the total number of transcripts (18,217) and differentially expressed genes (3630) to 182 differentially expressed pathways across all contrasts. Due to the extreme differences in microbiota composition between the two donor groups, pathway scores were correlated with genus-level bacterial abundance for each group individually. Correlation was performed using the Hierarchical All-against-All (HAllA) approach developed for multiomics data sets [30], having discarded pathways with log fold changes of < 0.5. Benjamini-Hochberg-adjusted P values of < 0.1 were retained from the results and were plotted using ggraph and igraph [66] (Additional file 10). Availability of data and materials Microarray data is available from the Gene Expression Omnibus (GEO) with accession number GSE131785, while microbiota sequencing data is available from the Sequence Read Archive (SRA) with accession number PRJNA544904. Files detailing the code used for the analysis are available as R knitr.html files. For the microbiota analysis, this is in Additional file 9, and for the transcriptomic analysis, this is in Additional file 10.
Ethics approval and consent to participate All procedures in relation to the care and use of laboratory animals were carried out in accordance with European guidelines, and the local ethics committee approved this study.

Consent for publication
Not applicable.


extfieldsofstudy: Chemistry Medicine

s2fieldsofstudy: Biology Chemistry

text:
Is MYND Domain-Mediated Assembly of SMYD3 Complexes Involved in Calcium Dependent Signaling?

Macromolecular complexes are essential to intracellular signal transduction by creating signaling niches and enabling a chain of reactions that transmit external signals into various cellular responses. Analysis of SMYD3 interactome indicates this protein lysine methyltransferase might be involved in calcium dependent signaling pathways through forming complexes with the phospholipase PLCB3, calcium/calmodulin dependent kinase CAMK2B, or calcineurin inhibitor RCAN3. SMYD3 is well-known as a histone H3K4 methyltransferase involved in epigenetic transcriptional regulation; however, any roles SMYD3 may play in signaling transduction remain unknown. KEGG pathway enrichment analysis reveals the SMYD3 interacting proteins are overrepresented in several signaling pathways such as estrogen signaling pathway, NOD-like receptor signaling pathway, and WNT signaling pathway. Sequence motif scanning reveals a significant enrichment of PXLXP motif in SMYD3 interacting proteins. The MYND domain of SMYD3 is known to bind to the PXLXP motif. The enrichment of the PXLXP motif suggests that the MYND domain is likely to be a key interaction module that mediates formation of some SMYD3 complexes. The presence of the PXLXP motifs in PLCB3 and CAMK2B indicates the potential role of the MYND domain in mediating complex formation in signaling. The structural basis of SMYD3 MYND domain-mediated interactions is unknown. The only available MYND-peptide complex structure suggests the MYND domain-mediated interaction is likely transient and dynamic. The transient nature will make this domain well-suited to mediate signaling transduction processes where it may allow rapid responses to cellular perturbations and changes in environment.

INTRODUCTION SMYD3 belongs to a special class of protein lysine methyltransferases that is classified by the presence of SET and MYND domains (Spellmon et al., 2015). SET domain is an evolutionarily conserved domain responsible for protein lysine methylation (Tian et al., 2013). MYND domain is a zinc finger motif involved in protein-protein interaction with transcriptional co-repressor complexes (Liu et al., 2007). SMYD3 is a ubiquitously expressed protein with the highest expression in brain and thyroid tissues (Uhlen et al., 2015). In cells, SMYD3 is found in both the nucleus and cytoplasm (Thul et al., 2017). Current research of SMYD3 focuses primarily on its roles in tumor cell growth and muscle development. SMYD3 expression is abnormally elevated in over 15 types of cancers (Spellmon et al., 2017). Overexpression of SMYD3 is associated with tumor cell growth (Hamamoto et al., 2004). While SMYD3 is required for the development of cardiac and skeletal muscles in zebrafish (Fujii et al., 2011), knockout of SMYD3 in mice shows no significant muscle phenotypes (Bagislar et al., 2016). However, in mice SMYD3 mediates the recruitment of transcriptional cofactors at the myostatin and c-Met genes and regulates skeletal muscle atrophy (Proserpio et al., 2013). Current understanding of the molecular mechanisms by which SMYD3 exerts its functions is largely limited to its roles in epigenetic control of gene expression. SMYD3 can upregulate expression of a myriad of oncogenes including c-MYC, STAT3, and matrix metalloproteinase-9 (Cock-Rada et al., 2012;Liu et al., 2015;Sarris et al., 2016). In almost all cases binding of SMYD3 to the promoter region of the target genes is associated with both H3K4 trimethylation and gene activation (Kim et al., 2009;Cock-Rada et al., 2012;Liu et al., 2012;Luo et al., 2014). However, the in vitro methyltransferase activity of SMYD3 toward H3K4 is much lower compared to other targets such as histone H4K5 (Van Aller et al., 2012). This has raised a question of whether H3K4 is a true in vivo epigenetic target (Van Aller et al., 2012;Mazur et al., 2014). SMYD3 also methylates non-histone proteins. Methylation of MAP3K2, a mitogen-activated protein kinase, links SMYD3 to the MAPK signaling pathway and Ras-driven cancer (Mazur et al., 2014). Methylation of vascular endothelial growth factor receptor-1 (VEGFR1) regulates VEGFR1 kinase activity potentially affecting VEGFR1-dependent angiogenesis (Kunizaki et al., 2007). While these studies have established the methyltransferase activity of SMYD3 being essential for its function, the catalytic-independent roles have been uncovered in several cases. Overexpression of a methylation-defective mutant of SMYD3 transactivated the promoter of MYC gene to a similar degree as that observed with the wild type protein (Sarris et al., 2016). SMYD3-mediated EMT (epithelial mesenchymal transition) gene regulation is also independent of SMYD3 methylation activity (Fenizia et al., 2019). Both the wild type and inactive SMYD3 mutant promoted upregulation of the mesenchymal genes N-cadherin and Vimentin and equally reduced E-cadherin and Occludin protein levels. Although the precise mechanism remains largely unknown, the proteinprotein interaction events have been proposed to contribute to the catalytic-independent function of SMYD3 (Sarris et al., 2016;Fenizia et al., 2019).

SMYD3 Interacts With Proteins Involved in Signaling
To date, a total of 43 human proteins were found to interact with SMYD3 (Stark et al., 2006;Fahey et al., 2011;Lopez et al., 2015;Szklarczyk et al., 2015). Among them 30 proteins were considered as high-confidence interactors based upon a reliability score for interaction (see Methods). We found these high-confidence interactors are enriched in two major functional groups, chromatin proteins that are involved in epigenetic transcriptional regulation and signaling molecules that are related to calcium dependent signaling ( Figure 1A). Interacting with the phospholipase PLCB3 may link SMYD3 to calcium-dependent lipid signaling pathways. PLCB3 belongs to the phospholipase C (PLC) subtype that plays an important role in many signal transduction processes (Cocco et al., 2015). Many extracellular signals including many hormones, neurotransmitters, and growth factors converge their physiological action on PLC enzymes which subsequently propagate the signals leading to many cellular responses. PLC enzymes propagate extracellular signals by catalyzing the hydrolysis of the membrane bound phosphatidylinositol 4,5-bisphosphate (PIP2) into the second messengers inositol 1,4,5-trisphosphate (IP3) and diacylglycerol (DAG). IP3 and DAG then further propagate the signals by inducing the calcium release from endoplasmic reticulum and activating protein kinase C (PKC), respectively. Calcium is an essential downstream effector in PLC signal transduction pathways required for activation of many protein kinases including PKC and calcium/calmodulin dependent protein kinases (CAMKs). Therefore, interacting with CAMK2B further suggests SMYD3 might be involved in calcium-dependent signaling pathways. CAMK2B belongs to a family of serine/threonine specific kinases whose activation is initially dependent on the binding of calcium/calmodulin (Swulius and Waxham, 2008). CAMKs play an integral role in translating intracellular calcium signals to a variety of downstream targets. Upon activation, CAMK2B has been shown to phosphorylate a range of targets involved in learning, memory, synaptic plasticity, and long-term potentiation (Lisman et al., 2012;Zalcman et al., 2018).
Interacting with RCAN3 indicates SMYD3 might modulate additional aspects of calcium-dependent signaling. RCAN3 belongs to a family of small evolutionarily conserved proteins that can directly bind and inhibit calcineurin (Mulero et al., 2007). The inhibition regulates calcineurin-dependent signaling events in which calcineurin is activated by calcium and calmodulin. SMYD3 could thus indirectly affect calcineurin activation and function by binding to RCAN3. The best known function of calcineurin is to control the transcription of genes involved in cardiac hypertrophy development and T cell activation by dephosphorylating and activating nuclear factor of activated T cells (NFAT) (Molkentin, 2004;Hogan, 2017).

Pathway Analysis of SMYD3 Interactome
Gene Ontology analysis reveals the SMYD3 interacting proteins are enriched in positive regulation of intracellular signal transduction, cell surface receptor signaling, cellular response to chemical stimulus, and regulation of MAPK cascade ( Figure 1B). KEGG (Kyoto Encyclopedia of Genes and Genome) pathway enrichment analysis also reveals these proteins are overrepresented in several signaling pathways including estrogen signaling pathway, GnRH signaling pathway, NOD-like receptor signaling pathway, and WNT signaling pathway ( Figure 1C). In estrogen signaling, SMYD3 may interact with the proteins in both nuclear-initiated and membrane-initiated signaling pathways (Figure 2). In the nuclear pathway, SMYD3 may interact with estrogen receptor α (ERα) and the heat shock protein 90 (HSP90). The interaction between SMYD3 and the FIGURE 1 | Pathway enrichment analysis of SMYD3 interactome. (A) SMYD3 interactors are enriched in two major functional groups: chromatin organization and signal transduction. Lines are drawn from SMYD3 to its interactors and between SMYD3 interactors if their interactions are documented in STRING with the interaction score larger than 0.3. A selected list of enriched GO biological processes (B) and KEGG pathways (C). Numbers at the right end of the bars indicate the false discovery rate (FDR) and fold enrichment in parenthesis.
ligand binding domain of ERα is required for recruitment of SMYD3 to the proximal promoter regions of ER target genes (Kim et al., 2009). The recruited SMYD3 is responsible for the accumulation of di-and trimethylation of H3K4 and subsequent ER-regulated gene transcription. The role of its interaction with HSP90 is unknown. HSP90-based chaperone FIGURE 2 | Enriched KEGG signaling pathways. SMYD3 interacting proteins are highlighted in orange. Dashed lines indicate the interactions between SMYD3 and its interactors. Arrow lines indicate signaling pathways with GnRH signaling, NOD-like signaling, ER signaling, and WNT signaling colored in red, blue, green, and yellow, respectively.
complex is thought to repress the ER transcriptional regulatory activities while maintaining the receptor in a conformation that is competent for high-affinity steroid binding (Knoblauch and Garabedian, 1999). However, HSP90 can enhance SMYD3 methyltransferase activity suggesting HSP90 could also positively regulate ER transcriptional activity through SMYD3 (Brown et al., 2015). In the membrane pathway, estrogen can exert its action through a subpopulation of ER at the plasma membrane (mER) or G-protein coupled E2 receptors (GPERs) (Song, 2007). SMYD3 may regulate mER signaling through binding to PLC. The PLC-mediated intracellular calcium release is central to the mechanism of mER signaling (Chaban et al., 2004). The released calcium can affect various downstream signaling pathways including calcium, cAMP and PI3K-AKT pathways and lead to activation of CREB family transcription factors and various ER cellular responses (Manavathi and Kumar, 2006).
Also through binding to the PLC SMYD3 may be involved in GnRH signaling and NOD-like receptor signaling (Figure 2). In GnRH signaling pathway, gonadotropin-releasing hormone (GnRH) acts upon its receptor in the anterior pituitary to regulate the production of gonadotropins (Perrett and McArdle, 2013). PLC is responsible for transducing the GnRH signal, which activates calcium signaling pathway, protein kinase C (PKC) pathway, and MAPK signaling pathway (Call and Wolfe, 1999). SMYD3 may regulate CAMK activity in calcium signaling and methylate MAP3K2 in MAPK signaling. In NOD-like receptor signaling pathway, PLCB-mediated calcium signaling affects NLRP3 inflammasome assembly and activation (Murakami et al., 2012). The activated NLRP3 induces caspase-1 activation, which then regulates maturation of the pro-inflammatory cytokines IL-1B and IL-18 and drives pyroptosis (He et al., 2016). SMYD3 may affect NOD2 signaling also through binding to HSP90. The proper functioning of NLR proteins (NOD2) is dependent on HSP90. It has been shown that NOD2 interacts with HSP90 via both CARD domains; the interaction enhances NOD2 stability (Boyle et al., 2014).
SMYD3 interacts with the proteins that play roles at several points in WNT signaling (Figure 2). WNT signaling is required for basic developmental processes such as cell-fate specification, progenitor-cell proliferation and the control of asymmetric cell division, in many different species and organs (Nayak et al., 2016). There are at least three different WNT pathways: the canonical pathway, the planar cell polarity (PCP) pathway and the WNT-Ca 2+ pathway. SMYD3 interacts with the proteins involved in WNT-Ca 2+ signaling and PCP signaling. WNT-Ca 2+ signaling is mediated through G proteins and phospholipases and leads to transient increases in cytoplasmic free calcium that subsequently activates the kinase PKC and CAMK2 and the phosphatase calcineurin (Kuhl et al., 2000). In this signaling, SMYD3 may interact with PLCB, CAMK2, and the calcineurin inhibitor RCAN3. SMYD3 may regulate PCP signaling through binding to VANGL1. Van Gogh-like (VANGL) proteins are core components of the planar cell polarity pathway that controls epithelial polarity and cell migration (Mentink et al., 2018). PCP signaling leads to the activation of the small GTPases RHOA and RAC1, which activate the stress kinase JNK and ROCK and lead to remodeling of the cytoskeleton and changes in cell adhesion and motility (Yang and Mlodzik, 2015).

Functional Associations Between SMYD3 and Its Interactors
Appreciation of SMYD3 involvement in signaling would broaden our view of the role of SMYD3 in cancer beyond epigenetic regulation. The diversity of cellular signaling exploited by cancer is immense. Oncogenic mutations have been found in numerous signaling pathways including WNT signaling, calcium signaling, and MAPK signaling cascade (Sever and Brugge, 2015). In each case, the balance of signaling is disrupted by the mutations, which is not subject to the normal control mechanisms. Disrupted cell signaling allows cells to overproliferate and causes changes in the tumor microenvironment, angiogenesis, inflammation, and cell metabolism, and consequently promotes cancer progression and metastasis and evades apoptosis and immune destruction. The interaction of SMYD3 with proteins involved in these signaling pathways would provide additional mechanistic links between SMYD3 overexpression and cancers.
The available data also emphasize potential functional associations between SMYD3 and its interactors in normal physiology in addition to cancer. SMYD3 is a ubiquitously expressed protein with the highest expression levels in brain, thyroid, and reproductive organs (Uhlen et al., 2015). However, the precise roles of SMYD3 in these tissues remain completely unknown. The interaction with GnRH signaling could suggest a possible functional integration of SMYD3 in brain and testis. GnRH signaling is the primary regulator of mammalian reproductive function in both males and females controlling gametogenesis and steroidogenesis (Perrett and McArdle, 2013). The interaction with this signaling pathway would suggest that SMYD3 might be involved in the endocrine control of spermatogenesis by regulating the neuroendocrine activity along the hypothalamic-pituitary-testicular axis.
There might be potential functional connections between SMYD3 and CAMK2B in muscle atrophy. SMYD3 regulates skeletal muscle atrophy through the recruitment of transcriptional cofactors to the myostatin and c-Met genes (Proserpio et al., 2013). Reducing SMYD3 decreased myostatin and c-Met transcription, which prevented muscle loss and protected from glucocorticoid-induced myotube atrophy. CAMK2B is also linked to muscle atrophy. CAMK2B was identified as a downstream target of p38α MAPK and positively regulated muscle atrophy (Clarke et al., 2007). The pharmacological inhibition of CAMK2B activity suppressed denervation-induced muscle atrophy. SMYD3 could contribute to p38α MAPK-mediated muscle atrophy through binding to CAMK2B. CAMK2B was found to regulate MuRF1 transcription (Clarke et al., 2007). MuRF1 is an E3 ubiquitin ligase that degrades myosin heavy chain causing the breakdown of skeletal muscle under atrophy conditions (Yuasa et al., 2018). MuRF1 has been shown to be upregulated during both denervation and glucocorticoid treatment (Furlow et al., 2013). It is possible that SMYD3 could play multiple transcriptional regulatory roles in glucocorticoid-induced myotube atrophy. It is not only involved in epigenetic regulation of myostatin and c-Met genes but also could affect CAMK2B-dependent MuRF1 gene expression through binding to CAMK2B.
The interaction with NOD-like receptor signaling or PLCB3 might link SMYD3 to inflammation and the immune response. NOD-like receptors are specific families of pattern recognition receptors responsible for detecting various pathogens and generating innate immune responses (Mogensen, 2009). PLCB3 is required for FceRI-mediated mast cell activation (Xiao et al., 2011). Loss of PLCB3 leads to defective mast cell-dependent immune responses. PLCB3 also regulates the cystic fibrosis (CF) inflammatory response by amplifying the expression and release of IL-8 and recruiting neutrophils to CF airways (Bezzerri et al., 2011). SMYD3 has been shown to play a role in regulatory T cell development. SMYD3 regulates Foxp3 expression via its histone H3K4 methyltransferase activity (Nagata et al., 2015). This epigenetic control is an important mechanism in regulating TGFβ-dependent regulatory T cell formation. SMYD3 abrogation affected regulatory T cell formation while allowing dysregulated interleukin-17 production. The interaction with NOD-like receptor signaling and PLCB3 could suggest that SMYD3 might also be involved in innate immune response in addition to regulation of adaptive immunity.

Structures of SMYD3 and Its Substrate Complex
Like other SMYD proteins (SMYD1 and 2), SMYD3 structure is bilobal ( Figure 3A) (Sirinupong et al., 2010Jiang et al., 2011;Spellmon et al., 2015). The N-terminal lobe contains four domains: SET, MYND, SET-I, and post-SET. The C-terminal lobe contains a tetratricopeptide repeat (TPR)-like domain that is organized into seven up-down helices. There are at least three scenarios in that SMYD3 can form complex with other proteins. The SET domain is required for the formation of enzymesubstrate complexes (Spellmon et al., 2015). The MYND domain is a putative protein-protein interaction module that interacts with a proline-rich sequence (Liu et al., 2007). The TPR-like CTD domain has been shown to interact with HSP90 but the exact nature of such an interaction is unknown (Brown et al., 2015).
The structure of SMYD3 in complex with a MAP3K2 peptide provided the structural basis for SMYD3-substrate complex formation (Fu et al., 2016). The MAP3K2 peptide was found to bind at the bottom of a deep cleft formed between the N-and C-lobes ( Figure 3A). Similar to the SMYD2-substrate complexes (Ferguson et al., 2011;Wang et al., 2011;Jiang et al., 2014), the bound peptide adopts a U shaped conformation. The base of the U shape is sandwiched between a β-hairpin from the SET domain and a loop preceding the post-SET domain. The descending arm of the U shape is bound by the SET-I domain, while the ascending arm interacts with the CTD. The target lysine is located within the U base. The side chain of the target lysine is recognized by a long narrow channel through which it is deposited into the active site. The structural basis of SMYD3-non-substrate complexes remains unknown.

MYND Domain-Mediated Protein Interaction
The MYND is likely an important module that contributes to the assembly of several SMYD3 complexes. MYND domain is present in about 90 human proteins (Hunter et al., 2009). The role of MYND is largely unexplored and limited to protein interaction and recruitment of transcriptional co-repressors (Liu et al., 2007). In the SMYD protein family, SMYD1 is known to interact with skNAC through its MYND domain (Sims et al., 2002). Cysteine to serine mutations in the first or second zinc fingers of the MYND domain abolished the ability of SMYD1 to bind skNAC in immunoprecipitation assays. The interaction of SMYD2 with EPB41L3 is also mediated by the MYND domain (Abu-Farha et al., 2008). The wild type but not the MYND-deleted form of SMYD2 was found to co-immunoprecipitate with EPB41L3. The MYND domain of SMYD3 is also required for its binding to the co-repressor N-CoR (Foreman et al., 2011). N-CoR coimmunoprecipitated with wild type SMYD3 but not with the MYND domain point mutant C49S. MYND domain appears to prefer interacting with a proline-rich sequence. The MYND of BS69 was found to bind to a PXLXP motif (Ansieau and Leutz, 2002). A leucine to alanine mutation in the PXLXP motif abolished the interaction. ETO and DEAF-1 MYND domains were found to interact with PPPLI motif (Liu et al., 2007;Kateb et al., 2013). Disrupting the interaction between the ETO MYND domain and the SMRT PPPLI motif attenuated AML1/ETO effects on proliferation, differentiation, and gene expression (Liu et al., 2007). The MYND domains of SMYD proteins are likely to bind to the PXLXP motif (Sims et al., 2002;Abu-Farha et al., 2008;Foreman et al., 2011). Co-immunoprecipitation experiments performed in transiently transfected 293T cells identified the PPLIP motif in skNAC as an important SMYD1 interaction motif as the substitution of the leucine to alanine failed to associate with SMYD1 (Sims et al., 2002).
MYND domain adopts a conserved β-β-α topology in which two anti-parallel β-strands and one small kinked α-helix organize around two zinc atoms. The zinc atoms are chelated by seven cysteine residues and one histidine in a C4-C2HC arrangement ( Figure 3B). The solution structure of ETO in complex with a PPPLI sequence from the co-repressor SMRT revealed the peptide binding mode and residues responsible for peptide recognition (Liu et al., 2007). The bound peptide adopts a stretched conformation (Figure 3C). It is connected to the MYND domain by forming a β strand that pairs with the β strands of the MYND domain. The residue leucine (L) in the peptide binds at a shallow surface pocket containing the residues tyrosine and glutamine. The first proline (P) is packed against by a tryptophan residue. There are no significant specific interactions between the rest positions of the peptide and MYND domain.
MYND domain-based interactions show relatively low binding affinities with K D in the millimolar range (Kateb et al., 2013). The low binding affinities are consistent with the structural finding that ETO and SMRT bind with much smaller interface areas (∼400 Å 2 ). The weak interaction suggests the MYND-mediated interactions are transient in nature. Transient interactions usually show a fast bound-unbound equilibrium and can readily switch between different binding partners (Perkins et al., 2010). However, in some cases effective MYND interaction requires additional contacts between the other regions of proteins. The N-terminal residues located in the Ssequence of SMYD1 are required for efficient interaction with the PXLXP motif in skNAC (Sims et al., 2002). A set of positively charged residues located at the C-terminus of BS69 MYND domain is crucial for the interaction with PXLXP ligands (Kateb et al., 2013).
The structural basis of SMYD MYND domains-mediated interactions is unknown. Based on the structural superposition (Figure 3D), the peptide binding mode in SMYD proteins could be largely similar to that of ETO. The leucine and first proline of the peptide could be similarly bound since the residues responsible for binding to these two positions are conserved and well-structurally aligned between SMYD proteins and ETO. In SMYD3, three conserved residues Trp80, Gln76, and Tyr70 may contribute to the binding of the PXLXP motif. The tryptophan residue may pack against the first proline, and the glutamine and tyrosine residues may form a hydrophobic pocket for the leucine to bind. However, due to a sequence variation between the third and fourth zinc-chelating cysteines there are structural differences in the peptide binding sites. Such differences might be associated with the different binding specificities. In ETO a glycine residue between the third and fourth zinc-chelating cysteines is replaced by a glutamine residue in SMYD proteins. As a result, their structures differ in the peptide binding pocket that binds to the peptide region between the leucine and first proline. Since the PXLXP motif bound by SMYD proteins is one residue shorter than the ETO-bound PPPLI sequence at this region, this pocket may contribute to binding specificity. The structural difference in this pocket may only allow the corresponding peptide motif to be accommodated. To date there are no structures available for PXLXP-complexed MYND domains and the exact binding mode for this motif is unknown.
In contrast to binding specificity, MYND domains also show promiscuous binding. Many MYND domains are capable of binding to multiple targets, and single peptide motifs can be recognized by different MYND domains. The MYND domain of BS69 is able to bind to MGA, E1A, and EBNA2 proteins (Ansieau and Leutz, 2002), while binding to SMRT and N-CoR appears to be a common feature of many MYND domains including those of ETO, RACK7, and DEAF-1 (Liu et al., 2007;Kateb et al., 2013). Additionally, the promiscuous binding can be due to the coexistence of several proline-rich motifs in one target. N-CoR contains the overlapping PPPLI and PLXLP motifs and can then be bound by different classes of MYND domains (TRPPPPLIPSSK). The structure of ETO-SMRT complex has revealed limited specific interactions between the MYND domain and peptide (Figure 3C). The less constrained structural features in the peptide binding site would allow for a broader binding specificity providing some structural explanation of promiscuous MYND binding.

Enrichment of PXLXP Sequences in SMYD3 Interactors
Motif scanning using ScanProsite (de Castro et al., 2006) reveals a significant enrichment of PXLXP motif in SMYD3 interacting proteins. A total of 16 matches was found in 10 proteins including PLCB3, CAMK2B, PKN1, ESR1, HELZ, KDM3B, KMT2E, CSTF2T, MEST, and MAP3K2. In other words, 33.3% of SMYD3 interacting proteins contain 1 or more PXLXP motifs, and the match percentage is 53.3% (Figure 4A). The chance for a random match for this motif is 10.3% (Nicodeme, 2001). This indicates more than 5.1-fold enrichment of PXLXP motif in SMYD3 interacting proteins. Sequence alignment of all 16 matches reveals a specific PQLSP motif occurs three times each in ESR1, HELZ, and MEST ( Figure 4B). This occurrence represents more than 330-fold enrichment. Of note, the flanking sequences of this motif are also conserved containing additional proline residues (PPPQLSPF).
The significant enrichment of PXLXP motif in SMYD3 interacting proteins suggests that the MYND domain is likely to be a key interaction module that mediates the formation of some SMYD3 complexes. The presence of PXLXP motifs in PLCB3 and CAMK2B indicates the potential role of the MYND domain in mediating complex formation in signaling including the aforementioned estrogen signaling pathway, GnRH signaling pathway, NOD-like receptor signaling pathway, and WNT signaling pathway. In PLCB3 there are two PXLXP motifs, PPLRP at residues 265-269 and PSLEP at residues 530-534 ( Figure 4B). The former motif is located within the EF-hand repeats, and the latter motif is located at the X-Y linker connecting the X and Y halves of the triose phosphate isomerase (TIM)-like barrel domain ( Figure 4C). As shown by the crystal structure of PLCB3, the latter motif appears to be more accessible for binding than the former motif ( Figure 4D). The central leucine (L) in the former motif is buried inside the EF-hand domain and involved in intramolecular interactions. With its current conformation it is essentially impossible for this motif to interact with MYND domains unless there is a significant conformational change. In contrast, the X-Y linker where the latter motif is located is unstructured and likely located on the top of the protein surface. With this flexibility and location, the latter motif should be more accessible for binding with minimum steric hindrance. Additionally, the latter motif appears to be electrostatically more compatible with the MYND domain of SMYD3. The latter motif, PSLEP is negatively charged, whereas the former motif, PPLRP positively charged. The MYND domains of SMYD proteins have a highly positively charged surface and likely prefer binding to a negatively charged peptide (Spellmon et al., 2015). If the MYND domain binds to the latter motif, this binding could potentially regulate PLCB3 phospholipase activity and subsequently signal transduction processes, since the X-Y linker has been shown to play an autoinhibitory role (Hicks et al., 2008). CAMK2B has only one PXLXP motif, PCLSP at residues 476-480. This motif is located at the variable segment of CAMK2B ( Figure 4C). The variable segment of CAMK2B plays a role in targeting the holoenzyme CAMK2A/2B to Factin cytoskeleton and localizing the holoenzyme to dendritic spines (Shen et al., 1998). The targeting occurs at basal Ca 2+ concentrations; as intracellular Ca 2+ concentrations rise, the holoenzyme is released from the cytoskeleton and enters the cytosol (Swulius and Waxham, 2008). Binding to the cytoskeleton limits protein diffusion and has been thought as a mechanism to restrict the activation of the holoenzyme and subsequently the phosphorylation of its substrates (Swulius and Waxham, 2008). Since the PCLSP motif is located within the variable segment of CAMK2B, if this motif can be bound by the MYND domain of SMYD3, such binding could potentially regulate CAMK2B cytoskeletal targeting and the activation of the holoenzyme. The PCLSP motif may have a physiological role in signaling, as the serine residue in this motif is phosphorylated in the signaling pathways downstream of SDF-1/CXCR4 in breast cancer stem cells (Yi et al., 2014). However, the PCLSP motif is not present in mouse and rat proteins and also missing in other CAMK2 isoforms (UniProt Consortium, 2019). This unique feature might be relevant to the role of CAMK2B in human-specific synaptic plasticity and the functional differences between CAMK2A and CAMK2B in F-actin binding (Shen et al., 1998).

DISCUSSION
Analysis of SMYD3 interactome reveals SMYD3 might be involved in signal transduction pathways in addition to epigenetic gene regulation. Most enriched KEGG pathways are associated with the signaling processes in cancer. Among them, the WNT signaling has been shown to affect the maintenance of stemness and metastasis of cancer cells (Zhan et al., 2017), and the MAPK signaling can enhance cell proliferation and angiogenesis (Dhillon et al., 2007). This indicates that SMYD3 might adopt an integrated model that combines signaling and epigenetic pathways to contribute to tumor cell proliferation and growth. One notable interaction is with the phospholipase PLCB3. PLC enzymes represent convergence points for many signal transduction pathways. PLCB3 is known to regulate opioid-dependent signaling and is required for opioid-induced calcium release in the nervous system (Xie et al., 1999). PLCB3 also plays a role in Pseudomonas aeruginosa-induced intracellular calcium signaling and regulates the activation of PKC and the nuclear transcription factor NF-κB in human bronchial epithelial cells (Bezzerri et al., 2011). The interaction with PLCB3 would suggest that SMYD3 might be involved in PLCB3-dependent signaling pathways to regulate free calcium transients. One key question needed to be addressed is where the interaction between SMYD3 and PLCB3 occurs in cell. SMYD3 is present in both cytosol and nucleus, whereas PLCB3 is primarily a peripheral membrane protein but also found in nucleus (Thul et al., 2017). The translocation of PLCB3 from the nucleus to the plasma membrane is a controlled process that regulates its activation and subsequently IP3-induced calcium signaling in Jurkat Tcells (Xie et al., 1999). If the interaction occurs in the nucleus, SMYD3 could potentially affect the nuclear export of PLCB3. SMYD3 has been shown to regulate nucleo-cytoplasmic shuttling of Human T-lymphotropic Virus 1 (HTLV-1) Tax protein by direct interaction with Tax and tethering Tax in the nucleus of T cells (Alefantis et al., 2005). Therefore, knowledge of spatial interaction between SMYD3 and PLCB3 at the subcellular level would elucidate more detailed mechanistic insights into the role of SMYD3 in PLCB3-dependent signaling pathways in terms of regulation of PLCB3 nuclear export or its phospholipase activity at the plasma membrane.
The MYND domain might be a key complex-forming module in SMYD3 complexes. The PXLXP motif to which the MYND domain may bind is significantly enriched in SMYD3 interacting proteins. However, there is no significant difference in the enrichment between the interactors involved in signaling transduction pathways and epigenetic transcriptional regulation suggesting that the MYND domain does not have any particular binding preferences in either of these processes. No MYND domain has been directly linked to signaling, but some MYND domain-containing proteins such as BRAM1 as well as its Caenorhabditis elegans orthologs, Bra-1 and Bra-2, have been shown to participate in BMP and TGFβ signaling pathways (Kurozumi et al., 1998;Morita et al., 2001). For SMYD3, its MYND domain could be involved in PLCB3mediated signaling pathways since PLCB3 contains the PXLXP motifs. The probable MYND binding motif in PLCB3, PSLEP is located within the X-Y linker, a sequence region highly variable across PLCB isozymes and also across species. The PSLEP motif is not present in other PLCB isozymes while it shows a modest degree of conservation across species (UniProt Consortium, 2019). The motif is conserved in most vertebrate animals including higher primates, horse, dog, sheep, and some rodents. The motif is not conserved in bird, fish, reptile, some rodents including mouse and rat, some lower primates such as galago, and all invertebrate animals. The modest conservation or divergence of the PSLEP motif might be related to the regulatory roles of the X-Y linkers. The very divergent primary sequences of the various X-Y linkers are likely the results of lineage-or species-specific evolution that creates distinct modes of autoregulation, therefore helping to respond to multiple extra-and intracellular inputs with appropriate enzymatic kinetics. The structural basis of SMYD3 MYND domain-mediated interactions is unknown. The only available MYND-peptide complex structure suggests the MYND domain-mediated interaction is likely transient and dynamic (Liu et al., 2007). Transient interactions are frequently involved in signaling that allow rapid responses to cellular perturbations and changes in environment (Stein et al., 2009;Perkins et al., 2010). The transient nature of MYND-domain mediated interactions would make this domain well-suited to mediate signaling transduction processes.

Selection of High-Confidence SMYD3
Interacting Proteins SMYD3 interacting proteins were retrieved from four proteinprotein interaction databases including GPS-Prot (Fahey et al., 2011), BioGRID (Stark et al., 2006), HitPredict (Lopez et al., 2015), and STRING (Szklarczyk et al., 2015). GPS-Prot archives 38 experimentally determined interactors from human, BioGRID 29, HitPredict 33, and STRING 15; a total of 43 unique interactions. The final list of proteins used in this study was generated by removing interactors with low-confidence reliability scores. We considered proteins as the low-confidence interactors if their interaction score is <0.3 in STRING or labeled as "low" in HitPredict. The final list contains 30 interactors.

Pathway Enrichment Analysis
Pathway enrichment analysis of SMYD3 interactors was performed using STRING web-based analyzer (Szklarczyk et al., 2015). Pathway resources used are the Gene Ontology (GO) knowledgebase (The Gene Ontology Consortium, 2019) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2019). All human proteins were used as the background population in the analysis. Evidence for enrichment in annotated pathways was evaluated using the Benjamini and Hochberg procedure that controls the false discovery rate (FDR) for multiple comparisons. An enrichment was considered as significant if the FDR is <0.05.

PXLXP Motif Enrichment Analysis
The presence of PXLXP motif in SMYD3 interacting proteins was scanned using ScanProsite (de Castro et al., 2006). Enrichment analysis was performed by comparing the number of matches to the expected random matches of the motif in the background, i.e., ∼100,000 sequences or 50,000,000 residues (Nicodeme, 2001). Fold enrichment was calculated as the ratio of the match percentage within our protein list to the percentage of the random matches against the background. The expected random matches are 10318 for PXLXP motif, and 30 for PQLSP motif (Nicodeme, 2001). Significance of enrichment was evaluated using the Benjamini and Hochberg procedure. FDR <0.05 was considered as significant.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
ZY and YZ contributed conception and design of the study, performed the database and statistical analyses, and prepared the figures. ZY wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was supported by the National Institutes of Health grant number 1R01HL128647-01 (to CL).


text:
Long non-coding RNA colon cancer-associated transcript 2 may promote esophageal cancer growth and metastasis by regulating the Wnt signaling pathway

The aim of this study was to investigate how long non-coding (lnc)RNA colon cancer-associated transcript 2 (CCAT2) regulates the proliferation, invasion and metastasis of esophageal cancer cells via the Wnt signaling pathway. The expression of lncRNA CCAT2 was quantified by reverse transcription-quantitative PCR in four esophageal cancer cell lines (Eca-109, EC9706, KYSE150 and TE-1) and normal human esophageal epithelial cells (HEECs). The effect of silencing CCAT2 (si-CCAT2) and inhibiting Wnt signaling (using the inhibitor FH535) on the proliferation, migration and invasion of Eca-109 cells was measured by MTT, wound-healing and Transwell invasion assays. Flow cytometry was used to evaluate apoptosis in si-CCAT2 Eca-109 cells. The expression of β-catenin and proliferating cell nuclear antigen (PCNA) proteins was detected by immunohistochemistry. The pro-apoptotic protein Bax, cyclin D1 and Wnt target proteins, including c-Myc and adenomatous polyposis coli (APC), were detected by western blotting. LncRNA CCAT2 was highly expressed in the four esophageal cancer cell lines compared with the HEEC cells. The expression of CCAT2 was significantly decreased in si-CCAT2 Eca-109 cells. Treatment with si-CCAT2 and FH535 alone or in combination significantly inhibited the proliferation, migration and invasion of Eca-109 cells. The treatments also promoted apoptosis, upregulated the expression of Bax and APC proteins, and downregulated β-catenin, PCNA, cyclin D1 and c-Myc proteins. In summary, lncRNA CCAT2 is upregulated in esophageal cancer cells and the knockdown of lncRNA CCAT2 inhibits their proliferation, migration and invasion via the Wnt signaling pathway.

Introduction
Esophageal cancer, including squamous cell carcinoma and adenocarcinoma, affected half a million people worldwide in 2014, and causes a fatal outcome in the majority of cases (1,2). The incidence rate of esophageal cancer has increased dramatically over the past 40 years (2). Alcohol consumption and tobacco use are major risk factors (3)(4)(5). Untreated gastroesophageal reflux disease usually progresses to Barrett's esophagus, which is 50-100 times more likely to develop into cancer (1). Obesity can also increase the risk of developing esophageal cancer, via low-grade inflammation (2). Due to the lack of serosa in the esophageal wall and extensive lymphatic drainage to intrathorax and intra-abdominal lymph nodes, esophageal cancer can spread rapidly into the neck and thorax (2). Advanced stages of esophageal cancer can also spread hematogenously (2).
LncRNA CCAT2 was found to be associated with ESCC by in silico analysis and subsequent clinical studies (4,19). These studies showed that CCAT2 was highly expressed in cancer tissue and was associated with smoking (4,19). Lymph node metastasis, advanced lymph node metastasis and Myc amplification were also associated with high lncRNA CCAT2 expression (19). Furthermore, the expression of this lncRNA was positively correlated with Myc amplification and progression of cancer (19).
To investigate the molecular mechanism of lncRNA CCAT2 function in esophageal cancer, the expression of lncRNA CCAT2 in esophageal cancer cells and its association with proliferation and metastasis were investigated in the present study. The results from the current study may provide a theoretical basis for new treatment options for esophageal cancer. Grouping. Based on the expression level of lncRNA CCAT2 in the four esophageal cancer cell lines, the Eca-109 cell line was selected for further experiments because it exhibited a higher expression level of lncRNA CCAT2 than the other three cell lines. The cells were divided into 5 groups: i) Control, no treatment; ii) negative control group (si-NC), scramble sequence transfected; iii) lncRNA CCAT2-silenced group (si-CCAT2), transfection with siRNA sequences targeting CCAT2; iv) Wnt pathway inhibitor (FH535) group, 10 µM FH535 treatment (cat. no. HY-15721; MedChemExpress); and v) CCAT2 silencing and inhibitor group (si-CCAT2 + FH535), transfection with siRNA sequence targeting CCAT2 and 10 µM FH535 treatment. The esophageal cancer Eca-109 cells were digested and passaged with 0.25% trypsin (Invitrogen; Thermo Fisher Scientific, Inc.). Cells (2x10 5 ) were inoculated on a 6-well plate. After 24 h, when the cells reached 30-50% confluence, 10 µl/250 µl siRNA was transfected with Lipofectamine ® 2000 Transfection Reagent (Invitrogen; Thermo Fisher Scientific, Inc.). A further ~24 h later, the medium containing the mixture was preplaced with fresh medium containing 10% FBS. Cell transfection was assessed by reverse transcription-quantitative PCR (RT-qPCR). The sequences of the control and lncRNA CCAT2 siRNAs were as follows: si-CCAT2 sense strand, 5'-GCC UGU AGG AAG AGU CAA ATT-3'; si-CCAT2 antisense strand, 5'-UUU GAC UCU UCC UAC AGG CTT-3'; si-NC sense strand, 5'-UUC UCC GA A CGU GUC ACG UT T-3'; and si-NC antisense strand, 5'-ACG UGA CAC GUU CGG AGA ATT-3'.

Materials and methods
RT-qPCR. Total RNA was extracted from each cell line using the Total RNA Extraction kit (Invitrogen; Thermo Fisher Scientific, Inc.), and RNA was reverse-transcribed into cDNA using SuperScript III Reverse Transcriptase (Thermo Fisher Scientific, Inc.). SYBR Green PCR kit (Qiagen, Inc.) and Mastercycler ® Ep Realplex2 (Eppendorf) were used for RT-qPCR. The RT-qPCR was performed using 2 µl cDNA as a template under the following conditions: 95˚C for 10 min, 95˚C for 15 sec and 60˚C for 1 min for 40 cycles. GAPDH was used as the internal reference and the relative expression of lncRNA CCAT2 was calculated by the 2 -ΔΔCq method (20). The following primer pairs were used: CCAT2 forward, 5'-AGA CAG TGC CAG CCA ACC-3' and reverse, 5'-TGC CAA ACC CTT CCC TTA-3'; and GAPDH forward, 5'-AAT GGA CAA CTG GTC GTG GAC-3' and reverse, 5'-CCC TCC AGG GGA TCT GTT TG-3'.
MTT assay. Cells (~100 µl; 1x10 4 cells/ml) were inoculated on a 96-well flat-bottomed culture plate and cultured for 24, 48, 72 and 96 h. MTT (100 µg) was added and incubated at 37˚C for 4 h. The supernatant was dissolved with 200 µl DMSO (Sigma-Aldrich; Merck KGaA). Optical density (OD) of each well was measured at a wavelength of 490 nm using a plate reader (Bio-Rad Laboratories). The relative absorbance of each group was calculated by the following equation: Wound-healing assay. Cells (~1.0 ml; 3x10 5 cells/ml) were plated in each well (6-well plate). When cells reached 100% confluence within 24 h, a 10-µl pipette tip was used to scratch the monolayer of cells. The cells were washed with PBS twice to remove suspended cells and serum-free medium was gently added. The cell-free wound area was photographed with an inverted microscope (Olympus Corporation) and the results were analyzed by ImageJ 1.46r software (National Institutes of Health). The percentage of the wound size at 24 h relative to the 0 h time point of each treatment group was calculated.
Transwell assay. Matrigel mixture (~50 µl) was added to the upper chamber of Transwell inserts (EMD Millipore) and incubated at 37˚C for 30 min. Cells (2x10 5 cells/ml) were resuspended in serum-free DMEM (Gibco; Thermo Fisher Scientific, Inc.). A total of ~100 µl cell suspension was added into the upper compartment of the chamber, and 600 µl of 10% FBS-DMEM culture media was added to the lower chamber. Membranes were collected after 24 h incubation, fixed in 5% glutaraldehyde at 4˚C for 30 min, stained with 0.5% crystal violet at room temperature for 20 min and observed under a light microscope (Olympus Corporation) at x200 magnification. Images of nine random fields were captured for analysis of the number of invading cells.
Flow cytometry. Cells were digested by trypsin and collected by centrifugation at 11,180 x g for 5 min at 4˚C. The cells were washed twice with pre-cooled sterile PBS solution at 4˚C. Cells [~195 µl; 1x10 6 cells in 250 µl 1X binding buffer (Beyotime Institute of Biotechnology)] were incubated with 5 µl of FITC-labeled annexin-V (Beyotime Institute of Biotechnology) for 3 min. Then, 10 µl propidium iodide solution (20 µg/ml) was added and incubated in the dark for 10 min at room temperature. Binding buffer (~400 µl 1X solution) was added prior to analysis using a flow cytometer (Gallios; Beckman Coulter, Inc.) and Cell Quest 5.1 software (BD Biosciences). The cell numbers at different stages were recorded and the ratio of apoptotic cells was defined as the ratio of cells in quadrants 2 and 4 to total cells.
Immunocytochemistry. Cells, cultured for 24 h, were washed with PBS two to three times and fixed with paraformaldehyde for 30 min. PBS was used to wash the cells two to three times prior to adding 0.5% Triton X-100 in PBS and incubated for 20 min. This was followed by two to three PBS washes before 3% H 2 O 2 was added for 15-20 min. Cells were then washed with PBS two to three times and 5% bovine serum albumin (Beyotime Institute of Biotechnology) was added and incubated at 37˚C for 30 min. The primary antibodies, including rabbit anti-human anti-β-catenin (1:200; cat. no. ab6302) or anti-PCNA (1:100; cat. no. ab15497) (both Abcam), were added and incubated at 37˚C for 30 min. After two to three PBS washes, horseradish peroxidase-conjugated goat anti-rabbit IgG (1:1,000; cat. no. ABIN101988; antibodies-online GmBH) was added at 37˚C for 30 min. Cells were washed with PBS two to three times before streptavidin-biotin complex (Beyotime Institute of Biotechnology) was added and incubated at 37˚C for 30 min. Following three washes with PBS, 3,3'-diaminobenzidine was added, followed by incubation with hematoxylin at room temperature for 2 min and washing with ethanol at increasing concentrations (70% alcohol for 1 min, 80% alcohol for 1 min, 95% alcohol for 2 min and absolute ethanol for 4 min). Observations were performed at x400 magnification under a light microscope (Olympus Corporation) and three to five fields of each section were randomly selected for analysis.
Western blotting. Cells were lysed in RIPA buffer (Beyotime Institute of Biotechnology) and the protein concentration was measured by the bicinchoninic acid (Pierce™ BCA Protein Assay kit; Thermo Fisher Scientific, Inc.). Each protein sample (~40 µg) was separated on a SDS-PAGE 10% gel (Mini-Protean-3 type; Bio-Rad Laboratories) and then transferred to a PVDF membrane (Merck KGaA) for 30 min. The membrane was blocked with 5% non-fat dry milk in TBST solution for 1 h. TBST with 3% bovine serum albumin was used to dilute each of the antibodies, including rabbit anti-human anti-Bax (cat. no. ab53154), anti-APC (ab15270, Abcam), anti-cyclin D1 (cat. no. ab226977), anti-c-Myc (cat. no. ab39688) and β-actin (cat. no. ab8227) (all 1:1,000) polyclonal antibodies, and goat anti-rabbit IgG (1:2,000; cat. no. ab6721) (all Abcam). The membrane  was incubated with primary antibody overnight at 4˚C. Following incubation with horseradish peroxidase-labeled goat anti-rabbit secondar y antibody (1:1,000; cat. no. ABIN101988; Antibodies Online) at room temperature for 1 h, the ECL system (Thermo Fisher Scientific, Inc.) was used to detect the signals. The expression levels were quantified by ImageJ 1.46r software, and β-actin was used as an internal control.
Statistical analysis. Each experiment was repeated at least three times. SPSS 19.0 (IBM Corp.) was used to analyze all the data, which are expressed as the mean ± SD. One-way analysis of variance was used for analysis of data among groups. Least-Significant Difference was used for post hoc analysis. P<0.05 was considered to indicate s statistically significant difference.

Results
Expression of lncRNA CCAT2 is higher in esophageal cancer cell lines. The expression of lncRNA CCAT2 in four esophageal cancer cell lines (Eca-109, EC9706, KYSE150 and TE-1) and one normal esophageal epithelial cell (HEEC) was examined by RT-qPCR. The expression was significantly higher in all esophageal cancer cell lines compared with that in HEEC (P<0.05; Fig. 1A). Subsequently, siRNA was used to silence lncRNA CCAT2 in Eca-109 cells. The RT-qPCR data showed no off-target effect occurring, and the knockdown of lncRNA CCAT2 was successfully performed (Fig. 1B).

Inhibition of lncRNA CCAT2 attenuates proliferation, migration and invasion of esophageal cancer cells.
After 24, 48, 72, and 96 h of culture following transfection with si-CCAT2, proliferation decreased significantly in the si-CCAT2 group compared with the control and si-NC groups (P<0.05; Fig. 2). To further investigate the effect of lncRNA CCAT2 on the metastasis of esophageal cancer cells, the wound-healing and Transwell assays were performed. As shown in Fig. 3, the cell migration and invasion in the si-CCAT2 group were significantly suppressed compared with the control and si-NC groups (P<0.05).

Inhibition of lncRNA CCAT2 promotes apoptosis in esophageal cancer cells.
Flow cytometry was used to analyze apoptosis in esophageal cancer cells, following annexin V staining. The number of apoptotic cells in the si-CCAT2 group was significantly higher compared with the control and si-NC groups (P<0.05; Fig. 4). Similarly, the immunohistochemistry staining showed significantly decreased expression of β-catenin and PCNA in the si-CCAT2 group (P<0.05; Fig. 5). Protein expression, determined by western blotting, showed significant upregulation of Bax and APC in the si-CCAT2 group compared with the control and si-NC groups, whereas cyclin D1 and c-Myc proteins were significantly downregulated (P<0.05; Fig. 6).
LncRNA CCAT2 and Wnt signaling pathway. To further analyze how lncRNA CCAT2 functions in esophageal cancer cells, the effect of the Wnt signaling pathway on the cell proliferation, invasion and metastasis was investigated, using a Wnt signaling inhibitor (FH535). No differences in the effects on apoptosis, proliferation, migration and invasion of Eca-109 cells were observed between the si-CCAT2 + FH535 group and the si-CCAT2 or FH535 groups (Figs. 7-9). Expression of β-catenin, PCNA, Bax, APC, cyclin D1 and c-Myc were also similar among these groups (Figs. 10 and 11).

Discussion
β-catenin can bind T-cell factor (TCF) family DNA-binding proteins and translocate to the nucleus. This complex regulates cellular differentiation and proliferation (21)(22)(23)(24). Wnt signaling molecules can activate the β-catenin/TCF complex and induce gene expression (21)(22)(23). This signaling pathway has been identified to serve a major role in in certain cancer types including highly malignant carcinoma such as esophageal cancer (21,(25)(26)(27). It is believed that this pathway plays a critical role in the maintenance and self-renewal of cancer stem cells, which leads to malignant tumors and poor prognosis (22,25). Wnt signaling is positively correlated with higher ratio of cancer stem cells and advanced clinical stages of esophageal carcinoma (25). Other genes, such as Dapper homolog 2, naked cuticle homolog 2, microRNA (miR)-942, transcription factor Sox17 and miR-141, can modulate the  Wnt signaling pathway and affect tumorigenesis in esophageal cells (22,23,26,28). Human papillomavirus can lead to esophageal cancer via the miR-125b and Wnt/β-catenin signaling pathway (27).
The present study demonstrated that esophageal cancer cell lines had higher expression levels of lncRNA CCAT2 compared with normal esophageal cells. To further explore the function of lncRNA CCAT2, si-CCAT2 was constructed and tested in esophageal cancer Eca-109 cells. The results showed that si-CCAT2 decreased the cell proliferation. Furthermore, the population of apoptotic cells, quantified by annexin V staining, increased significantly following knockdown of lncRNA CCAT2. This downregulation also significantly decreased cell migration and invasion. To elucidate the molecular mechanism of lncRNA CCAT2, a Wnt inhibitor (FH535) was used to study the associated signaling pathway. FH535 and lncRNA CCAT2 equally elicited suppressive effects on proliferation, migration and invasion, suggesting that lncRNA may function via the Wnt pathway.
LncRNA CCAT2 has been demonstrated to be involved in the Wnt pathway in a number of cancer types, and one of its functions is to modulate the Wnt pathway (35)(36)(37)(38). For example, downregulating CCAT2 suppresses the transcriptional activity of Wnt/β-catenin signaling pathway and the combination of siCCAT2 and FH535 synergistically inhibits Wnt signaling in breast cancer (39). The present study demonstrates an association between lncRNA CCAT2 and Wnt signaling pathway in esophageal cancer.
In the present study, only one esophageal cancer cell line (Eca-109) was used to study the effects and mechanism of lncRNA CCAT2. Therefore, a limitation in the present findings is that they may not apply to other esophageal cancer cell lines or different esophageal cancer types. In addition, the present results showed that the expression level of lncRNA CCAT2 was different in the four esophageal cancer cell lines, which was not investigated any further in this study. Furthermore, more detailed studies of the regulatory effect of CCAT2 on the Wnt pathway should be conducted.  In the future, the correlation between the expression level of lncRNA CCAT2 and the proliferation and metastasis of the four esophageal cancer cell lines will be explored. In addition, the hypothesis tested in the present study will be applied in an  in vivo model. These studies will provide more critical insights into the mechanism, for the development of novel treatments for esophageal cancer in the future.
In summary, the present study demonstrated that lncRNA CCAT2 was highly expressed in esophageal cancer cells, and its downregulation led to the inhibition of cell proliferation,