Are all zeros of $\Gamma(s) \pm \Gamma(1-s)$ on a line with real part = $\frac12$ ?
This is a continuation of the argument above, which completes the argument.
Let $C_n$ denote the square with vertices $[n \pm 1/2, \pm 4 I]$ for a positive integer $n$. We have the following inequalities for $z \in C_n$ and $n \ge 15$: $$|\sin(\pi z)| \ge 1, \quad z \in C_n.$$ $$|\Gamma(z)| \ge \frac{1}{2} \Gamma(n - 1/2),$$ $$|\Gamma(1-z)| \le \frac{\pi}{\Gamma(n - 1/2)} \le 1,$$ $$|\psi(1-z)|, |\psi(z)| \le 2 \log(n), $$
The first is easy, the second follows from Stirling's formula (this requires $n$ to be big enough, and also requires $z$ to have imaginary part at most $4$), the third follows from the previous two by the reflection formula for $\Gamma(z)$, the last follows by induction and by the formula $\psi(z+1) = \psi(z) + 1/z$. It follows that $$\left| \frac{1}{2 \pi i} \oint_{C_n} \frac{\Gamma'(z)}{\Gamma(z)} - \frac{d/dz (\Gamma(z) + \theta \cdot \Gamma(1-z))}{\Gamma(z) + \theta\cdot \Gamma(1-z)} \right|$$ $$= \left| \frac{1}{2 \pi i} \oint_{C_n} \frac{\theta \Gamma(1-z) (\psi(1-z) + \psi(z))} {\Gamma(z) + \theta \cdot \Gamma(1-z)} \right|$$ $$ \le \frac{8 |\theta| \cdot \log(n) \pi}{2 \pi \cdot \Gamma(n - 1/2)} \oint_{C_n} \frac{1} {|\Gamma(z) + \theta \cdot \Gamma(1-z)|}$$ $$ \le \frac{8 |\theta| \cdot \log(n) \pi}{2 \pi \cdot \Gamma(n - 1/2)} \cdot \frac{1}{1/2 \Gamma(n - 1/2) + 1} \ll 1,$$ where $\theta = \pm 1$ (or anything small) and $n \ge 15$, where the final inequality holds by a huuuge margin. It follows that $\Gamma(z) + \theta \cdot\Gamma(1-z)$ and $\Gamma(z)$ have the same number of zeros minus the number of poles in $C_n$. Since $\Gamma(z)$ has no zeros and poles in $C_n$, it follows that $\Gamma(z) + \theta\cdot\Gamma(1-z)$ has the same number of zeros and poles. It has exactly one pole, and thus exactly one zero. If $\theta = \pm 1$ (and so in particular is real), by the Schwarz reflection principle, this zero is forced to be real. By symmetry, the same argument applies in the region $z = s + i t$ with $|t| \le 4$ and $s \le -15$. Combined with the above argument, this reduces the claim to $z = s + i t$ with $|s| \le 15$ and $|t| \le 4$ where the claim can be checked directly.
Hence all the zeros are either in $\mathbf{R}$, or lie on the line $1/2 + i \mathbf{R}$.
EDIT To clarify, I didn't actually check that there were no ``exceptional'' zeros in the box $\pm 15 \pm 4 I$, since I presumed that the original poster had done so. If $F(z) = \Gamma(z) - \Gamma(1-z)$, then computing the integral $$\frac{1}{2 \pi i} \oint \frac{F'(z)}{F(z)} dz$$ around that box, one obtains (numerically, and thus exactly) $1$. There are (assuming the OP at least computed the critical line zeros correctly) $2$ zeros in that range on the critical line. Along the real line in that range, there are $30$ poles and $25$ zeros. This means that there must be $1 + 30 - 25 = 6$ unaccounted for zeros. For such a zero $\rho$ off the line, by symmetry one also has $\overline{\rho}$, $1 - \rho$ and $1 - \overline{\rho}$ as zeros. Hence there must be either $1$ or $3$ pairs of zeros on the critical line, and either $1$ or $0$ quadruples of roots off the line. Varying the parameters of the integral, one can confirm there is a zero with $\rho \sim 2.7 + 0.3 i$, which is one of the four conjugates of the root found by joro. A similar argument applies for $\Gamma(z)+\Gamma(1-z)$. Hence:
Any zero of $\Gamma(z) - \Gamma(1-z)$ is either in $\mathbf{R}$, on the line $1/2 + i \mathbf{R}$, or is one of the four exceptional zeros $\{\rho,1-\rho,\overline{\rho},1-\overline{\rho}\}$. A similar calculation implies the same for $\Gamma(z) + \Gamma(1-z)$, except now with an exceptional set $\{\mu,1-\mu,\overline{\mu},1-\overline{\mu}\}$.
In the first part, we show that there are no zeros for $z = s + i t$ with $|t| \ge 4$ .
Let $\psi(z):= \Gamma'(z)/\Gamma(z)$ be the digamma function. If $z = s + i t$, then $$\frac{d}{ds} |\Gamma(z)|^2 = \frac{d}{ds} \Gamma(z) \Gamma(\overline{z}) = |\Gamma(z)|^2 \left(\psi(z) + \psi(\overline{z})\right).$$ (Both $\Gamma(z)$ and $\psi(z)$ are real for real $z$, and so satisfy the Schwartz reflection principle.) The product formula for the Gamma function implies that there is an identity $$\psi(z) = - \ \gamma + \sum_{n=1}^{\infty} \left(\frac{1}{n} - \frac{1}{z + n} \right) = 1 - \gamma + \sum_{n=1}^{\infty} \left(\frac{1}{n + 1} - \frac{1}{z + n} \right),$$ and hence $$\psi(z) + \psi(\overline{z}) = 2(1 - \gamma) + \sum_{n=1}^{\infty} \left(\frac{2}{n + 1} - \frac{1}{z + n} - \frac{1}{\overline{z} + n} \right).$$ Suppose that $z = s + i t$, and that $s \in [0,1]$. Then $$ \frac{2}{n + 1} - \frac{1}{s + i t + n} - \frac{1}{s - i t + n} = \frac{2(s^2 + t^2 + n s - s - n)}{(1+n)(n^2 + 2 n s + s^2 + t^2)} \ge \frac{-2}{(n^2 + t^2)}.$$ (The last inequality comes from ignoring all the positive terms in the numerator, and then setting $s = 0$ in the denominator.) It follows that $$\psi(z) + \psi(\overline{z}) \ge 2(1 - \gamma) - \sum_{n=1}^{\infty} \frac{2}{n^2 + t^2},$$ which is positive for $t$ big enough, e.g. $|t| \ge 4$. On the other hand, $$\psi(z + 1) + \psi(\overline{z} + 1) = \psi(z) + \psi(\overline{z}) + \frac{1}{z} + \frac{1}{\overline{z}} = \psi(z) + \psi(\overline{z}) + \frac{2s}{|z|^2}.$$ In particular, if $\psi(z) + \psi(\overline{z})$ is positive for $s \in [0,1]$ for some particular $t$, it is positive for all $s$ and that particular $t$. It follows that, if $|t| > 4$, that $|\Gamma(s + it)|^2$ is increasing as a function of $s$. In particular, if $|t| > 4$, then any equality $$|\Gamma(s + i t)| = |\Gamma(1 - (s + i t))| = |\Gamma(1 - s + i t)|$$ implies that $s = 1/2$.
The second part is a continuation of the argument above, which completes the argument. (merged from a different answer.)
Let $C_n$ denote the square with vertices $[n \pm 1/2, \pm 4 I]$ for a positive integer $n$. We have the following inequalities for $z \in C_n$ and $n \ge 15$: $$|\sin(\pi z)| \ge 1, \quad z \in C_n.$$ $$|\Gamma(z)| \ge \frac{1}{2} \Gamma(n - 1/2),$$ $$|\Gamma(1-z)| \le \frac{\pi}{\Gamma(n - 1/2)} \le 1,$$ $$|\psi(1-z)|, |\psi(z)| \le 2 \log(n), $$
The first is easy, the second follows from Stirling's formula (this requires $n$ to be big enough, and also requires $z$ to have imaginary part at most $4$), the third follows from the previous two by the reflection formula for $\Gamma(z)$, the last follows by induction and by the formula $\psi(z+1) = \psi(z) + 1/z$. It follows that $$\left| \frac{1}{2 \pi i} \oint_{C_n} \frac{\Gamma'(z)}{\Gamma(z)} - \frac{d/dz (\Gamma(z) + \theta \cdot \Gamma(1-z))}{\Gamma(z) + \theta\cdot \Gamma(1-z)} \right|$$ $$= \left| \frac{1}{2 \pi i} \oint_{C_n} \frac{\theta \Gamma(1-z) (\psi(1-z) + \psi(z))} {\Gamma(z) + \theta \cdot \Gamma(1-z)} \right|$$ $$ \le \frac{8 |\theta| \cdot \log(n) \pi}{2 \pi \cdot \Gamma(n - 1/2)} \oint_{C_n} \frac{1} {|\Gamma(z) + \theta \cdot \Gamma(1-z)|}$$ $$ \le \frac{8 |\theta| \cdot \log(n) \pi}{2 \pi \cdot \Gamma(n - 1/2)} \cdot \frac{1}{1/2 \Gamma(n - 1/2) + 1} \ll 1,$$ where $\theta = \pm 1$ (or anything small) and $n \ge 15$, where the final inequality holds by a huuuge margin. It follows that $\Gamma(z) + \theta \cdot\Gamma(1-z)$ and $\Gamma(z)$ have the same number of zeros minus the number of poles in $C_n$. Since $\Gamma(z)$ has no zeros and poles in $C_n$, it follows that $\Gamma(z) + \theta\cdot\Gamma(1-z)$ has the same number of zeros and poles. It has exactly one pole, and thus exactly one zero. If $\theta = \pm 1$ (and so in particular is real), by the Schwarz reflection principle, this zero is forced to be real. By symmetry, the same argument applies in the region $z = s + i t$ with $|t| \le 4$ and $s \le -15$. Combined with the above argument, this reduces the claim to $z = s + i t$ with $|s| \le 15$ and $|t| \le 4$ where the claim can be checked directly.
Hence all the zeros outside the box $z = s + it$ with $|t| \le 4$ and $|s| \le 15$ are either in $\mathbf{R}$, or lie on the line $1/2 + i \mathbf{R}$.
EDIT To clarify, I didn't actually check that there were no ``exceptional'' zeros in the box $\pm 15 \pm 4 I$, since I presumed that the original poster had done so. If $F(z) = \Gamma(z) - \Gamma(1-z)$, then computing the integral $$\frac{1}{2 \pi i} \oint \frac{F'(z)}{F(z)} dz$$ around that box, one obtains (numerically, and thus exactly) $1$. There are (assuming the OP at least computed the critical line zeros correctly) $2$ zeros in that range on the critical line. Along the real line in that range, there are $30$ poles and $25$ zeros. This means that there must be $1 + 30 - 25 = 6$ unaccounted for zeros. For such a zero $\rho$ off the line, by symmetry one also has $\overline{\rho}$, $1 - \rho$ and $1 - \overline{\rho}$ as zeros. Hence there must be either $1$ or $3$ pairs of zeros on the critical line, and either $1$ or $0$ quadruples of roots off the line. Varying the parameters of the integral, one can confirm there is a zero with $\rho \sim 2.7 + 0.3 i$, which is one of the four conjugates of the root found by joro. A similar argument applies for $\Gamma(z)+\Gamma(1-z)$. Hence:
Any zero of $\Gamma(z) - \Gamma(1-z)$ is either in $\mathbf{R}$, on the line $1/2 + i \mathbf{R}$, or is one of the four exceptional zeros $\{\rho,1-\rho,\overline{\rho},1-\overline{\rho}\}$. A similar calculation implies the same for $\Gamma(z) + \Gamma(1-z)$, except now with an exceptional set $\{\mu,1-\mu,\overline{\mu},1-\overline{\mu}\}$.
I would like to expand on Guild of Pepperers's answer by noting that the zeros are essentially uniformly spaced and may easily be approximated to a high degree of accuracy. Using Stirling approximation, I obtained the formula $$ \Gamma\left(\frac12+it\right) = \sqrt{\frac{2\pi}{1+e^{-2\pi|t|}}}\exp\left(-\frac\pi2|t|+i(t\log|t|-t+\varepsilon(t))\right), $$ valid for real $t$, where the error $\varepsilon(t)$ is an odd, bounded, real-valued function asymptotically equal to $\frac{1}{24t}$. (Indeed, $\varepsilon(t)$ has asymptotic and convergent expansions coming from the asymptotic and convergent versions of Stirling approximation, respectively.) We then have, for $s = \frac12+it$ on the critical line, $$ \Gamma(s)+\Gamma(1-s) = 2\sqrt{\frac{2\pi}{1+e^{-2\pi|t|}}}e^{-\frac\pi2|t|}\cos\left(t\log|t|-t+\varepsilon(t)\right), $$ $$ \Gamma(s)-\Gamma(1-s) = 2\sqrt{\frac{2\pi}{1+e^{-2\pi|t|}}}e^{-\frac\pi2|t|}\sin\left(t\log|t|-t+\varepsilon(t)\right). $$ One may show by means fair or foul that $t\log|t|-t+\varepsilon(t)$ is monotonically increasing for $|t|\geq1.05$, is bounded between $-0.96$ and $0.96$ for $|t|<1.05$, and is only zero when t = 0. Therefore, the zeros of $\Gamma(s)+\Gamma(1-s)$ on the critical line occur, with multiplicity one, very near those $t$ for which $t\log|t|-t$ is an odd integer multiple of $\frac{\pi}{2}$, and similarly for $\Gamma(s)-\Gamma(1-s)$ and the even integer multiples of $\frac{\pi}{2}$.
It's interesting that the number of zeros up to a given height $T$ is of the same order of magnitude, $T \log(T)$, as for the Riemann zeta function, but that these zeros have (essentially) uniform spacings rather than GUE spacings.