It seems that it´s fixed by adding this line in my training code:
hidden = tuple([each.data for each in hidden])

The resulting training code is now:
for epoch in range(1, n_epochs):
    hidden = lstm.init_hidden(batch_size)

    #now looping through data
    for batch_i, (features, targets) in enumerate(trainloader, 1):
        #watch out for incomplete batches
        print(batch_i)
        n_batches = len(trainloader.dataset) // batch_size
        if (batch_i > n_batches):
            break
        print(features)
        print(features.shape)
        print(targets)
        print(targets.shape)
        #defining a forward pass
        
        hidden = tuple([each.data for each in hidden])
        
        optimizer.zero_grad()
        output, hidden = lstm(features, hidden)
        loss = criterion(output, targets)
        loss.backward(retain_graph=True)
        optimizer.step()