Ok, I figured out the answer to my own question, though without the $\ln(ep K^2/\sigma^2)$ in the denominator. I'll see whether I can get that later (or maybe someone else sees it?). We can assume without loss of generality that $\mathbb{E} X_i = 0$ for each $i$ (otherwise apply what follows with the random variables replaced by $X_i - \mathbb{E} X_i$). The trick is to use that $\|\langle g,x\rangle\|_p$$\|\langle g,x \rangle\|_p$ is not justmonotonically increasing as a function of $O(\sqrt{p}\|x\|_2)$ but is$\|x\|_2$ $\Theta(\sqrt{p}\|x\|_2)$(this follows by 2-stability of the gaussian). Let $\sigma_i$ be i.i.d. Rademacher and $g_i$ be i.i.d. normal. \begin{align} \|\sum_i X_i\|_p &= \|\sum_i X_i - \mathbb{E} X_i\|_p\\ &\le \|\sum_i (X_i - X_i')\|_p\\ &\le 2\|\sum_i \sigma_i X_i\|_p\\ &= 2\sqrt{\frac{\pi}{2}}\|\mathbb{E}_g\sum_i \sigma_i |g_i| X_i\|_p\\ &\lesssim \|\sum_i \sigma_i |g_i| X_i\|_p\\ &= \|\sum_i g_i X_i\|_p (1)\\ &\lesssim \sqrt{p} \|(\sum_i X_i^2)^{1/2}\|_p\\ &\le \sqrt{p} \|\sum_i X_i^2\|_p^{1/2}\\ &= \sqrt{p} \|\sum_i X_i^2 - \mathbb{E}\sum_i X_i^2 + \sum_i X_i^2\|_p^{1/2}\\ &\le \sigma\sqrt{p} + \sqrt{p}\|\sum_i X_i^2 - \mathbb{E}\sum_i X_i^2\|_p^{1/2}\\ &\lesssim \sigma\sqrt{p} + \sqrt{p}\|\sum_i g_i X_i^2\|_p^{1/2}\\ &= \sigma\sqrt{p} + \sqrt{p}\left(\mathbb{E}_X\|X\|_\infty^p\left(\mathbb{E}_g|\sum_i g_i X_i \frac{X_i}{\|X\|_\infty}|^p\right)\right)^{1/2p}\\ &\le \sigma\sqrt{p} + \sqrt{pK}\left(\mathbb{E}_X\mathbb{E}_g|\sum_i g_i X_i \frac{X_i}{\|X\|_\infty}|^p\right)^{1/2p}\\ &\lesssim \sigma\sqrt{p} + \sqrt{pK}\left(\mathbb{E}_X\mathbb{E}_g|\sum_i g_i X_i|^p\right)^{1/2p} (2)\\ &= \sigma\sqrt{p} + \sqrt{pK}\|\sum_i g_i X_i\|_p^{1/2} (3) \end{align}
Inequality (2) used that $\|\langle g,x\rangle\|_p$ is $\Theta(\sqrt{p}\|x\|_2)$ for anymonotonically increasing as a function of $x$$\|x\|_2$, as mentioned above. Now notice we have bounded (1) in terms of its square root in (3). Thus $\|\sum_i g_i X_i\|_p^{1/2}$ is at most the larger root of the associated quadratic equation, and this gives us the bound $$ \|\sum_i g_i X_i\|_p \lesssim \sigma\sqrt{p} + Kp $$ as desired.
I tried to be careful and write all the details to make it clear, but much of the above argument can be compressed once you're familiar with the right bag of tricks (symmetrization, gaussian comparison principle, and triangle inequality in the right places).