z-Transform Methods: Definition vs. Integration Rule
Just to avoid a misunderstanding: the $\mathcal{Z}$-transform is a transform defined for sequences, comparable to the Laplace transform for continuous functions. What you are talking about is not the $\mathcal{Z}$-transform, but methods for converting analog to digital (actually, discrete-time) systems. [And it doesn't help that one of those conversion methods is called the matched Z-transform, which is just an unfortunately chosen name.]

There are several methods to achieve such a conversion, and all of them have certain properties. There is no ideal or "correct" mapping, and some methods may be more well-suited for certain applications than others.

The methods mentioned in your question are direct mapping methods, where the variable $s$ in the continuous-time transfer function is simply replaced by a (rational) function of $z$. Apart from the two methods mentioned in your question (1. Tustin's method - also called bilinear transform, and 2. forward Euler), there is also the backward Euler method, which replaces $s$ by $(1-z^{-1})/T$. Note that Tustin's method (bilinear transform) is equivalent to the trapezoidal integration rule. The important advantage of the bilinear transform for the conversion of analog systems is that the frequency axis ($j\omega$-axis) of the $s$-plane is transformed to the unit circle of the $z$-plane, and, consequently, the frequency response of the resulting discrete-time system is just a compressed (warped) version of the analog frequency response. One consequence of this is that analog frequency selective filters with optimal magnitudes (Chebysev, Butterworth, elliptic) are mapped to discrete-time filters satisfying the same optimality criteria. Furthermore, since the left half-plane is mapped to the inside of the unit circle in the $z$-plane, stability (and the minimum-phase property) are preserved by the transformation.

This is different for the forward and backward Euler methods. With those methods, the $j\omega$-axis is not transformed to the unit circle of the $z$-plane, and the left half-plane is not mapped to the region $|z|&lt;1$. In case of the forward Euler method, the $j\omega$-axis is mapped to the line $\textrm{Re}\{z\}=1$, and the left half-plane is mapped to the half-plane $\textrm{Re}\{z\}&lt;1$. This means that a stable analog prototype filter might be transformed to an unstable discrete-time filter. The backward Euler method transforms the $j\omega$-axis to the circle $|z-\frac12|=\frac12$, and the left half-plane is mapped to the region inside that circle: $|z-\frac12|&lt;\frac12$. Consequently, stability and the minimum phase property are preserved, but the analog frequency axis is not mapped to the unit circle of the complex plane, so the resulting frequency response is not just a warped version of the original analog frequency response.

Note that there are also other very common methods for converting analog to digital systems, which are not based on direct mappings from the $s$-plane to the $z$-plane. These methods are based on time domain criteria. The two most common of these methods are the impulse-invariant method and the step-invariant method. They preserve the shape of the impulse (or step) response of the analog system by defining the desired digital system as the one whose impulse (step) response is a sampled version of the corresponding analog impulse (step) response. These methods introduce aliasing in the frequency domain, but - unlike the bilinear transform - they do not introduce frequency warping.