The tl;dr for your particular example is that your symmetry assumption turns out to not help you at all in this setting (interestingly, this non-result is a bit of a knife-edge case: it would have helped you if you were interested instead in the mean value of $aE[X]+bE[Y]$ for any $a\neq b$, just not when $a=b$).
To answer questions like this more generally though, you might find the literature on semi-parametric efficiency bounds to be illuminating. In particular, Newey (1990) is a classic reference. There are a number of technical definitions to formalize all of the concepts below, and you can find the excruciating details in the linked paper and the references contained therein, but let me try to sketch some intuition for what can be said.
Consider some parameter of interest $\theta$, which is a functional of $F$, i.e. $\theta = \mu(F)$ where $\mu(\cdot)$ maps from some space of distributions to the parameter of interest (so in your case, $\theta = E[X+Y]$. We would like to be able to find a lower bound The basic idea here is that if $X,Y$ were to belong to a parametric distribution, parameterized by a parameter $\eta$, the Cramer-Rao lower bound would imply that provided $\mu$ is a sufficiently well behaved mapping, $\hat\theta = \mu(F(\cdot;\hat\eta_{MLE}))$ is the most efficient possible estimator of $\theta$. Then (subject to some more regularity conditions on the class of estimators under consideration) the lowest possible variance one could possibly hope to achieve for an estimator is the $\sup$ among the Cramer-Rao lower bounds of all parametric models containing the true distribution.
For parameters which are based on expected values or conditional expectations, the computation of this supremum is typically relatively simple. The details are somewhat difficult to state, but if you work through similar arguments as Newey, you end up finding that the following intuitive estimator turns out to be the most efficient:
$$\hat\mu_{X} = \hat\mu_Y = \frac1{2n}\sum_{i=1}^n X+Y$$
Then if you are interested in estimating the parameter $\theta = E[aX+bY]$, you would take
$$\hat\theta = a\hat\mu_X + b\hat\mu_Y = (a+b)\hat\mu_X$$
which has variance
$$n\mathrm{Var}(\hat\theta) = \frac{(a+b)^2}{2} \mathrm{Var}(X)$$
By contrast, when you do not use the information that $X=Y$, the best you can do is
$$\tilde\theta = a\underbrace{\frac1n\sum_{i=1}^n X_i}_{\equiv \tilde\mu_X} + b\underbrace{\frac1n\sum_{i=1}^n Y_i}_{\equiv \tilde\mu_Y}$$
In this case, you would have
$$n\mathrm{Var}(\tilde\theta) = a^2 \mathrm{Var}(X) + b^2\mathrm{Var}(Y) = (a^2+b^2)\mathrm{Var}(X)$$
Since $x\mapsto x^2$ is a strictly convex function, we have that
$$\frac{a^2 + b^2}{2} \geq \left(\frac{a+b}{2}\right)^2 \implies \frac{(a+b)^2}2 \leq a^2+b^2$$
with strict inequality when $a\neq b$. Then in general, the information that $X,Y\overset{i.i.d.}\sim F$ leads to a lower variance estimator of $\theta = E[aX+bY]$, but as chance would have it, this information is not helpful in your particular case where $a=b=1$.
