If $M,N$ are two Riemann surfaces with boundary, then we can glue them along one of each of their boundary component, which is $S$, to form a new Riemann surface with boundary, but for different gluing we may form different Riemann surfaces with boundary, for example, there may be a $S$ twist, intuitively, it is just we rotate one $S$ an angle then glue it with another surfaces, but my question is how can we show the resulting two surfaces (twisted gluing and untwisted gluing) are different Riemann surface with boundary (their differential structures are the same because we can regard it as a kind of connected sum)? Thanks!

Here is a particularly simple example. We will take a cylinder and glue its ends together in two different ways. The cylinder will be $$C:=\{ z: 0 \leq \mathrm{Im}(z) \leq 1 \} / \mathbb{Z},$$ where $\mathbb{Z}$ acts on the horizontal strip by translation. Let $\theta$ be a real number between $0$ and $1$. We will glue the ends of $C$ together by identifying $x+\mathbb{Z}$ with $(x+i+\theta) + \mathbb{Z}$, for $x$ real. Call the resulting torus $T_{\theta}$. We claim that all the $T_{\theta}$ are nonisomorphic as Riemann surfaces. Here is a sketch of a proof: Let $\Lambda_{\theta}$ be the lattice in $\mathbb{C}$ generated by $1$ and $i+\theta$. Clearly, $T_{\theta} \cong \mathbb{C} / \Lambda_{\theta}$. Suppose we had holomorphic isomorphism $f:T_{\theta} \to T_{\theta'}$. Then $f$ would lift to a map on the universal covers, $g: \mathbb{C} \to \mathbb{C}$, and there would be some additive map $h: \Lambda_{\theta} \to \Lambda_{\theta'}$ such that $$g(z+\lambda) = g(z) + h(\lambda)\ \mbox{for} \ \lambda \in \Lambda_{\theta}. \quad (*)$$ Let $D$ be a fundamental domain for $\Lambda_{\theta}$. The function $g$ is $O(1)$ on $D$ (compactness), so, applying $(*)$, we can deduce that $g(z) = O(\max(z, 1)$. Since $g$ is analytic, this implies that $g$ is of the form $z \mapsto az+b$. In particular, multiplication by $a$ must carry $\Lambda_{\theta}$ to $\Lambda_{\theta'}$. But there is no complex number $a$ such that $a \Lambda_{\theta} = \Lambda_{\theta'}$, a contradiction. 


If your compact Riemann surfaces $M$ and $N$ have a hyperbolic metric in which the boundary curves are totally geodesic of the same length, then this follows from the FenchelNielsen coordinate parametrization of Teichmuller space. Any compact Riemann surface with boundary has a uniformization to a hyperbolic surface with totally geodesic boundary (one may see this by doubling, and applying the uniformization theorem). If you perform this uniformization, and then compare the two sides, your identification between boundary components may differ from the parametrization given by FenchelNielsen coordinates (for example, the lengths of the boundary components may differ), and so it's more complicated to see how the conformal structure changes. In fact, if $M$ or $N$ is a disk or annulus, then the conformal structure may not be changing at all. If $M$ and $N$ both have Euler characteristic less than zero, then large twists should change the conformal structure. For example, a twist by $2\pi$ will change the Riemann surface by Dehn twist, and therefore gives a different conformal structure (up to isotopy). There's also a slight issue of moduli, in that a rotation by $2\pi$ changes the point in Teichmuller space, but not in moduli space. I'm implicitly assuming you're asking for the Teichmuller parameter. If you're asking for the parameter in moduli space, then the twist will take you along a nontrivial closed loop in moduli space, so at least there will be uncountably many different conformal moduli as you twist. 


It shall amount to understand precisely the parametrization of Teichmüller space. The key word is "pant decomposition". If you are willing to admit the (length, twist) parameters of Teichmüller space, then simply decompose the two initial surfaces in pants. The union of the two decomposition is a pant decomposition of both glued surfaces, with different parameters, so that they must be different surfaces. 


When you glue $M$ and $N$, for example, along some boundary component $C$. Choose a closed geodesic $\gamma$ on $S$ intersect with $C$. It can be shown that the hyperbolic length of $\gamma$ is changed when you twist along $C$. This means that different twist makes different hyperbolic surface. 

