This is an attempt to understand the very nice argument by Peter Mueller (completed by OP, Sean Eberhard, in the comments). I am not completely happy with it, and would like to encourage others to look for more conceptual explanations.

Denote $\Omega:=P^1\mathbb{F}_q=\mathbb{F}_q\cup \infty$, so that $|\Omega|=q+1$ and the group $G=PSL(2,q)$ acts on $\Omega$ by projective transformations. Assume that $S\subset G=PSL(2,q)$, $|S|=q+1$ is chosen so that
$$\sum_{s\in S} \mathbb{1}_{s(i)=j}=1,\quad\forall i,j\in \Omega.$$
Then the same holds for the set $g_0S$ for any $g_0\in G$, thus we may suppose that $id\in S$, therefore other elements of $S$ do not have fixed points. Denote by $G^\star$ the set of fixed-point-free elements from 𝐺 together with 1, we have $S\subset G^\star$. Consider the function $M(i,j)$ on $\Omega\times \Omega$ defined as
$$
M(i,j)=\begin{cases}1,\, \text{if}\,\, i=\infty,\, j\ne \infty\\
1,\, \text{if}\,\, i,j\in \mathbb{F}_q,\chi(i-j)=1\\
0,\, \text{otherwise}.\end{cases}
$$
Here $\chi$ is a quadratic character of $\mathbb{F}_q$ (Legendre symbol if $q$ is prime).
We get
$$
\sum_{i,j\in \Omega,s\in S} M(i,j)\mathbb{1}_{s(i)=j}=\sum_{i,j\in \Omega} M(i,j)=q+q(q-1)/2
$$
is odd. Thus to get a contradiction it suffices to prove that
$$
\sum_{i,j\in \Omega} M(i,j)\mathbb{1}_{s(i)=j}
$$
is even for any fixed element $s\in G^\star$. For $s=id$ all summands are just zeroes. If $s$ does not have fixed point, than
$$
\sum_{i,j\in \Omega} M(i,j)\mathbb{1}_{s(i)=j}=1+\sum_{i\in \mathbb{F}_q} M(i,s(i))=1+\left|i\in \mathbb{F}_q:\chi(s(i)-i)=1\right|.
$$
Note that for fixed $\alpha\in \mathbb{F}_q$, the equation $$s(x)-x=\alpha\quad \quad (1)$$ is quadratic with respect to $x$, thus it has even number of roots unless it has a double root. We prove that there exists unique quadratic residue $\alpha$ for which (1) has a double root, the result would follow.

It may be proved by straightforward calculations, but let me give an argument if not conceptual, but at least almost calculations-free.

First of all, we use the following characterization of $PSL(2,q)$ (probably well known, but I did not see it before.)

**Lemma.** A projective transformation $s(x)=\frac{ax+b}{cx+d}, c\ne 0$, belongs to $PSL(2,q)$ if and only if the equation $s'(x)=1$ has two roots in $\mathbb{F}_q$.

**Proof.** Note that both properties are preserved when we replace $s(x)$ to $s(x)+C$ and $s(x+C)$ for constant $C\in \mathbb{F}_q$. So we may suppose that $a=0$ (subtract constant $a/c$ from $s(x)$) and also that $d=0$ (replace $s(x)$ by $s(x-d/c)$). So $s(x)=M/x$ for some constant $M$, and both conditions say that $-M$ is a non-zero square.

Next observation concerns the discriminants of the quadratic equation and critical values of corresponding functions. For the equation $f(x)=-x^2+2ax-b=0$, we all know that its discriminant $a^2-b$ equals $f(x_0)$, where $x_0=a$ is the critical value. But we deal with a slightly different equation $h(x):=s(x)-x=0$. After shifting the variable we get a function of the form $h_0(x)=-A/x+B-x$. Its discriminant is $D=B^2-4A$, and the critical values are $x_1=\sqrt{A},x_2=-\sqrt{A}$, so in our situation $A$ is a square by Lemma. Also we get $$D=h_0(x_1)h_0(x_2).$$

Since $h(x)=s(x)-x$ does not have roots, $D$ is not a square. Thus exactly one of two values of $h$ at critical points is a square. This value is the unique appropriate $\alpha$, as desired.

1more comment