The answers given to the question whether all zeros in the critical strip of $\zeta(s)\pm\zeta(1-s)$ lie on the critical line, suggest that this can indeed be proven, however only for those zeros where $s \ne \rho$ (to be more precise; those zeros occur when $\chi(s)=2^s \pi^{s-1} \sin(\pi s/2) \phantom. \Gamma(1-s) = \pm 1$).

Now assume $s \in \mathbb{C}$, $\Re(s) \ge 0$ and take the known expression:

$$\zeta(s) = \dfrac{s}{s-1} - \frac12+s \int_1^\infty \frac{1/2-\{x\}}{x^{s+1}}\,\mathrm{d}x$$

and substitute the fractional part of $\{x\}$ by a closed form (derived here):

$$\displaystyle \{x\} = x - \lfloor x \rfloor = \frac12 + \frac{i}{2 \pi} \ln \left(-\mathrm{e}^{-2 \pi i x} \right)$$

which gives:

$$\displaystyle \zeta(s) = \dfrac{s}{s-1} - \frac12-\frac{s i}{2 \pi} \int_1^\infty \frac{\ln \left(-\mathrm{e}^{-2\pi i x}\right)\,}{x^{s+1}}\mathrm{d}x $$

From the discussion in the comments section below, it has become clear that the various CAS-packages give different outcomes when evaluating this integral. When the integral is finite from 1 to $N$, both Maple and Mathematica give the correct outcome, but Sage seems to struggle. The difference probably can be explained from CAS picking the correct (principal?) branch of the multi-valued $\ln(-e)$ element in the integral. However, despite 2 CAS results continuously improving in accuracy with increasing $N$, in all CAS the integral is yielding a very wrong outcome at $\infty$. This is not only the case for the $\ln(-e)$ integral, but also for the integral with $\{x\}$ (which is a proven formula for $\zeta(s)$). I now wonder if this has something to do with how the various CAS evaluate the fractional part at infinity. In any case, CAS are not going to give us the answer and some real pen and paper math is required to assess what exactly happens at $\infty$. Any thoughts are welcome.

Since I used a finite integral in Maple to test my conjecture below (which still holds and no other zeros in the strip have been found so far), I decided to be more precise in the OP and replaced $\infty$ by an as large as you like $N$ in the integral below.

Isolate the integral part,

$$I(s) =\frac{s i}{2 \pi} \int_1^N \frac{\ln \left(-\mathrm{e}^{-2\pi ix}\right)\,}{x^{s+1}}\mathrm{d}x $$

and I like to conjecture, that in the critical strip, all zeros of:

$$I(s) \pm I(1-s)$$

are on the critical line $\Re(s)=\frac12$, however now with the certainty that when $s= \rho$ then $I(s) \ne 0$.

Via the reflective relation $\zeta(s) = \chi(s)\phantom . \zeta(1-s)$, this can be simplified into the following relation:

$$I(s) \pm I(1-s) =0 \text{ when } \displaystyle \chi(s)= \frac{\frac{s}{1-s} + \frac12 \pm -I(1-s)}{\frac{1-s}{s} + \frac12 + \hspace{3 mm} I(1-s)}$$

Appreciate any thoughts on possible approaches to proof this conjecture, e.g. based on the symmetry between the two integrals $I(s)$ and $I(1-s)$ or the symmetry around $\chi(s)$.

Thanks.