Here are just a few observations to discard some trivial cases. I'll talk directly about the Istratescu formulation and will prefer the probabilistic language, but that shouldn't be a problem after everything that has been posted already.
First of all, the problem is equivalent to asking what is the best constant in the inequality $E[|X-X'|^p]^{1/p}\le C_{p,q}E[|X|^q]^{1/q}$ where $X$ is an arbitrary real random variable and $X'$ is an independent copy of $X$.
Indeed, our constant cannot be any better because we can just take a sequence of independent copies of $X$ normalized by $E[|X|^q]=1$. On the other hand, suppose that $X_j$ is any sequence of random variables satisfying $E[|X_j|^q]\le 1$. Take large integer $N$ and consider the random variable $X$ taking real values $x_j$ ($j=1,\dots,N$) with probability $1/N$ each.
Then we conclude that
$$
\left[\frac 1{N^2}\sum_{1\le i,j\le N}|x_i-x_j|^p\right]^{q/p}
\le C_{p,q}^q\frac 1N\sum_{1\le j\le N}|x_j|^q\,.
$$
In particular, for each fixed $\omega$ in the probability space, we have
$$
\left[\frac 1{N^2}\sum_{1\le i,j\le N}|X_i(\omega)-X_j(\omega)|^p\right]^{q/p}
\le C_{p,q}^q\frac 1N\sum_{1\le j\le N}|X_j(\omega)|^q
$$
Taking the expectations of both sides, we conclude that
$$
\left[\frac 1{N^2}\sum_{1\le i,j\le N, i\ne j}E|X_i(\omega)-X_j(\omega)|^p\right]^{q/p}\\ \le
E\left[\frac 1{N^2}\sum_{1\le i,j\le N, i\ne j}|X_i(\omega)-X_j(\omega)|^p\right]^{q/p}
\le C_{p,q}^q\,,
$$
so
$$
\left[\frac {N-1}{N}\min_{1\le i,j\le N, i\ne j}E|X_i(\omega)-X_j(\omega)|^p\right]^{1/p}
\le C_{p,q}\,.
$$
This shows that $C_{p,q}$ can be used in the original problem (by Fedor Petrov's remark about infinite Ramsey, say).
The most trivial case now is $p=2$. It boils down to the observation that
$$
E[|X-X'|^2]^{1/2}=\sqrt 2(E[|X|^2]-E[X]^2)^{1/2}\le\sqrt 2E[|X|^2]^{1/2}\le \sqrt 2E[|X|^q]^{1/q}
$$
for all $q\ge 2$ and independent Rademachers give the identity.
This can be easily generalized to $p\ge 2$. All we need is the following elementary inequality:
$$
|x-y|^p\le 2^{p-2}\left||x|^{p/2}\operatorname{sgn} x-|y|^{p/2}\operatorname{sgn} y\right|^2\,,\quad x,y\in\mathbb R\,,
$$
which reduces to 2 inequalities
$$
(1-t)^{p/2}\le 2^{\frac p2-1}(1-t^{p/2})\quad\text{and}\quad
(1+t)^{p/2}\le 2^{\frac p2-1}(1+t^{p/2})
$$
for $t\in[0,1]$, the first of which holds because the LHS is convex, the RHS is concave, and the endpoints $0$ and $1$ are fine, while the second one is just the convexity of $s\mapsto s^{p/2}$.
Once we have it, we just write
$$
E[|X-X'|^p]^{1/p}\le
2^{1-\frac 2p}E\left[\left||X|^{p/2}\operatorname{sgn}X-|X'|^{p/2}\operatorname{sgn}X'\right|^2\right]^{1/p}=
\\
2^{1-\frac 1p}(E[|X|^p]-E[|X|^{p/2}\operatorname{sgn}X]^2)^{1/p}\le 2^{1-\frac 1p}E[|X|^p]^{1/p}\le 2^{1-\frac 1p}E[|X|^q]^{1/q}\,,
$$
and the Rademachers give an equality again, so
$$
C_{p,q}=2^{1-\frac 1p}
$$
for any $q\ge p\ge 2$.
The interesting case is thus $p\in[1,2)$. I will do $p=1$ now, which can be pulled through to a final closed formula.
The heuristic is simple. Let $\lambda$ be the density of $X$. Then we want to maximize $\iint |x-y|\lambda(x)\lambda(y)\,dx\,dy$ subject to the constraints $\lambda\ge 0$, $\int\lambda(x)\,dx=1$, $\int|x|^q\lambda(x)\,dx=\operatorname{const}$. It is tempting to use the Lagrange multipliers technique and to write
$$
I(x)=\int |x-y|\lambda(y)\,dy=\alpha+\beta|x|^q\,,
$$
from where $\lambda(x)\asymp |x|^{q-2}$ (the second derivative of the LHS is $2\lambda(x)$).
This doesn't seem to make much sense until you realize that the Lagrange equation doesn't really need to hold on the entire $\mathbb R$, just on the support of $\lambda$. This prompts one to consider the probability densities
$$
\lambda_M(x)=\frac{(q-1)|x|^{q-2}}{2M^{q-1}}\chi_{[-M,M]}\,.
$$
The corresponding value of the potential $I_M(x)$ can be determined from
$$
I_M(0)=\int|y|\lambda_M(y)\,dy=\frac{q-1}q M\,,
$$
so $\alpha=\frac{q-1}{q}M$,
and $I_M''(x)=2\lambda_M(x)$, so $\beta=\frac {1}{q}M^{1-q}$.
Note also that for $|x|>M$, $I_M$ is just linear with correct slope at $\pm M$, so we have $I_M(x)\le \alpha+\beta|x|^q$ there.
Now we have
$$
\int|x|^q\lambda_M(x)\,dx=\frac{q-1}{2q-1}M^q
$$
and
$$
\iint |x-y|\lambda_M(x)\lambda_M(y)\,dx\,dy=\int I_M(x)\lambda_M(x)\,dx
\\
=M\left[\frac{q-1}q+\frac1q\frac{q-1}{2q-1}\right]=\frac{2(q-1)}{2q-1}M\,.
$$
Thus,
$$
C_{1,q}\ge 2\left(\frac{q-1}{2q-1}\right)^{1-\frac 1q}\,.
$$
Now we want to show that $\lambda_M$ is, indeed, optimal. Let $\lambda$ be any compactly supported probability density (I'll skip the mumbo-jumbo about how to pass to the limit to get arbitrary probability measures from here). Choose $M$ so that $\int|x|^q\lambda(x)\,dx=\int|x|^q\lambda_M(x)\,dx$. Then we can write
$$
\iint|x-y|\lambda(x)\lambda_M(y)\,dx\,dy=\int I_M(x)\lambda(x)\,dx\le
\int(\alpha+\beta|x|^q)\lambda(x)\,dx
\\
=\int(\alpha+\beta|x|^q)\lambda_M(x)\,dx=\int I_M(x)\lambda_M(x)\,dx=\iint|x-y|\lambda_M(x)\lambda_M(y)\,dx\,dy\,.
$$
On the other hand, $|x|$ is negative definite in the sense that for every zero integral compactly supported real function $\psi$,
$$
\iint|x-y|\psi(x)\psi(y)\,dx\,dy\le 0\,.
$$
Applying it with $\psi=\lambda-\lambda_M$, we conclude that
$$
\iint |x-y|\lambda_M(x)\lambda_M(y)\,dx\,dy+\iint |x-y|\lambda(x)\lambda(y)\,dx\,dy
\\
\le 2\iint |x-y|\lambda(x)\lambda_M(y)\,dx\,dy\le 2\iint |x-y|\lambda_M(x)\lambda_M(y)\,dx\,dy\,,
$$
so
$$
\iint |x-y|\lambda(x)\lambda(y)\,dx\,dy\le \iint |x-y|\lambda_M(x)\lambda_M(y)\,dx\,dy\,,
$$
and we are done.
What happens for $1<p<2$? One can try to run the same argument, but to carry it out one would need to show that we can find a probability density $\lambda_1$ on $[-1,1]$ (going to $M$ is just trivial scaling) such that
$$
I(x)=\int_{[-1,1]}|x-y|^p\lambda_1(y)\,dy=\alpha+\beta|x|^q
$$
on $[-1,1]$ and $I(x)\le \alpha+\beta|x|^q$ for $|x|>1$. Then this density will be the extremizer ($|x|^p$ is still negative definite, so the rest of the argument will sail through). However, currently I cannot even prove the existence of $\lambda_1$, forget about a decent expression for it and related quantities. Any ideas?