Thanks @Amandeep!
I had a look, and this is to be expected. It is because PennyLane and Qiskit use different conventions for the underlying statevector ordering.
For more details, see here:
In this particular case, if you want to get the Qiskit results to match the PennyLane results, you can simply reshape and transpose the statevector:
def f_statevector(theta):
qc = construct_circuit(theta, 0)
statevector = backend_statevector.run(qc).result().get_statevector(qc)
statevector = statevector.reshape([2] * qubits).T.flatten()