The evolution operator of a general state $\psi(S,E)\in \mathscr{H}_S\otimes \mathscr{H}_E$ is indeed linear, but it involves both the variables of $\mathscr{H}_S$ (that I denoted with $S$), and of $\mathscr{H}_E$ (denoted with $E$). Thus, given a self-adjoint operator $H(S,E)$ on $\mathscr{H}_S\otimes \mathscr{H}_E$, you can associate the linear unitary evolution group $e^{-it H(S,E)}$, and $$\rho_{S+E}(t)= e^{-it H(S,E)}\rho_{S+E}e^{it H(S,E)} \; .$$ However in general $H(S,E)$ is not a tensor product of two operators acting on $\mathscr{H}_S$ and $\mathscr{H}_E$. Suppose that it has the folowing form: $H(S,E)=H(S)\otimes 1 + 1\otimes H(E)$ (the sum of two parts that act only on one space). Then you have
$$\rho_{S+E}(t)=e^{-itH(S)}e^{-itH(E)}\rho_{S+E}e^{it H(E)}e^{it H(S)}\; ;$$
and you can take the partial trace and obtain
$$\rho_{S}^E(t)=e^{-itH(S)}\rho^E_{S+E}(t)e^{it H(S)}\; ;$$
where $\rho^E_{S+E}(t)=\mathrm{Tr}^E[e^{-itH(E)}\rho_{S+E}e^{it H(E)}]$. You see that you do not have unitary evolution on $\mathscr{H}_S$, because $\rho^E_{S+E}(t)$ depends on time (the partial trace is not one, is a function of the $S$ variables and is in general modified by the time evolution). However if $\rho_{S+E}=\rho_S\otimes\rho_E$, everything factors and you get
$$\rho_{S+E}(t) = e^{-itH(S)}\rho_S e^{itH(S)}e^{-itH(E)}\rho_E e^{itH(E)}$$
and taking the trace on $E$, since your state is normalized to have trace 1 (i.e. both $\rho_S$ and $\rho_E$ have trace one), you get the evolution $\rho^E_{S}(t) = e^{-itH(S)}\rho_S e^{itH(S)}$.
If the system interacts with the environment, i.e. the Hamiltonian is not of the form above, or if the state is not factorized, the time evolution "mixes" between the two spaces, and you can't get a linear evolution in just $\mathscr{H}_S$.
EDIT: Maybe it is better to explain what is meant by linearity, because it may not be clear. An operator $A$ on a Hilbert space $\mathscr{H}$ is linear if $\forall \psi,\phi\in \mathscr{H}$ and $\lambda\in \mathbb{R}$:
$$A(\lambda(\psi+\phi))=\lambda (A(\psi)+A(\phi))\; .$$
Now the solution of
$$i\partial_t\rho=[H,\rho]$$
is a linear operator acting on the initial state, because both hands of the equation are linear in $\rho$. Now derive in time the expression of $\rho^E_S(t)$ above (the partial trace of your time evolved general state):
$$i\partial_t\rho^E_S = [H(S),\rho^E_{S}]+ e^{-it H(S)}(i\partial_t \rho^E_{S+E}(t))e^{itH(S)}\; .$$
The second term on the right hands side is in general not linear in $\rho^E_S$, even if it may be in some particular case (for sure factorized states $\rho_S\otimes\rho_E$, it may be difficult to give characterization of other ones). This example may be another possibility: let $\psi(x)\in \mathscr{H}_x$ and $\phi(y)\in\mathscr{H}_y$ and define $\Psi(x,y)=\psi(x)\phi(y)e^{ixy}$. This state is almost factorized. Let $H(E)=-\Delta_y$ and $\rho_{S+E}=\lvert\Psi(x,y)\rangle\langle \Psi(x,y)\rvert$. Then the last equation becomes
$$i\partial_t\rho^E_S = [H(x),\rho^E_{S}]+ e^{-it H(x)}[x^2,\rho_{S+E}^E(t)]e^{itH(x)}-e^{-it H(x)}(\mathrm{Tr}^y[e^{it\Delta_y}(\lvert \psi(x)e^{ixy}(\Delta_y \phi(y))\rangle\langle \Psi(x,y)\rvert + \lvert \Psi(x,y)\rangle\langle \psi(x)e^{ixy}(\Delta_y \phi(y))\rvert)e^{it\Delta_y}])e^{itH(x)}\; .$$
Suppose now that $\Delta_y\phi_y=0$, then we get
$$i\partial_t\rho^E_S = [H(x),\rho^E_{S}]+ e^{-it H(x)}[x^2,\rho_{S+E}^E(t)]e^{itH(x)}\; .$$
Suppose also that the commutator $e^{-itH(x)}x^2e^{it H(x)}=F(t,x)$, for some operator $F(t,x)$. Then
$$i\partial_t\rho^E_S = [H(x),\rho^E_{S}]+ F(t,x)\rho^E_S - \rho^E_S(F(t,x))^\dagger\; .$$
This would eventually lead to a non-unitary evolution with time-dependent generator, but linear ;-). Sometimes the last equation is called master equation.
