Laboratory Evolution of GH11 Endoxylanase Through DNA Shuffling: Effects of Distal Residue Substitution on Catalytic Activity and Active Site Architecture

Endoxylanase with high specific activity, thermostability, and broad pH adaptability is in huge demand. The mutant library of GH11 endoxylanase was constructed via DNA shuffling by using the catalytic domain of Bacillus amyloliquefaciens xylanase A (BaxA) and Thermomonospora fusca TF xylanase A (TfxA) as parents. A total of 2,250 colonies were collected and 756 of them were sequenced. Three novel mutants (DS153: N29S, DS241: S31R and DS428: I51V) were identified and characterized in detail. For these mutants, three residues of BaxA were substituted by the corresponding one of TfxA_CD. The specific activity of DS153, DS241, and DS428 in the optimal condition was 4.54, 4.35, and 3.9 times compared with the recombinant BaxA (reBaxA), respectively. The optimum temperature of the three mutants was 50°C. The optimum pH for DS153, DS241, and DS428 was 6.0, 7.0, and 6.0, respectively. The catalytic efficiency of DS153, DS241, and DS428 enhanced as well, while their sensitivity to recombinant rice xylanase inhibitor (RIXI) was lower than that of reBaxA. Three mutants have identical hydrolytic function as reBaxA, which released xylobiose–xylopentaose from oat spelt, birchwood, and beechwood xylan. Furthermore, molecular dynamics simulations were performed on BaxA and three mutants to explore the precise impact of gain-of-function on xylanase activity. The tertiary structure of BaxA was not altered under the substitution of distal residues (N29S, S31R, and I51V); it induced slightly changes in active site architecture. The distal impact rescued the BaxA from native conformation (“closed state”) through weakening interactions between “gate” residues (R112, N35 in DS241 and DS428; W9, P116 in DS153) and active site residues (E78, E172, Y69, and Y80), favoring conformations with an “open state” and providing improved activity. The current findings would provide a better and more in-depth understanding of how distal single residue substitution improved the catalytic activity of xylanase at the atomic level.

