I come up with an elementary proof.
Let me first sketch the strategy. Given two real numbers $a<b$, a function $f\colon[a,b]\to\mathbb R$, and an interval $I\subseteq[a,b]$, we denote by $\DeclareMathOperator\osc{osc}\osc(f,I)=\sup f(I)-\inf f(I)$ the oscillation of $f$ on $I$.
We start with an elementary characterization of Riemann integrability:
Theorem (du Bois-Reymond) A function $f\colon[a,b]\to\mathbb R$ is Riemann integrable if and only if it is bounded, and for each $\epsilon,\delta>0$, there exists a partition $P=(a=x_0\le x_1\le\dots\le x_n=b)$ (allowing endpoints equal does not affect anything) of $[a,b]$ such that the total length of the intervals on which the oscillation of the function $\ge\epsilon$ is less than $\delta$, that is to say, $\sum_{\osc(f,[x_i,x_{i+1}])\le\epsilon}(x_{i+1}-x_i)<\delta$.
As indicated in OP, we denote by $C(f,P)$ the Cauchy sum of a function $f$ associated to a partition $P$, i.e. $\sum_{i=0}^{n-1}f(x_i)(x_{i+1}-x_i)$. We prove the contraposition: suppose that $f$ is not Riemann integrable, then by the du Bois-Reymond's characterization, we can pick up $\epsilon,\delta>0$ such that for each partition $P$ of $[a,b]$, the total length of the intervals on which the oscillation of the function $\ge\epsilon$ is at least $\delta$. Then for each partition $P$ of $[a,b]$, we construct two refinements $P',P''$ of $P$ such that the difference of Cauchy sums $\lvert C(f,P')-C(f,P'')\rvert$ is at least a constant which only depends on $a,b,f,\epsilon,\delta$ but not on the choice of $P$, then $f$ is not Cauchy-integrable.
Observations of the simplest case
To do this, we start with the simplest case: consider the trivial partition $P=(a=x_0<x_1=b)$. By assumption, $\osc(f,[a,b])\ge\epsilon$ and that $b-a\ge\delta$. It is natural to pick $u,v\in[a,b]$ such that $f(u)\approx\inf f([a,b])$ and $f(v)\approx\sup f([a,b])$. If $u<v$, then we consider two refinements $S=(a\le u\le b)$ and $T=(a\le u<v\le b)$ of $P$. Then $C(f,T)-C(f,S)=(b-v)(f(v)-f(u))$. Similarly, if $u>v$, we consider two refinements $S=(a\le v\le b)$ and $T=(a\le v<u\le b)$ of $P$. Then $C(f,S)-C(f,T)=(b-u)(f(v)-f(u))$. In both case, we can choose two refinements so that the difference of Cauchy sums is $(b-\max(u,v))(f(v)-f(u))\approx(b-\max(u,v))\osc(f,[a,b])\ge(b-\max(u,v))\epsilon$.
It would be a good choice if $b-\max(u,v)$ is large enough. In fact, if we modify the choice a bit, namely taking $f(u)\approx\inf f([a,(a+b)/2])$ and $f(v)\approx\sup f([a,(a+b)/2])$, then the same argument gives us two refinements so that the difference of Cauchy sums is
\[
(b-\max(u,v))(f(v)-f(u))\ge\frac{b-a}2(f(v)-f(u))\ge\frac\delta2(f(v)-f(u))\approx\frac\delta2\osc\left(f,\left[a,\frac{a+b}2\right]\right)
\]
If $\osc(f,[a,(a+b)/2])$ is large enough, that is to say, $\ge\epsilon$, then this difference is approximately no less than $2^{-1}\delta\epsilon$, we are done. However, this is of course not guarenteed.
The key idea is to bisect $[a,b]$ to obtain a new partition $(a<(a+b)/2<b)$. If $\osc(f,[a,(a+b)/2])\ge\epsilon$, then we are done as above. Otherwise, the du Bois-Reymond's characterization forces $b-(a+b)/2\ge\delta$. We replace $[a,b]$ by $[(a+b)/2,b]$ and continue the argument as above. If, for example, $\osc(f,[(a+b)/2,(a+3b)/4]<\epsilon$, we again bisect $[(a+b)/2,b]$ and restrict to the sub-interval to proceed. This process eventually terminates because the lengths of the intervals we get $\ge\delta$, so the process of bisection cannot be interminable.
In the general case, we also repeatedly bisect intervals to refine a partition and then pick out intervals, on the left half of which the oscillation of the function is at least $\epsilon$. The total length of the intervals we pick out is at least a constant times $\delta$. Then in each of these intervals, we pick out two points approximating infinum and supremum respectively. There is another techical point here: there are two cases, either inf lives to the left of sup, or inf lives to the right of sup. We pick out those intervals such that inf always lives to the left of sup, or always to the right of sup, so that the accumulations of differences on those intervals do not "cancel". By pidgeon-hole principle, either the first choice or the second choice will work. See the following sketch for details.
Sketch of details of the general proof
For an (closed) interval $I$, denote by $L(I)$ the left half of $I$, i.e. the closed interval from the left endpoint to the midpoint of $I$.
Denote by $\mathcal I(P)$ the intervals (cut by a partition $P$) on which the oscillation of the function $f$ is at least $\epsilon$, i.e. $\osc(f,I)\ge\epsilon$, and denote by $\mathcal J(P)\subseteq\mathcal I(P)$ the intervals $I\in\mathcal I(P)$ such that $\osc(f,L(I))\ge\epsilon$.
For each partition $P$, set $P_0=P$. We recursively construct a (finite) sequence of $(P_n)$: by the du Bois-Reymond characterization and the assumption that $f$ is not Riemann-integrable, the total length of intervals in $\mathcal I(P_n)$ is no less than $\delta$, i.e. $l_n:=\sum_{I\in\mathcal I(P_n)}m(I)\ge\delta$, where $m(I)$ is the length (i.e. the measure) of $I$. If $g_n:=\sum_{I\in\mathcal J(P_n)}m(I)\ge\delta/2$, then we terminate. Otherwise, we insert all midpoints of intervals $I\in\mathcal I(P_n)\setminus\mathcal J(P_n)$ into $P_n$, get a new partition $P_{n+1}$. Then we proceed from $n$ to $n+1$. Note that on one hand, $l_{n+1}:=\sum_{I\in\mathcal I(P_n)}m(I)\ge\delta$; on the other hand, by construction, we have bisected each $I\in\mathcal I(P_n)\setminus\mathcal J(P_n)$, the oscillation of the left half of which is less than $\epsilon$, $l_n-l_{n+1}=(l_n-g_n)/2$. Since by assumption, $l_n\ge\delta$ and $g_n<\delta/2$ (otherwise we terminate), we have $l_n-l_{n+1}>\delta/4$. We conclude that this procedure terminates eventually. We denote by $Q$ the $P_n$ we get when we terminate.
By construction, $\sum_{I\in\mathcal J(Q)}m(I)\ge\delta/2$. For each $I\in\mathcal J(Q)$, we pick two points $s_I,t_I\in L(I)$ such that $f(s_I)\approx\inf f(L(I))$ and $f(t_I)\approx\sup f(L(I))$. Since $I\in\mathcal J(Q)$, we have $f(t_I)-f(s_I)\approx\sup f(L(I))-\inf f(L(I))\ge\epsilon$. For each $I\in\mathcal J(Q)$, either $s_I<t_I$ or $s_I>t_I$. Denote by $\mathcal J_1\subset\mathcal J(Q)$ the subset of intervals $I\in\mathcal J(Q)$ such that $s_I<t_I$, and by $\mathcal J_2\subset\mathcal J(Q)$ the subset of intervals $I\in\mathcal J(Q)$ such that $s_I>t_I$. We note that $\sum_{I\in\mathcal J_1}m(I)+\sum_{I\in\mathcal J_2}m(I)=\sum_{I\in\mathcal J(Q)}m(I)\ge\delta/2$. Hence, the total length of intervals in either $\mathcal J_j$ is at least $\delta/4$, $j=1$ or $j=2$. We insert $\min(s_I,t_I)$ for each $I\in\mathcal J_j$ into $Q$, obtaining $Q'$ as a refinement of $Q$. We insert both $s_I$ and $t_I$ into $Q$, obtaining $Q''$ as a refinement of $Q'$. Now let's calculate $\lvert C(Q',f)-C(Q'',f)\rvert$. As in the simplest case of $P=(a<b)$, we deduce that
\[
\lvert C(Q',f)-C(Q'',f)\rvert\ge\frac{\sum_{I\in\mathcal J_j}m(I)}2\epsilon\ge\frac{\delta\epsilon}8
\]
where the first inequality is deduced from the fact that $f(t_I)-f(s_I)$ is (approximately) no less than $\epsilon$ for $I\in\mathcal J_j\subseteq\mathcal J(Q)$. Q.E.D.
