Arduino Nicla Sense ME IMU orientation sensing issues

I'm currently using an Arduino Nicla Sense ME to gather orientation data from the internal IMU and send it to a PC for further processing. However, I've run into a few strange issues with the orientation data that I get from the board.

When requesting data using the SensorOrientation class with the full rotation vector, I seem to get some non-linearities in the sensor output. After I rotate the nicla board 90 degrees (in a single axis), the board reads around a 120 degree rotation. However, after another 90 degree rotation it reads the expected value of 180 degrees. Due to this nonlinearity I don't believe this issue is caused by an incorrect scaling factor (I have the scaling factor set to 360/2^-15 as suggested in the datasheet), as linear scaling wouldn't cause an issue like this.

I have also tested retrieving quaternion outputs from the board, and although these do scale correctly, I've run into other issues here too. 2 of the axes seem to work as expected, however, when I rotate the board purely in the 3rd axis (when both other axes are at 0 degrees), the resulting sensor readings imply movement in all 3 axes.

I thought this may be a magnetometer calibration issue, however it is still present when using the game rotation vector (which doesn't use the magnetometer). I have also tested 2 boards and found the same issues present on both.

Code for reading euler angles:

#include "Arduino_BHY2.h"
#include "ArduinoBLE.h" // Required for Arduino_BHY2 to compile

// Get the full rotation vector (using all 9 axes of the IMU), as euler angles
SensorOrientation orientation(SENSOR_ID_RV);

void print_orientation()
{
    Serial.print(orientation.pitch());
    Serial.print(",");
    Serial.print(orientation.roll());
    Serial.print(",");
    Serial.println(orientation.heading());
}

// Scaling to convert from fixed point 16 bit value, to an angle in degrees
// See: https://www.bosch-sensortec.com/media/boschsensortec/downloads/datasheets/bst-bhi260ap-ds000.pdf
// Page 122, Table 990
constexpr float euler_scaling_factor = 360.0f/((float) (1 << 15));

constexpr float sample_rate_hz = 800;
constexpr unsigned int latency_ms = 0;

void setup()
{
    Serial.begin(1e6);

    BHY2.begin();

    orientation.configure(sample_rate_hz, latency_ms);
    orientation.setFactor(euler_scaling_factor);
}

unsigned long last_loop_time_us = 0;
constexpr unsigned int delay_val_us = 1e6/sample_rate_hz;

void loop()
{
    BHY2.update();
    if ((micros() - last_loop_time_us) > delay_val_us || (micros() - last_loop_time_us) < 0) { // x < 0 case is checking for overflow on micros()
      print_orientation();    
      last_loop_time_us = micros();
    }
}

Code for reading quaternions:

#include "Arduino_BHY2.h"
#include "ArduinoBLE.h" // Required for Arduino_BHY2 to compile

// Get the full rotation vector (using all 9 axes of the IMU), as quaternion angles
SensorQuaternion quaternion(SENSOR_ID_RV);

void print_orientation()
{
    Serial.print(quaternion.x());
    Serial.print(",");
    Serial.print(quaternion.y());
    Serial.print(",");
    Serial.print(quaternion.z());
    Serial.print(",");
    Serial.println(quaternion.w());
}

constexpr float sample_rate_hz = 800;
constexpr unsigned int latency_ms = 0;

void setup()
{
    Serial.begin(1e6);

    BHY2.begin();

    quaternion.configure(sample_rate_hz, latency_ms);
}

unsigned long last_loop_time_us = 0;
constexpr unsigned int delay_val_us = 1e6/sample_rate_hz;

void loop()
{
    BHY2.update();
    if ((micros() - last_loop_time_us) > delay_val_us || (micros() - last_loop_time_us) < 0) { // x < 0 case is checking for overflow on micros()
      print_orientation();    
      last_loop_time_us = micros();
    }
}

Any advice or thoughts would be appreciated.