As a biological catalyst, the catalytic activity and efficiency of xylanases are susceptible to be influenced by both external and internal factors. The most common external factors include temperature, pH, substrate, and inhibitors, while the internal factors include primary structure and crystal structure. Xylanases with improved properties (i.e., higher activity and stability) are urgently welcomed and needed for environmental, economic and industrial reasons (Pandey et al., 2016;Bala and Singh, 2019). It was reported that numerous strategies were applied to modify enzyme based on the structure-function relationship and included rational, semi-rational, and irrational design (directed evolution) (Farinas et al., 2001;Wong and Schwaneberg, 2003;Chica et al., 2005;Sarmiento et al., 2015, Cheng et al., 2015. Therefore, lots of researches confirmed that these strategies mentioned above could be adopted to enhance the catalytic activity and enzymatic properties of xylanase (Shibuya et al., 2000;Miyazaki et al., 2006;Stephens et al., 2007;Wang et al., 2013;Wahab et al., 2016;Prajapati et al., 2018;Damis et al., 2019).
amyloliquefaciens Xylanase A (BaxA) is mesophilic and less thermostable. The activity of TfxA is lower than that of plenty of GH11 xylanases such as BaxA and Aspergillus niger xylanase A (AnxA) (Xu et al., 2016). The sequence similarity index results revealed that the catalytic domain of TfxA is nearly 40% identical with that of BaxA. In our previous study, we have modified the baxA gene (KM624029), encoding the BaxA by using error-prone PCR (epPCR) (Xu et al., 2016). The mutant reBaxA50 (S138T, including signal peptide) with improved catalytic activity was screened and characterized.
DNA shuffling, as one of rapid and powerful techniques for directed evolution of enzymes, was discovered by Stemmer in 1994(Stemmer, 1994a. In the study, the mutant library was constructed via DNA shuffling by using the catalytic domain (CD) of TfxA and BaxA as parents. Moreover, three mutants (DS153, DS241, and DS428) were identified, and the variation in the enzymatic properties between reBaxA and mutants were analyzed in detail. Additionally, we have explored the activation mechanism for the first time, induced by GOF-mutations using essential dynamics simulation approaches. Subsequently, noncovalent interactions (NCI) analysis was also carried out to deepen understanding of modeling distinctions at the atomic level. Exploring effects of distal residue substitution on active site architecture may assist in providing valuable information for better clarifying the mechanism of enhanced activity.

Construction of DNA Shuffling Mutant Library
The mutant library was generated as described by Shibuya with slight modifications (Shibuya et al., 2000). DNA fragments (baxA and tfxA_cd, about 600 bp) were amplified from pCold TF-baxA and pCold I-tfxa_cd plasmids with primers DNAScold1: 5 ′ -TCGGTACTCTCGAAGGTTTCGAATTC−3 ′ and DNAScold2: 5 ′ -GTCCTTTTAAGCAGAGCTTACTATCTAGA-3 ′ , which contained EcoRI and XbaI recognition sites (underlined), respectively. The PCR parameters were as followed: 95 • C, 2 min; 32 cycles of (94 • C, 50 s; 61 • C, 50 s with a 0.2 • C decline per cycle; and 72 • C, 1 min); and 72 • C, 10 min. The PCR products were purified and recovered. Two genes (baxA, 3 µg; tfxA_cd, 3 µg) were mixed and randomly digested by DNaseI (0.2 U) at 37 • C for 15 min. The digested products were analyzed in 2.0% (w/v) agarose gel, and then DNA fragments with the molecular mass <100 bp were recovered and purified. The gel-recovered fragments were reassembled into full-length genes via primerless PCR under the following conditions: 94 • C, 150 s; 40 cycles of (30 s at 94 • C, 47.5 s at 61 • C, 4 min at 72 • C); 72 • C, 10 min. The primerless PCR product was analyzed by electrophoresis and used as the template to amplify the single product of full-length (about 600 bp) by using the primers DNAScold1 and DNAScold2. Moreover, the full-length gene was purified and digested by EcoRI and XbaI. The digested full-length gene was inserted into pCold TF vector and then transformed into E. coli BL21 (DE3) competent cells. All of the transformed products were plated on lysogenic broth (LB) agar plates (ampicillin, 100 µg/mL) and cultured at 37 • C for 12 h.

Screening of the Mutant Library and Activity Assay
The colonies on LB agar plates (ampicillin, 100 µg/mL) were inoculated in 5 mL of LB medium (ampicillin, 100 µg/mL) at 37 • C for 12 h and then 1 mL cultured cells were transferred into shake flask (50 mL of LB medium with 100 µg/mL ampicillin) continuously cultured at 37 • C for 3 h. The shake flask was stored at 15 • C for 30 min. In order to induce expression of xylanase, isopropyl-β-D-thiogalactopyranoside (IPTG) was added to the medium with a final concentration of 1.0 µM, and the cells were cultured under shaking (150 rpm) at 15 • C for 24 h.
The fermentation supernatant was collected and used for enzyme activity assay in a 96-well plate as follows. First, 50 µL of 0.5% beechwood xylan and 20 µL of fermentation supernatant were co-incubated at 50 • C for 5 min. Then, 50 µL of 3,5-dinitrosalicylic acid (DNS) was added to stop the reaction and placed in a boiling water bath for 5 min. After 100 µL of deionized water was added, absorbance was determined spectrophotometrically at 540 nm by using a multiscan spectrum microplate spectrophotometer (Xu et al., 2016). These colonies that showed xylanase activity were sequenced. Mutant DS153, DS241, and DS428 were characterized in details.
One unit of xylanase activity was defined as the amount that released 1 µmol of reducing sugar equivalents (D-xylose as the standard) per minute from beechwood xylan under optimal conditions (Miller, 1959;Miller et al., 1959). Protein concentration was measured by Bradford method with bovine serum albumin as the standard (Bradford, 1976). The Michaelis-Menten constant (K m ), the maximum reaction rate (V max ), and molar catalytic activity (k cat ) of xylanase were determined from nonlinear regression fitting to the Michaelis-Menten equation by using beechwood xylan as substrate (1, 3, 5, 7.5, and 10 mg/mL). All assays were performed in triplicate to obtain the mean and standard deviation.

SDS-PAGE and Western Blotting Analysis
The xylanase with high activity was purified using high-affinity Ni-charged resin from the culture supernatant. The samples were analyzed by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) with the stacking and separating gels consisting of 5% and 12% polyacrylamide, respectively (Laemmli, 1970). Proteins in the gel were stained with Coomassie brilliant blue R-250. Anti-His monoclonal and horseradish peroxidase-labeled goat anti-mouse IgG antibodies were used for Western blotting analysis (Huo et al., 2018).

Optimum Temperature and Thermostability
The effect of temperature on the xylanase activity was determined from 30 to 90 • C. For thermal stability assay, the xylanase was heat-treated from 30 to 90 • C for 15 min and then placed in an ice-water bath for 5 min. Residual activity was determined under optimal conditions.

Optimum pH and pH Stability
The effect of pH on xylanase activity was evaluated from pH 3.0 to 9.0 at 50 • C. The highest activity was taken as 100%. For pH stability, xylanase was incubated in buffer solutions (pH 3.0-9.0) at 25 • C for 1 h. Residual activity was measured under optimal conditions. For assays of pH stability and thermostability, the activity of the untreated xylanase under optimal conditions was taken as 100%.

Inhibitory Activity of Recombinant Rice Xylanase Inhibitor on Xylanase
Rice xylanase inhibitor (RIXI) is one of the XIP-type xylanaseinhibiting proteins in rice, which involved in protecting rice from attacks by plant-pathogenic microorganisms and inhibiting most exogenous GH10 and 11 xylanases (Goesaert et al., 2005). The content of xylanase-inhibiting proteins in rice and wheat is relatively high, and theses inhibitors reduced the application performance of xylanase in the feed and food industries (Dornez et al., 2010). According to our previous study (Huo et al., 2018), we determined inhibitory activity of recombinant rice xylanase inhibitor (rePRIXI) on reBaxA, TfxA_CD, and three mutants (DS153, DS241, and DS428). Xylanases (0.4 U) were incubated with rePRIXI (0.2 mg) at 50 • C for 40 min, and then the residual activity and the inhibitory rate were calculated.

All-Atom Molecular Dynamics Simulations
Appling the three-dimensional (3D) structure of GH11 B. subtilis 1A1 (PDB ID 1XXN) as the template, the structural coordinates of BaxA and mutants (DS153, DS241, and DS429) were obtained via homology modeling using the SWISS-MODEL server (www. swissmodel.expasy.org) (Waterhouse et al., 2018), followed by energy minimization using Discovery Studio 3.5. The sequence identity of BaxA with the template was 96.22%.
The stability of wild-type and mutant xylanases were evaluated through extensive molecular dynamics simulations using the Amber14 package with the ff14SB force field (Case et al., 2018;Rehman et al., 2019). The LEAP module was used to add hydrogen atoms to the crystal structures. Counterions (Na + and Cl − ) were added for maintenance of system neutrality.
The system was first subjected to 10,000 steps of steepest descent energy minimization followed by 1,000 cycles of conjugate gradient minimization with bonds involving hydrogen constrained by SHAKE algorithm (Ryckaert et al., 1977), and then another 10,000 steps of steepest descent energy minimization followed by 5,000 cycles of conjugate gradient minimization with no constraint exerted. Then, the system was then gradually heated from 0 to 300 K through 25,000 iterations. After a 200ps-equilibrium in NPT ensemble, five 50-ns simulations (300 K, 1 atm) with different random seeds were conducted. The VDW interactions were cut off at 10 Å, and long-range electrostatic interactions were calculated with particle mesh Ewald (PME) method (Darden et al., 1993). The GPU supported pmemd code performed the production step of molecular dynamics simulation for each system (Gotz et al., 2012). Analyses of trajectories were performed using cpptraj in Ambertools14.

Non-covalent Interactions (NCI) Analysis of Substituted Residue and Active Sites
In this study, the quantum theory of atoms in molecules (QTAIM) (Bader, 1990) analysis, independent gradient model (IGM) (Lefebvre et al., 2017) analysis, reduced density gradient (RDG) (Johnson et al., 2010) analysis, and molecular electrostatic potential (MESP) were carried out with the Multiwfn 3.6 program (Lu and Chen, 2012). Molecular plots were visualized by the VMD 1.9.3 program (Humphrey et al., 1996). In atoms in molecules (AIM) analysis, the points at where the gradient norm of function value is zero (except at infinity) were called critical points (CPs). The appearance of a (3, −1) type of critical point usually appeared on a bond path or between the atoms which had attractive interaction, hence called as bond critical point (BCP). According to the AIM theory, the electron density (ρ) and its Laplacian value (∇ 2 ρ) were used to indicate the strength and nature, respectively. RDG analysis is a prevalent method to reveal weak interlayer interactions. The definition of the RDG function is shown below; it is essentially a dimensionless form of electron density gradient norm function: It is based on the RDG as a function of sign(λ 2 )ρ, where the sign(λ 2 ) is the sign of the second eigenvalue of the electron density Hessian matrix. Negative values of sign(λ 2 )ρ indicate attractive interactions, whereas positive values suggest repulsive interaction. In the RDG analysis, the isosurfaces were colored according to the value of sign(λ 2 )ρ, where a BGR (blue-greenred) color scale was adopted. Blue color represents attractive or bonding interaction, green weak van der Waals interaction, and red repulsive interaction. All isosurfaces are colored according to a BGR scheme over the electron density range −0.05 a.u. < sign(λ 2 )ρ < 0.05 a.u.. Lefebvre et al. proposed a handy way of visually studying interfragment and intrafragment interactions, named IGM. The method employs, in addition to ρ, quantities related to the first and second derivatives of the density. The IGM descriptor δg inter is given by the difference between the first derivatives of the charge densities for the total system and the fragments:

Construction and Screening of Mutant Library
The digested products with a size of <100 bp were assembled through primerless PCR, and ∼600 bp DNA was obtained (Figure 1). The following product of primerless PCR was further used as a template for amplification using both DNAScold1 and DNAScold2 as a primer. The PCR products were ligated into pCold TF vector and transformed into competent E. coli BL21 (DE3) cells. A total of 2,250 colonies were collected from the LB plates (ampicillin, 100 µg/mL). The mutant library was screened via 96-well plate method, and 756 colonies were sequenced and analyzed. The specific activities of DS153, DS241, and DS428 were 11.95, 11.45, and 10.26 U/mg, respectively. And they were 4.54, 4.35, and 3.9 times as high as that of reBaxA, respectively. The sequencing results for DS153, DS241, and DS428 revealed that the residues Asn29, Ser31, and Ile51 of BaxA were substituted by Ser38, Arg40, and Val59 of TfxA_CD correspondingly ( Figure S1). There are not only many single site mutations in these sequenced mutants, but some fragment replacement mutations.

SDS-PAGE, Western Blotting, and Kinetic Constant Analysis
The pCold TF vector system provides a cold shock technology and a folding chaperone trigger factor for high-level expression of soluble recombinant proteins. The cold shock technology, which requires decreasing the culture temperature from 37 to 15 • C during the induction period can suppress the expression of other cellular proteins and temporarily halt overall cell growth (Esposito and Chatterjee, 2006;Bjerga and Williamson, 2015). Xylanases were both secreted into the culture medium and remained in the cytoplasm of the three mutants after being induced by IPTG at 15 • C. The results of SDS-PAGE and Western blotting unveiled a major band existed in the supernatant of the culture, and the proportion of the recombinant protein was high. The DS153, DSA241, and DS428 were further purified by Ni-affinity resin from culture supernatant with molecular weight approximately equal to 74.2 kDa (Figure 2). (B) Digestion of baxa and tfxa_cd gene by DNaseI. The digestion products in the yellow box were recovered. (C) Primerless/primer PCR product. Lane 1 was primerless PCR product by using the recovered digestion product as template. Lane 2 was PCR product by using primerless PCR product (lane 1) as template and DNAScold1, DNAScold2 as primers.
Moreover, the reaction rates of the three mutants and reBaxA on beechwood xylan with different concentrations were plotted on Lineweaver-Burk graphs. Michaelis-Menten constant (K m ) of DS428 was lower than that of reBaxA, whereas the K m of DS153 and DS241 was slightly higher than that of reBaxA (Table 1). An increase in the K m indicates the low affinity of xylanase for the substrate; it can be assumed that the alteration in affinity happens due to the substitutions of residues. The k cat /K m of DS153, DSA241, and DS428 was higher than that of reBaxA, suggesting the enhanced catalytic efficiency. Shibuya et al. constructed the xylanase mutant library by using S. lividans xylanase B (SlxB-cat) and T. fusca xylanase A (TfxAcat) as parents through DNA shuffling technique. Mutant Stx15 and Stx18 exhibited significant enhancement in thermostability at 70 • C and their optimum temperature were equal to that of TfxA-cat. Moreover, their catalytic activities were improved in comparison with the parental xylanases. Sequencing results revealed that the 33 amino acids in N-terminus of Stx15 were from TfxA-cat and the other residues were from SlxB-cat, and the 62 amino acids in N-terminus of Stx18 were from TfxAcat. Another mutant, Stx2 that had two residue substitutions (Q24P-Y170D) showed an enhancement in thermostability at 70 • C (Shibuya et al., 2000).

Optimum Temperature and Thermostability
The optimum temperatures of reBaxA and the three mutants (DS153, DS241, and DS428) were 55 and 50 • C, respectively ( Figure 3A). The reBaxA and three mutants showed high activity at low temperature (30-40 • C), suggesting that these enzymes might be suitable for monogastric animal feed additives. The DS153 and DS241 were more thermostable compared with reBaxA. The residual activities for reBaxA, DS153, DS241, and DS428 were 69.42, 69.80, 87.41, and 45.78% respectively ( Figure 3B) while treating them at temperature 50 • C for 15 min. The strengthened thermostability was probably attributed to the mutation of S31R in DS241. The previous research found that the mutation from Ser to Arg and Asn to Ser were the most common mutations. So far as the consequences of these mutations, GH11 xylanases enhanced their catalytic activity and thermostability (Satyanarayana, 2013;Ayadi et al., 2015). Apart from that, the mutations (T11Y, N12H, N13D, F15Y, Y16F, S62T, S138T, S144C, N198D, A217V, and L49P) were discovered to be involved in conferring thermostability and catalytic activity on xylanases, including SoxB, SlxB, AnxA, AnxB, and BaxA (Zhang et al., 2010;Wang et al., 2012Wang et al., , 2015Xu et al., 2016). Furthermore, incorporating the results of molecular dynamics simulations and visual inspection of protein modeling hierarchy revealed that the β-sheets of GH11 xylanases had high Arg contents and were more thermostable. The xylanase mutant MxylM4 with four Arg residues was obtained by substituting Ser and Thr with Arg residue. The T 1/2 of the mutant MxylM4 at 80 • C with the substrate increased from 130 to 150 min, and the T m increased to 6 • C (Satyanarayana, 2013).

Optimum pH and pH Stability
The optimum pH of DS153 and DS428 was 6.0, which was the same as that of reBaxA. The optimum pH of DS241 was 7.0 ( Figure 3C). GH11 xylanase catalyzes the hydrolysis step of xylan via the double-displacement mechanism, and two glutamic acid residues (E78 & E172 in BaxA) serve as a nucleophile and an acid/base catalyst. The optimum pH of GH11 xylanase was closely related to the residue (Asp or Asn) that formed hydrogen bonds with the acid/base catalyst residue (E172). Normally, the residue in xylanases with an optimum pH below 5.0 is Asp, whereas the residue in xylanases that are active under more alkaline conditions is Asn (Krengel and Dijkstra, 1996). The A. niger xylanase I and T. fusca TF xylanase A have optimum pH of 3.0 and 6.0, and their residues are D35 and N44, respectively. The B. circulans xylanase (BCX) and its mutant (N35D) have optimum pH of 5.7 and 4.6, respectively (Joshi et al., 2000). In BaxA, the N35 is the optimum pH-dependent residue, which is adjacent to the acid/base catalyst E172. The substitution of N29S in DS153 and S31R in DS241 did not affect the optimum pH; on the other hand, the substitution of I51V in DS428 changed its optimum pH from 6.0 to 7.0. Thus, it indicated that the DS153 and DS241 were more stable than reBaxA in the range from pH 3.0 to 9.0. The pH stability of DS428 at pH 3.0-9.0 was similar to that of reBaxA ( Figure 3D).
The loop of RIXI (residues 144-153) inserted in the catalytic groove of GH11 xylanase and the residue Y147 and R148 directly interacted with the two catalytic residues, including E78 and E172 in BaxA. Our previous studies revealed that the inhibition of RIXI on BaxA belonged to competitive inhibition , and the molecular structural basis for the inhibition of BaxA by RIXI was similar to that of wheat XIP-I (Payan et al., 2004). Effect of temperature on activity of xylanases. 0.5% beechwood xylan was used as substrate. The highest activity was taken as 100% and the corresponding temperature was the optimum temperature. (B) Thermostability of xylanases. 0.5% beechwood xylan was used as substrate. The activity of xylanase under the optimum conditions was taken as 100% in assay of thermostability, and then the residual activity was assayed.
(C) Effect of pH on activity of xylanases. 0.5% beechwood xylan in different pH buffer was used as substrate to explore effect of pH on activity. The highest activity was taken as 100% and the corresponding pH was the optimum pH. (D) pH stability. The activity of xylanase under the optimum conditions was taken as 100% in assay of pH stability, and then the residual activity was assayed. (E) Inhibitory activity of rePRIXI on xylanase. The rePRIXI (0.3 mg) was mixture with 5 types of xylanases (0.3 U) at 50 • C for 40 min, respectively. And then the residual activity of xylanase was measured under their corresponding optimal condition. The inhibitory rate was obtained. Triplicate experiments were conducted for each assay.

Hydrolytic Products Released From Xylans by Xylanase
The hydrolytic products released from three xylans (beechwood, birchwood and oat spelt xylan) by the three mutants were determined by HPLC (Figure S2). The hydrolysates from beechwood, birchwood and oat spelt xylan by DS153, DS241, and DS428 were xylobiose (X2), xylotriose (X3), xylotetraose (X4), and xylopentaose (X5), which showed resemblance to those of reBaxA. The main product released from beechwood xylan and birchwood xylan by the three mutants was X3 (Table 2). Interestingly, the major product released by DS153, DS241, and DS428 from oat spelt xylan was X5. After hydrolysis at 50 • C for 15 min, the concentrations of X3 from beechwood xylan by DS153, DS241, and DS428 were 1.729, 1.767, and 1.431 mg/mL, respectively ( Table 2), which were higher than that of reBaxA (1.168 mg/mL). The concentrations of X5 from oat spelt xylan by DS153, DS241, and DS428 were 0.678, 1.018, and 0.459 mg/mL, respectively ( Table 2). Among hydrolysates of oat spelt xylan, the concentration of X3 by DS241 was higher than that of DS153 and DS428. The main product from oat spelt xylan by reBaxA was X4, and the concentrations of X4 and X5 were 1.764 and 0.943 mg/mL, respectively.
The single substitution in the three mutants changed the major product of oat spelt and birchwood xylan, but did not change the kind of hydrolysates from the three xylans. The affinity of the three mutants toward beechwood and birchwood xylan was higher than that toward oat spelt xylan. Our previous study (Xu et al., 2016) revealed that reBaxA50 (S138T) obtained from epPCR mutant library has higher catalytic efficiency on beechwood and birchwood xylans than on oat spelt xylan.
XOs are used as functional additives in food and feed industries (Achary and Prapulla, 2011). X2-X5 are the key functional components of XOs. An in-vivo test on male Wistar rats has revealed that dietary XOs had positive influence on the development of aberrant crypt focal formation and metabolic abnormalities, which were related to colon cancer (Aachary et al., 2015). In an in-vivo test on laying hens, dietary XOs supplementation considerably increased egg production, improved calcium digestibility and eggshell quality, and decreased the levels of plasma glutamic-pyruvic transaminase, cholesterol and high-density lipoprotein (Li et al., 2017). The clinical observations in a pilot study have demonstrated that XOs could modify gut microbiota in healthy and prediabetic (Pre-DM) subjects. The XOs could remarkably decrease or reverse the increase in the abundance of Howardella, Enterorhabdus, and Slackia in healthy and Pre-DM subjects.
In the Pre-DM group, the XOs reduced the two hrs-insulin levels in the oral glucose tolerance test but did not affect body composition, serum glucose, and satiety hormones (Yang et al., 2015). For this reason, the XOs can effectively maintain the health of human and animals, respectively.

Structural Dynamics Features and Analysis of Underling Mechanisms
Molecular dynamics (MD) simulations have been widely used to study the conformational changes of protein at the atomic level. The conformational features of all the systems were analyzed using MD trajectories along with the incorporation of various statistical parameters (Lobanov et al., 2008). The compactness of protein was measured by the radius of gyration (Rg). If a protein is stably folded, it will likely maintain a relatively steady value of Rg, whereas it will change over time for unfolded proteins (Rogerson and Arteca, 2011). The Rg results for wild-type and three mutants showed a slight difference in their steadiness (Figures 4A-C), with similar average Rg value (around 15.2 Å) during the simulation time, suggesting relatively stable conformation. Root-mean-square deviation (RMSD) analysis gave insights into its structural conformation during the simulations, providing an indication of the stability of the protein and whether the simulation has equilibrated. The average backbone RMSD for WT and mutants were between 0.6 and 1.6 Å and remained stable through the entire MD simulations period (Figures 4D-F). Stable RMSD of the protein till the end of the simulation, suggested that the simulations were suitable for further rigorous analysis. In addition, the root-meansquare fluctuations (RMSF) were analyzed to figure out the impact of individual residues in both WT and MT xylanases (Figures 4G-I). There is slight residue fluctuation in regions adjacent to the mutation site. The increased variations in fluctuation in mutants indicated that point substitution (S31R, N29S, and I51V) involved in slight local structure change of xylanase. The three mutants and BaxA shared similar 3D model of GH11 xylanase, which resembled a β-jelly-roll fold or a partially closed right hand ( Figure S3).
In DS153, DS241, and DS428, the conformations of local residues showed variation, i.e., in the wild type, van der Waals interaction was observed between S31 and N32, but in DS241, S31 mutating to R31 vanished the interaction between R31 and N32, instead introduced the interaction between the guanidinium group of R31 and N29 (Figures S4A-C). We also plotted electrostatic potential colored van der Waals surface map and penetration graph of van der Waals surfaces. It could be clearly seen that due to formation of non-covalent interaction, there was inter-penetration between the van der Waals surfaces of the R31 and N29 (Figure S4D). The NCI caused the van der Waals surfaces of the molecules to penetrate each other. Between R31 and N29, non-bonded radius of NH2 R31 , OD1 N29 and the mutual penetration distance were 1.833, 1.705, and 0.453 Å, respectively. This was a non-trivial value, indicating the van der Waals force is obvious. However, in the I51V and N29S, there was no significant change in the environment around the mutated amino acid residues.
From the angle of protein structure, it was noticed that a residue substitution might have a dramatic effect on protein structure and function, especially in active site or cryptic site which occupied near to substrate-binding site (Worth et al., 2009). Distal single or several amino acid mutation/mutations may induce the change in molecular structure of enzyme, and the consequences could lead to opening or closing of the active site cleft (Shibuya et al., 2000;Zeymer and Hilvert, 2018). Hence, probing the impact of the substitution of distal residue on the active site architecture can provide valuable information for better understanding the altered protein structure and the function phenomenon.
The substitution of N29S (DS153) and S31R (DS241) occurred in A3 β-strand which was in the "finger" structural region, while the substitution of I51V (DS428) occupied in A5 β-strand which was on the "palm" structural region ( Figure S3). However, they could affect the architecture of the catalytic cleft based on the atomistic pair-wise distance analysis (Figures 5A-E). The NH1 atom of R112, OH atom of E78, E172, Y69, and Y80 were chosen as the representative atom. To elucidate a clear picture of the dynamic behavior along MD trajectory, the distribution of the atomic distance between R112-NH1 and E78-OH, E172-OH, Y69-OH, Y80-OH in WT and MT (DS241 and DS428) was calculated (Figures 5F-H). As a consequence of post-MD extracted protein structures, it was observed that the distance between both groups remained steady (from 3.5Å to 5.5Å) in the BaxA, showing the native conformation ("closed state") of wild-type enzyme. However, in S31R and I51V system, the interface distance was far (from 6 to 8 Å), in which the residue R112 and N35 were flipped open for favoring the protein conformation in an "open state" and ultimately lost the strong interface interactions.
The solvent-accessible surface area (SASA) analysis also revealed the same mechanism too. The SASA is a parameter which measures the fraction of the protein surface interacting with the solvent molecules and could be used for predicting the extent of the conformational change. The SASA was analyzed . Key residues were shown as a stick model. R112 and N35 belonged to the "gate" residues. E78, Y69, Y80, and E172 belonged to the active sites. The plots (F, BaxA; G, DS241; H, DS428) showed the dramatic change in distance along MD simulations time, and the frequency distribution of distances between residue R112 and other key residues (E78, E172, Y69, Y80, and N35). Protein structures were generated and viewed by PyMOL 2.2.0 (DeLano, 2002) and Chimera 1.13 (Pettersen et al., 2004).
using the cppptraj module in Amber Tools implemented linear combination of pairwise overlaps (LCPO) algorithm to explore how the distal substitution concealed or increased the accessibility of the catalytic site. Compared with BaxA, R112 slightly moved away from catalytic residues in DS241 and DS428, hence further boosting the SASA of mutants (Figure S5), indicating that once contacts of interface residues diminish or weaken, the SASA of the interface region will increase. Protein structure and dynamics may be beneficially changed by mutations that are distant from the active sites (Zeymer and Hilvert, 2018). These results delineate that the "open state" further promoted the catalytic efficiency to beechwood xylan.

NCI Between "Gate" Residues and Active Site Residues
Additionally, we have analyzed the NCI to elucidate the interface residues interaction between "gate" residues (R112 and N35) and the active site architecture residues (E78, E172, Y69, and Y80). Compared to the popular NCI method, such as RDG, the inter-fragment interactions could be studied individually. Therefore, mutual interference was avoided by IGM method. The value of δg function in the interaction region directly reflected the strength of the interaction (Lefebvre et al., 2017). The ellipsoids were intermolecular interaction, the green and blue part generally represented the dispersion force and electrostatic force, respectively. It was intuitive that there was interaction between the "gate" residues and the active site architecture residues in BaxA. There was a vast green oval region in R112/E78 and R112/Y80 (Figure 6A). But in DS241 and DS428, single site substitution (S31R and I51V) induced a dramatically rotameric conformational change of "gate" residues, in which most of the residues completely lost the atomic interactions (R112/E78 and R112/Y80) as shown in Figures 6B,C, and there was not revealing area of a green oval. Additionally, it was intuitive that the residue pair R112/N35 was still involved in van der Waals interactions in DS428.
Topological analysis of the interface residues has been carried out to study the NCI among the residue pairs. Table 3 showed that there were topological properties at the bond critical points (BCPs) including the parameters of electron density (ρ), Laplacian of electron density (∇ 2 ρ), kinetic energy density (G), and potential energy density (V) at the representative BCPs. FIGURE 6 | Intermolecular interactions between "gate" residues and active site residues. (A) BaxA; (B) DS241; (C) DS428. The blue represented a strong attraction, and red denoted a strong repulsion. All isosurfaces were colored according to a BGR (blue-green-red) scheme over the electron density range −0.05 a.u. < sign(λ 2 )ρ < 0.05 a.u. Bond critical points (BCPs, small orange spheres) using AIM analysis for different models (D) BaxA; (E) DS241; (F) DS428.
Several BCPs among R112/E78, R112/Y80, R112/N35 can be observed in Figures 6D-F. In Table 3, it was noticed that the residue R112 and E78 attempted to form NCI in BaxA. The ρ value between residue R112-NH1 and Y80-OH was 0.7256 a.u., but in DS241 and DS428, the ρ and ∇ 2 ρ-values were not in the range where weak interaction could be possible to form. Although the atoms of HH11, HD2, and HB2 in residue R112 and the OE1 in E78 formed a weak van der Waal's force in DS241, the value of ρ was minimal. In addition, in DS428, the HD3 atom of residue R112 formed van der Waals interaction with the OE1 of residues E78, where the value of ρ was 3.48E-06 a.u. It can be showed that this residue pair's interaction was still weaker, compared with the BCP intensity formed by residue R112 and E78 in BaxA. Therefore, the results of AIM further confirmed and suggested that the residue R112 in BaxA attempted to form an extensive interaction with residues in the catalytic cleft which did not exist in the DS241 and DS428.
To elucidate the impact of the substitution of N29S on the interaction of "gate" residues (W9 and P116), distribution of the residue pair's distance was calculated. For BaxA, W9 and P116 were in proximity to each other (close ensemble), remaining steady around 4 Å. However, the distance tremendously increased to ∼7 Å and dramatically reduced the NCI (open ensemble) in DS153 (Figures 7A,D,E). Change in SASA (Figure S5) also demonstrated substitution of N29S could inevitably cause the hydrophobic residues (W9 and P116) to be exposed to the aqueous surrounding. The 3D RDG isosurfaces with BGR color scales representing sign(λ 2 )ρ values were showed in Figures 7B,B ′ . It could be seen that the intercommunications between residue W9 and P116 were prevailed by van der Waals FIGURE 7 | Post-MD modeling presentation and frequency distribution of distances between "gate" residues in DS153. (A) the cartoon representation to illustrate the variation in "gate" residues (W9 and P116) in BaxA and DS153. (B,B ′ ) RDG isosurfaces for the WT and DS153 were individually plotted by using VMD 1.9.3. The green indicated van del Waals-type interactions. (C,C ′ ) color-mapped scatter map was used to facilitate identification of correspondence between spikes and RDG isosurfaces using gnuplot program. (D) The stick models represented the location of crucial "gate" residues. W9 and P116 were substrate-binding subsites and flipped open allowing the substrate access the catalytic cleft. (E) The frequency distribution of distances between P116 and W9 in DS153 and BaxA.
interaction (represented in green) in BaxA. However, in DS153, the interaction was completely lost. In order to quantitatively understand the interlayer interactions, the dependence of RDG on sign(λ 2 )ρ was also calculated with a color-scale from −0.03 to 0.02 a.u. as shown in Figures 7C,C ′ . Accordingly, it could be seen that many RDG spikes in scatter graph, were shown as a function of sign(λ 2 )ρ from negative to positive values, as well as the spikes corresponding to H-bonding interactions, van der Waals interactions and steric effects. Compared with BaxA, DS153 mutant lost the spike corresponding to W9-P116 ( Figures 7C,C ′ ). Possibly, GH11 xylanases had six sugar-binding subsites from −3 through +3 and preferentially bound more extended chain substrate (Brusa et al., 2018). The positive and negative numbers represented reducing and non-reducing end of the substrate, and the cleavage site was between the subsite −1 and +1, respectively (Wu et al., 2018). The residue R112 and N35 are the essential residues which attempt to form the −2 subsite (Brusa et al., 2018). And the N35 is the optimum pH-dependent residue of many GH11 xylanase. For the B. circulans xylanase (BCX), mutation of N35D shifted its pH optimum from 5.7 to 4.6, which further associated with increased activity (Joshi et al., 2000). The W9 and P116 were residues at subsite −2 and −1 of BaxA, respectively. The Crystal structure of BCX-xylotetraose complex revealed that the C5 of xylose residue 1 (xyl-1) bound with the CE2 and NE1 of W9 through hydrogen bond. The O3 atom of xylose residue 2 (xyl-2) interacted with NH2 atom of R112 and O of P116 with short hydrogen bond. Therefore, according to observations mentioned above, it could be delineated that the residue R112 played an essential role in the BCX, while mutation of N112 caused gain-ofloss (GOL) catalytic activity at a level of 65%, respectively (Wakarchuk et al., 1994).

DATA AVAILABILITY STATEMENT
All data generated or analyzed during this study are included in this published article and its Additional files.

AUTHOR CONTRIBUTIONS
M-QL designed the xylan hydrolysis experiments and supervised the study. J-YL is a collaborator in the MD simulations and structure analysis. AR helped in writing the manuscript and refining the graphs. XX performed the mutant library and enzymatic properties. Z-JG and R-CW performed the purification and inhibitory activity assay.