Here is an attempt to answer my own question. Please feel free to comment and point out mistakes if you see some.
(Edit: I am learning this stuff right now, so I would really appreciate if an expert could confirm or refute my arguments.)
Lemma: Let $S,T$ be normal operator on a finite dimensional inner-product space, and suppose that $ST=TS$. Then, there is an orthonormal basis for which both $S$ and $T$ are diagonal.
Proof: Since $S$ is normal, $S$ is diagonalizable, so its minimal polynomial factors into distinct linear terms and thus, by the Primary Decomposition Theorem,
$$V=E_{\lambda_1}\oplus \cdots\oplus E_{\lambda_k}$$
where $\lambda_i$ are the distinct eigenvalues of $S$ and $E_{\lambda_i}=\ker(S-\lambda_i I)$. Now, if $v\in E_{\lambda_i}$, then
$$(S-\lambda_i I)T(v)=ST(v)-\lambda_i T(v)=TS(v)-\lambda_i T(v)=\lambda_i T(v)-\lambda_i T(v)=0$$
so $E_{\lambda_k}$ is $T$-invariant, and thus $T_i:=T|_{E_{\lambda_i}}$ is diagonalizable. Hence,
$$E_{\lambda_i}=E_{\mu^i_1}\oplus\cdots\oplus E_{\mu^i_{n_i}}$$
where $\mu^i_1,\ldots,\mu^i_{n_i}$ are the distinct eigenvalues of $T_i$, and $E_{\mu^i_j}=\ker(T_i-\mu^i_j I)$.
Let $B^i_j$ be an orthonormal basis for $E_{\mu^i_j}$. Since $T$ is normal and, for fixed $i$, $\mu^i_1,\ldots,\mu^i_{n_i}$ are distinct, we have $E_{\mu^i_1}\bot\cdots\bot E_{\mu^i_{n_i}}$, so $B^i_1\cup\cdots\cup B^i_{n_i}$ is orthonormal. Moreover, since $S$ is normal and $\lambda_1,\ldots,\lambda_k$ are distinct, $B:=\cup B^i_j$ is still orthonormal. Now, $B$ is an orthonormal basis of eigenvectors for both $S$ and $T$.
Q.E.D.
Now the result is clear. We have
$$[ST^*]_B=[S]_B[T]^*_B=[T]^*_B[S]_B=[T^*S]_B$$
so $ST^*=T^*S$, and similarly, $S^*T=TS^*$.
Thus,
$$(S+T)(S+T)^*=SS^*+TS^*+ST^*+TT^*=S^*S+S^*T+T^*S+T^*T=(S+T)^*(S+T)$$
and
$$(ST)(ST)^*=STT^*S^*=T^*S^*ST=(ST)^*(ST)$$