The assumption behind the Saint Venant-Kirchhoff constitutive model is that the energy per unit volume of the material in the reference configuration (i.e. strain energy density) can be adequately approximated as a quadratic function of the Lagrangian Green strain $E$. (The logic is similar to the approximation of a point particle's potential energy being quadratic near a stable equilibrium, but is ultimately guided by experimental validation.)
Therefore, the generic form of the strain energy density in the Saint Venant-Kirchhoff model, in index notation, is of the following form:
$$W = \frac{1}{2}C_{ijkl}E_{ij}E_{kl}$$
where $C = C_{ijkl}$ is a fourth-order tensor of elastic constants.
Note that, due to the symmetry of the Lagrangian Green strain tensor and the fact that $W$ is a scalar, $C$ also possesses quite a few internal symmetries as well; $C_{ijkl} = C_{jikl} = C_{ijlk} = C_{lkij}$.
If we take the partial derivative of this general expression with respect to the Lagrangian Green strain tensor, we find:
$$\frac{\partial W}{\partial E_{mn}} = \frac{1}{2}C_{ijkl}\frac{\partial E_{ij}}{\partial E_{mn}}E_{kl} + \frac{1}{2}C_{ijkl}E_{ij}\frac{\partial E_{kl}}{\partial E_{mn}}$$
This results in a few Kronecker deltas:
$$\frac{\partial W}{\partial E_{mn}} = \frac{1}{2}C_{ijkl}\delta_{im}\delta_{jn}E_{kl} + \frac{1}{2}C_{ijkl}E_{ij}\delta_{km}\delta_{ln}$$
$$\frac{\partial W}{\partial E_{mn}} = \frac{1}{2}C_{mnkl}E_{kl} + \frac{1}{2}C_{ijmn}E_{ij}$$
Since the Lagrangian green strain tensor is symmetric ($E_{ij} = E_{ji}$) and $C$ possesses the symmetries listed above, the two terms combine to make:
$$\frac{\partial W}{\partial E_{mn}} = S_{mn} = C_{mnkl}E_{kl}$$
Or, in index-free notation,
$$S = C:E$$
Now suppose that you assume the material is also isotropic; that is, $C$ is an isotropic tensor. As it turns out, the most generic form for a fourth-order isotropic tensor $C$ would be:
$$C_{ijkl} = \lambda\delta_{ij}\delta_{kl} + \nu_1 \delta_{ik}\delta_{jl} + \nu_2 \delta_{il}\delta_{jk}$$
where $\lambda$, $\nu_1$ and $\nu_2$ are unknown parameters.
Inputting this into the expression for $S$ above, noting again that $E_{ij} = E_{ji}$, and defining $\mu = \nu_1 + \nu_2$, we find that:
$$S_{ij} = \lambda (E_{kk}) \delta_{ij} + \nu_1 E_{ij} + \nu_2 E_{ji}$$
$$S_{ij} = \lambda (E_{kk}) \delta_{ij} + 2 \mu E_{ij}$$
Or, in index-free notation,
$$S = \lambda\ \text{tr}E\ I + 2 \mu E$$
Instead of making the assertion that $C$ is isotropic after deriving a general expression for $S$, you can simply take the statement that the Saint Venant-Kirchhoff model for the strain energy density is accurate based on empirical modeling and then derive an expression for $S$ from that. If so, consider the strain-energy density function for the Saint Venant-Kirchhoff model in index notation, which you can either take as a postulate or can be derived from the general expression for $W$ above but substituting in the isotropic form of $C$:
$$W(E_{ij}) = \frac{\lambda}{2} (E_{kk})^2 + \mu (E_{lm}E_{ml})$$
The second Piola-Kirchhoff stress tensor $S_{ij}$ is equal to the partial derivative of $W$ with respect to $E_{ij}$:
$$S_{ij} = \frac{\partial W}{\partial E_{ij}}$$
Taking the partial derivative, we find the following, very messy, expression:
$$S_{ij} = \lambda (E_{kk}) \frac{\partial E_{kk}}{\partial E_{ij}} + \mu (E_{lm}\frac{\partial E_{ml}}{\partial E_{ij}}) + \mu (\frac{\partial E_{ml}}{\partial E_{ij}}E_{lm})$$
This actually has everything we need; we just need to flip back to an index-free notation. The first identity to spot here is that the derivative of the trace of a matrix with respect to each matrix element is simply the identity matrix:
$$\frac{\partial E_{kk}}{\partial E_{ij}} = \delta_{ij}$$
You will also find that similar identity matrices pop up in the derivatives of each matrix element with respect to themselves:
$$\frac{\partial E_{ml}}{\partial E_{ij}} = \delta_{mi}\delta_{lj}$$
$$\frac{\partial E_{lm}}{\partial E_{ij}} = \delta_{li}\delta_{mj}$$
Inserting all of this back in, you'll find the following, far more simplified, expression:
$$S_{ij} = \lambda (E_{kk}) \delta_{ij} + \mu (E_{lm}\delta_{li}\delta_{mj}) + \mu (\delta_{mi}\delta_{lj} E_{lm})$$
Using index substitution on the last two terms:
$$S_{ij} = \lambda (E_{kk}) \delta_{ij} + \mu (E_{ij} + E_{ji})$$
The key remaining insight is to observe that the Lagrangian Green strain tensor is symmetric: $E_{ij} = E_{ji}$. Therefore,
$$S_{ij} = \lambda (E_{kk}) \delta_{ij} + 2 \mu (E_{ij})$$
which, in index-free notation, is the desired answer:
$$S = \lambda\ \text{tr}E\ I + 2 \mu E$$